Helicity Modulus

cmake/Sources.cmake

@@ -63,6 +63,7 @@ set(JAMS_SOURCES_CXX
         monitors/topological_charge_finite_diff.cc
         monitors/topological_charge_geometrical_def.cc
         monitors/hdf5.cc
+        monitors/helicity_modulus.cc
         monitors/magnetisation.cc
         monitors/magnetisation_layers.cc
         monitors/magnetisation_rate.cc

src/jams/core/hamiltonian.cc

@@ -95,6 +95,8 @@ Hamiltonian * Hamiltonian::create(const libconfig::Setting &settings, const unsi
 Hamiltonian::Hamiltonian(const libconfig::Setting &settings, const unsigned int size)
         : Base(settings),
           energy_(size),
+          helicity_internal_(size),
+          helicity_entropy_(size),
           field_(size, 3)
 {
 

src/jams/core/hamiltonian.h

@@ -42,6 +42,26 @@ class Hamiltonian : public Base {
     // calculate the energy difference of spin i in initial and final states
     virtual double calculate_energy_difference(int i, const Vec3 &spin_initial, const Vec3 &spin_final) = 0;
 
+    // calculate the internal energy difference between energy minima for spin i
+    virtual double calculate_internal_energy_difference(int i) {
+        return 0.0;
+    }
+
+    // calculate the total internal energy difference between energy minima
+    virtual double calculate_total_internal_energy_difference() {
+        return 0.0;
+    }
+
+    // calculate the entropy contribution to the helicity modulus for spin io
+    virtual double calculate_entropy(int i) {
+        return 0.0;
+    }
+
+    // calculate the total entropy contribution to the helicity modulus
+    virtual double calculate_total_entropy() {
+        return 0.0;
+    }
+
     inline double energy(const int i) const {
       assert(i < energy_.elements());
       return energy_(i);
@@ -76,6 +96,8 @@ class Hamiltonian : public Base {
 
     jams::MultiArray<double, 1> energy_; // energy of every spin for this Hamiltonian
     jams::MultiArray<double, 2> field_; // field at every spin for this Hamiltonianl
+    jams::MultiArray<double, 1> helicity_internal_; // internal energy difference between two minima for every spin
+    jams::MultiArray<double, 1> helicity_entropy_; // entropy contribution to the helicity modulus for each spin
 };
 
 

src/jams/core/monitor.cc

@@ -9,6 +9,7 @@
 #include "jams/monitors/topological_charge_finite_diff.h"
 #include "jams/monitors/topological_charge_geometrical_def.h"
 #include "jams/monitors/hdf5.h"
+#include "jams/monitors/helicity_modulus.h"
 #include "jams/monitors/magnetisation.h"
 #include "jams/monitors/magnetisation_layers.h"
 #include "jams/monitors/magnetisation_rate.h"
@@ -81,6 +82,7 @@ Monitor* Monitor::create(const Setting &settings) {
   DEFINED_MONITOR("topological-charge-finite-diff", TopologicalFiniteDiffChargeMonitor, settings);
   DEFINED_MONITOR("topological-charge-geometrical-def", TopologicalGeometricalDefMonitor, settings);
   DEFINED_MONITOR("hdf5", Hdf5Monitor, settings);
+  DEFINED_MONITOR("helicity-modulus", HelicityModulusMonitor, settings);
   DEFINED_MONITOR("magnetisation", MagnetisationMonitor, settings);
   DEFINED_MONITOR("magnetisation-layers", MagnetisationLayersMonitor, settings);
   DEFINED_MONITOR("magnetisation-rate", MagnetisationRateMonitor, settings);

src/jams/hamiltonian/cuda_uniaxial_anisotropy.cu

@@ -10,6 +10,8 @@
 CudaUniaxialHamiltonian::CudaUniaxialHamiltonian(const libconfig::Setting &settings, const unsigned int num_spins)
         : UniaxialHamiltonian(settings, num_spins)
 {
+    helicity_internal_.resize(globals::num_spins);
+    helicity_entropy_.resize(globals::num_spins);
 }
 
 double CudaUniaxialHamiltonian::calculate_total_energy() {
@@ -32,3 +34,33 @@ void CudaUniaxialHamiltonian::calculate_fields() {
             (globals::num_spins, power_, magnitude_.device_data(), axis_.device_data(), globals::s.device_data(), field_.device_data());
   DEBUG_CHECK_CUDA_ASYNC_STATUS;
 }
+
+void CudaUniaxialHamiltonian::calculate_internal_energy_differences() {
+    cuda_uniaxial_helicity_energy_kernel<<<(globals::num_spins+dev_blocksize_-1)/dev_blocksize_, dev_blocksize_, 0, dev_stream_.get()>>>
+            (globals::num_spins, power_, magnitude_.device_data(), axis_.device_data(), globals::s.device_data(), helicity_internal_.device_data());
+    DEBUG_CHECK_CUDA_ASYNC_STATUS;
+}
+
+void CudaUniaxialHamiltonian::calculate_entropies() {
+    cuda_uniaxial_entropy_kernel<<<(globals::num_spins+dev_blocksize_-1)/dev_blocksize_, dev_blocksize_, 0, dev_stream_.get()>>>
+            (globals::num_spins, power_, magnitude_.device_data(), axis_.device_data(), globals::s.device_data(), helicity_entropy_.device_data());
+    DEBUG_CHECK_CUDA_ASYNC_STATUS;
+}
+
+double CudaUniaxialHamiltonian::calculate_total_entropy() {
+    calculate_internal_energy_differences();
+    double TS_total = 0.0;
+    for (auto i = 0; i < energy_.size(); ++i) {
+        TS_total += helicity_entropy_(i);
+    }
+    return TS_total;
+}
+
+double CudaUniaxialHamiltonian::calculate_total_internal_energy_difference() {
+    calculate_entropies();
+    double dU_total = 0.0;
+    for (auto i = 0; i < energy_.size(); ++i) {
+        dU_total += helicity_internal_(i);
+    }
+    return dU_total;
+}

src/jams/hamiltonian/cuda_uniaxial_anisotropy.h

@@ -18,10 +18,17 @@ class CudaUniaxialHamiltonian : public UniaxialHamiltonian {
     double calculate_total_energy() override;
     void   calculate_energies() override;
     void   calculate_fields() override;
+    double calculate_total_internal_energy_difference() override;
+    double calculate_total_entropy() override;
+    void calculate_internal_energy_differences();
+    void calculate_entropies();
+
 private:
 
     CudaStream dev_stream_;
     unsigned int dev_blocksize_ = 64;
+    jams::MultiArray<double, 1> helicity_internal_;
+    jams::MultiArray<double, 1> helicity_entropy_;
 };
 
 #endif //JAMS_CUDA_UNIAXIAL_ANISOTROPY_H

src/jams/hamiltonian/cuda_uniaxial_anisotropy_kernel.cuh

@@ -24,3 +24,31 @@ __global__ void cuda_uniaxial_field_kernel(const int num_spins, const int power,
     dev_h[3 * idx + 2] = pre * a[2];
   }
 }
+
+__global__ void cuda_uniaxial_helicity_energy_kernel(const int num_spins, const int power,
+  const double * magnitude, const double * axis, const double * dev_s, double * dev_dU) {
+    const int idx = blockIdx.x*blockDim.x+threadIdx.x;
+    if (idx < num_spins) {
+        const double s[3] = {dev_s[3*idx], dev_s[3*idx+1], dev_s[3*idx+2]};
+        const double a[3] = {axis[3*idx], axis[3*idx+1], axis[3*idx+2]};
+
+        const double cross[3] = {cross_product_x(s, a), cross_product_y(s, a), cross_product_z(s, a)};
+        const double cross_sq = dot(cross, cross);
+        const double dot_prod = dot(s, a);
+
+        dev_dU[idx] = magnitude[idx] * power * (pow(dot_prod, power) - (power - 1)*cross_sq*pow(dot_prod, power-2) );
+    }
+}
+__global__ void cuda_uniaxial_entropy_kernel(const int num_spins, const int power,
+                                                     const double * magnitude, const double * axis, const double * dev_s, double * dev_TS) {
+    const int idx = blockIdx.x*blockDim.x+threadIdx.x;
+    if (idx < num_spins) {
+        const double s[3] = {dev_s[3*idx], dev_s[3*idx+1], dev_s[3*idx+2]};
+        const double a[3] = {axis[3*idx], axis[3*idx+1], axis[3*idx+2]};
+
+        const double cross[3] = {cross_product_x(s, a), cross_product_y(s, a), cross_product_z(s, a)};
+        const double cross_norm = sqrt(dot(cross, cross));
+
+        dev_TS[idx] = magnitude[idx] * power * pow(dot(s, a), power-1) * cross_norm;
+    }
+}

src/jams/hamiltonian/exchange.cc

@@ -150,4 +150,8 @@ ExchangeHamiltonian::ExchangeHamiltonian(const libconfig::Setting &settings, con
 
 const jams::InteractionList<Mat3,2> &ExchangeHamiltonian::neighbour_list() const {
   return neighbour_list_;
+}
+
+const double &ExchangeHamiltonian::get_input_conversion() const {
+    return input_energy_unit_conversion_;
 }

src/jams/hamiltonian/exchange.h

@@ -18,6 +18,8 @@ class ExchangeHamiltonian : public SparseInteractionHamiltonian {
 
     const jams::InteractionList<Mat3, 2> &neighbour_list() const;
 
+    const double &get_input_conversion() const;
+
 private:
     jams::InteractionList<Mat3, 2> neighbour_list_; // neighbour information
     double energy_cutoff_; // abs cutoff energy for interaction

src/jams/hamiltonian/uniaxial_anisotropy.cc

@@ -173,3 +173,60 @@ void UniaxialHamiltonian::calculate_fields() {
     }
   }
 }
+
+double UniaxialHamiltonian::calculate_internal_energy_difference(const int i) {
+    using namespace globals;
+    double dU = 0.0;
+
+    const auto dot = (axis_(i,0) * s(i,0) + axis_(i,1) * s(i,1) + axis_(i,2) * s(i,2));
+
+    const Vec3 cross = {s(i,1)*axis_(i,2) - s(i,2)*axis_(i,1), s(i,2)*axis_(i,0) - s(i,0)*axis_(i,2), s(i,0)*axis_(i,1) - s(i,1)*axis_(i,0)};
+    const auto mod_cross_sq = norm_sq(cross);
+
+    dU += power_*magnitude_(i)*( pow(dot, power_) - (power_ - 1) * mod_cross_sq*pow(dot, power_-2) );
+
+    return dU;
+}
+
+double UniaxialHamiltonian::calculate_total_internal_energy_difference() {
+    using namespace globals;
+    double dU_total = 0.0;
+
+    #if HAS_OMP
+    #pragma omp parallel for default(none) shared(num_spins) reduction(+:dU_total)
+    #endif
+    for (int i = 0; i < num_spins; ++i) {
+        dU_total += calculate_internal_energy_difference(i);
+    }
+    return dU_total;
+}
+
+double UniaxialHamiltonian::calculate_entropy(int i) {
+    using namespace globals;
+    double TS = 0.0;
+
+    const Vec3 spin = {s(i,0), s(i,1), s(i,2)};
+    const Vec3 axis = {axis_(i,0), axis_(i,1), axis_(i,2)};
+
+    const auto dot_product = dot(spin, axis);
+
+    const auto cross_product_norm = norm(cross(spin, axis));
+
+    TS += power_*magnitude_(i) * pow(dot_product, power_-1) * cross_product_norm;
+
+    return TS;
+}
+
+double UniaxialHamiltonian::calculate_total_entropy() {
+    using namespace globals;
+    double TS_total = 0.0;
+
+    #if HAS_OMP
+    #pragma omp parallel for default(none) shared(num_spins) reduction(+:TS_total)
+    #endif
+    for (int i = 0; i < num_spins; ++i) {
+        TS_total += calculate_entropy(i);
+    }
+
+     return TS_total;
+}

src/jams/hamiltonian/uniaxial_anisotropy.h

@@ -27,6 +27,14 @@ class UniaxialHamiltonian : public Hamiltonian {
 
     double calculate_energy_difference(int i, const Vec3 &spin_initial, const Vec3 &spin_final) override;
 
+    double calculate_internal_energy_difference(int i) override;
+
+    double calculate_total_internal_energy_difference() override;
+
+    double calculate_entropy(int i) override;
+
+    double calculate_total_entropy() override;
+
 private:
     int power_; // anisotropy power exponent
     jams::MultiArray<double, 2> axis_; // local uniaxial anisotropy axis