diff --git a/cmake/Sources.cmake b/cmake/Sources.cmake index ea3ea3d8..276dc98e 100644 --- a/cmake/Sources.cmake +++ b/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 diff --git a/src/jams/core/hamiltonian.cc b/src/jams/core/hamiltonian.cc index c2fa19ba..1f2b4b6f 100755 --- a/src/jams/core/hamiltonian.cc +++ b/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) { diff --git a/src/jams/core/hamiltonian.h b/src/jams/core/hamiltonian.h index f435e448..d87047f5 100644 --- a/src/jams/core/hamiltonian.h +++ b/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 energy_; // energy of every spin for this Hamiltonian jams::MultiArray field_; // field at every spin for this Hamiltonianl + jams::MultiArray helicity_internal_; // internal energy difference between two minima for every spin + jams::MultiArray helicity_entropy_; // entropy contribution to the helicity modulus for each spin }; diff --git a/src/jams/core/monitor.cc b/src/jams/core/monitor.cc index 2ac78672..cb2b062e 100644 --- a/src/jams/core/monitor.cc +++ b/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); diff --git a/src/jams/hamiltonian/cuda_uniaxial_anisotropy.cu b/src/jams/hamiltonian/cuda_uniaxial_anisotropy.cu index 86a84a73..86034569 100644 --- a/src/jams/hamiltonian/cuda_uniaxial_anisotropy.cu +++ b/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; +} diff --git a/src/jams/hamiltonian/cuda_uniaxial_anisotropy.h b/src/jams/hamiltonian/cuda_uniaxial_anisotropy.h index 02a07c6b..62e462da 100644 --- a/src/jams/hamiltonian/cuda_uniaxial_anisotropy.h +++ b/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 helicity_internal_; + jams::MultiArray helicity_entropy_; }; #endif //JAMS_CUDA_UNIAXIAL_ANISOTROPY_H diff --git a/src/jams/hamiltonian/cuda_uniaxial_anisotropy_kernel.cuh b/src/jams/hamiltonian/cuda_uniaxial_anisotropy_kernel.cuh index 89b34348..a790b1c8 100644 --- a/src/jams/hamiltonian/cuda_uniaxial_anisotropy_kernel.cuh +++ b/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; + } +} diff --git a/src/jams/hamiltonian/exchange.cc b/src/jams/hamiltonian/exchange.cc index e255b12f..bee07fad 100755 --- a/src/jams/hamiltonian/exchange.cc +++ b/src/jams/hamiltonian/exchange.cc @@ -150,4 +150,8 @@ ExchangeHamiltonian::ExchangeHamiltonian(const libconfig::Setting &settings, con const jams::InteractionList &ExchangeHamiltonian::neighbour_list() const { return neighbour_list_; +} + +const double &ExchangeHamiltonian::get_input_conversion() const { + return input_energy_unit_conversion_; } \ No newline at end of file diff --git a/src/jams/hamiltonian/exchange.h b/src/jams/hamiltonian/exchange.h index b45b71e5..a2d720d8 100644 --- a/src/jams/hamiltonian/exchange.h +++ b/src/jams/hamiltonian/exchange.h @@ -18,6 +18,8 @@ class ExchangeHamiltonian : public SparseInteractionHamiltonian { const jams::InteractionList &neighbour_list() const; + const double &get_input_conversion() const; + private: jams::InteractionList neighbour_list_; // neighbour information double energy_cutoff_; // abs cutoff energy for interaction diff --git a/src/jams/hamiltonian/uniaxial_anisotropy.cc b/src/jams/hamiltonian/uniaxial_anisotropy.cc index 5d6e317f..50ec8bbc 100755 --- a/src/jams/hamiltonian/uniaxial_anisotropy.cc +++ b/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; +} \ No newline at end of file diff --git a/src/jams/hamiltonian/uniaxial_anisotropy.h b/src/jams/hamiltonian/uniaxial_anisotropy.h index 32c2662b..e15d3e3f 100755 --- a/src/jams/hamiltonian/uniaxial_anisotropy.h +++ b/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 axis_; // local uniaxial anisotropy axis diff --git a/src/jams/monitors/helicity_modulus.cc b/src/jams/monitors/helicity_modulus.cc new file mode 100644 index 00000000..7f24a171 --- /dev/null +++ b/src/jams/monitors/helicity_modulus.cc @@ -0,0 +1,322 @@ +// +// Created by Sean Stansill [ll14s26s] on 26/09/2022. +// + +#include +#include +#include + +#include "jams/helpers/consts.h" +#include "jams/core/lattice.h" +#include "jams/core/physics.h" +#include "jams/core/solver.h" +#include "jams/core/globals.h" +#include "jams/core/hamiltonian.h" +#include "jams/helpers/output.h" +#include "jams/containers/sparse_matrix_builder.h" +#include "jams/hamiltonian/exchange.h" + +#include "helicity_modulus.h" + +using namespace std; + +HelicityModulusMonitor::HelicityModulusMonitor(const libconfig::Setting &settings) + : Monitor(settings), tsv_file(jams::output::full_path_filename("free_eng.tsv")) { + tsv_file.setf(std::ios::right); + tsv_file << tsv_header(); + + const auto& exchange_hamiltonian = find_hamiltonian(::solver->hamiltonians()); + + jams::SparseMatrix::Builder sparse_matrix_builder_JRR_xx(3 *globals::num_spins, 3 * globals::num_spins); + jams::SparseMatrix::Builder sparse_matrix_builder_JRR_xy(3 *globals::num_spins, 3 * globals::num_spins); + jams::SparseMatrix::Builder sparse_matrix_builder_JRR_xz(3 *globals::num_spins, 3 * globals::num_spins); + jams::SparseMatrix::Builder sparse_matrix_builder_JRR_yy(3 *globals::num_spins, 3 * globals::num_spins); + jams::SparseMatrix::Builder sparse_matrix_builder_JRR_yz(3 *globals::num_spins, 3 * globals::num_spins); + jams::SparseMatrix::Builder sparse_matrix_builder_JRR_zz(3 *globals::num_spins, 3 * globals::num_spins); + + jams::SparseMatrix::Builder sparse_matrix_builder_JR_x(3 *globals::num_spins, 3 * globals::num_spins); + jams::SparseMatrix::Builder sparse_matrix_builder_JR_y(3 *globals::num_spins, 3 * globals::num_spins); + jams::SparseMatrix::Builder sparse_matrix_builder_JR_z(3 *globals::num_spins, 3 * globals::num_spins); + + const auto& input_energy_conversion = exchange_hamiltonian.get_input_conversion(); + const auto& nbr_list = exchange_hamiltonian.neighbour_list(); + for (auto n = 0; n < nbr_list.size(); ++n) { + auto i = nbr_list[n].first[0]; + auto j = nbr_list[n].first[1]; + auto Jij = input_energy_conversion * nbr_list[n].second; + auto rij = lattice->displacement(i, j); + rij = abs(rij); + + Mat3 JijRij_x = rij[0] * Jij; + Mat3 JijRij_y = rij[1] * Jij; + Mat3 JijRij_z = rij[2] * Jij; + + Mat3 JijRijRij_xx = rij[0]*rij[0] * Jij; + Mat3 JijRijRij_xy = rij[0]*rij[1] * Jij; + Mat3 JijRijRij_xz = rij[0]*rij[2] * Jij; + Mat3 JijRijRij_yy = rij[1]*rij[1] * Jij; + Mat3 JijRijRij_yz = rij[1]*rij[2] * Jij; + Mat3 JijRijRij_zz = rij[2]*rij[2] * Jij; + + for (auto m = 0; m < 3; ++m) { + for (auto l = 0; l < 3; ++l) { + if (JijRijRij_xx[m][l] != 0.0) { + sparse_matrix_builder_JRR_xx.insert(3 * i + m, 3 * j + l, JijRijRij_xx[m][l]); + } + if (JijRijRij_xy[m][l] != 0.0) { + sparse_matrix_builder_JRR_xy.insert(3 * i + m, 3 * j + l, JijRijRij_xy[m][l]); + } + if (JijRijRij_xz[m][l] != 0.0) { + sparse_matrix_builder_JRR_xz.insert(3 * i + m, 3 * j + l, JijRijRij_xz[m][l]); + } + if (JijRijRij_yy[m][l] != 0.0) { + sparse_matrix_builder_JRR_yy.insert(3 * i + m, 3 * j + l, JijRijRij_yy[m][l]); + } + if (JijRijRij_yz[m][l] != 0.0) { + sparse_matrix_builder_JRR_yz.insert(3 * i + m, 3 * j + l, JijRijRij_yz[m][l]); + } + if (JijRijRij_zz[m][l] != 0.0) { + sparse_matrix_builder_JRR_zz.insert(3 * i + m, 3 * j + l, JijRijRij_zz[m][l]); + } + + + if (JijRij_x[m][l] != 0.0) { + sparse_matrix_builder_JR_x.insert(3 * i + m, 3 * j + l, JijRij_x[m][l]); + } + if (JijRij_y[m][l] != 0.0) { + sparse_matrix_builder_JR_y.insert(3 * i + m, 3 * j + l, JijRij_y[m][l]); + } + if (JijRij_z[m][l] != 0.0) { + sparse_matrix_builder_JR_z.insert(3 * i + m, 3 * j + l, JijRij_z[m][l]); + } + } + } + } + cout << " dipole sparse matrix builder memory " << sparse_matrix_builder_JRR_xx.memory() / kBytesToMegaBytes << "(MB)\n"; + cout << " building CSR matrix\n"; + interaction_JRR_xx_ = sparse_matrix_builder_JRR_xx + .set_format(jams::SparseMatrixFormat::CSR) + .build(); + interaction_JRR_xy_ = sparse_matrix_builder_JRR_xy + .set_format(jams::SparseMatrixFormat::CSR) + .build(); + interaction_JRR_xz_ = sparse_matrix_builder_JRR_xz + .set_format(jams::SparseMatrixFormat::CSR) + .build(); + interaction_JRR_yy_ = sparse_matrix_builder_JRR_yy + .set_format(jams::SparseMatrixFormat::CSR) + .build(); + interaction_JRR_yz_ = sparse_matrix_builder_JRR_yz + .set_format(jams::SparseMatrixFormat::CSR) + .build(); + interaction_JRR_zz_ = sparse_matrix_builder_JRR_zz + .set_format(jams::SparseMatrixFormat::CSR) + .build(); + + interaction_JR_x_ = sparse_matrix_builder_JR_x + .set_format(jams::SparseMatrixFormat::CSR) + .build(); + interaction_JR_y_ = sparse_matrix_builder_JR_y + .set_format(jams::SparseMatrixFormat::CSR) + .build(); + interaction_JR_z_ = sparse_matrix_builder_JR_z + .set_format(jams::SparseMatrixFormat::CSR) + .build(); + cout << " exchange sparse matrix memory (CSR): " << sparse_matrix_builder_JRR_xx.memory() / kBytesToMegaBytes << " (MB)\n"; + + + helicity_field_rxrx_.resize(globals::num_spins, 3); + helicity_field_rxry_.resize(globals::num_spins, 3); + helicity_field_rxrz_.resize(globals::num_spins, 3); + helicity_field_ryry_.resize(globals::num_spins, 3); + helicity_field_ryrz_.resize(globals::num_spins, 3); + helicity_field_rzrz_.resize(globals::num_spins, 3); + + entropy_field_x_.resize(globals::num_spins, 3); + entropy_field_y_.resize(globals::num_spins, 3); + entropy_field_z_.resize(globals::num_spins, 3); + +} + +void HelicityModulusMonitor::update(Solver * solver) { + using namespace globals; + + tsv_file.width(12); + + tsv_file << std::scientific << solver->time() << "\t"; + + for (auto &hamiltonian : solver->hamiltonians()) { + if (hamiltonian->name() == "exchange") { + + calculate_helicity_fields(); + calculate_entropy_fields(); + + Mat3 energy_difference = exchange_total_internal_energy_difference(); + Mat3 entropy = exchange_total_entropy(); + + for (auto i = 0; i < 3; ++i) { + for (auto j = 0; j < 3; ++j) { + tsv_file << std::scientific << std::setprecision(15) << energy_difference[i][j] << "\t"; + tsv_file << std::scientific << std::setprecision(15) << entropy[i][j] << "\t"; + } + } + } + + else { + auto energy_difference = hamiltonian->calculate_total_internal_energy_difference(); + auto entropy = hamiltonian->calculate_total_entropy(); + + tsv_file << std::scientific << std::setprecision(15) << energy_difference << "\t"; + tsv_file << std::scientific << std::setprecision(15) << entropy << "\t"; + } + } + + tsv_file << std::endl; +} + +std::string HelicityModulusMonitor::tsv_header() { + std::stringstream ss; + ss.width(12); + + ss << "time\t"; + for (auto &hamiltonian : solver->hamiltonians()) { + if (hamiltonian->name() == "exchange") { + string unit_vecs[3] = {"x", "y", "z"}; + for (auto i = 0; i < 3; ++i){ + + + for (auto j = 0; j < 3; ++j){ + ss << hamiltonian->name() << "_dU_" << unit_vecs[i] << unit_vecs[j] <<"_E_meV\t"; + ss << hamiltonian->name() << "_TS_" << unit_vecs[i] << unit_vecs[j] <<"_E_meV\t"; + } + } + } + + else { + ss << hamiltonian->name() << "_dU_E_meV\t"; + ss << hamiltonian->name() << "_TS_E_meV\t"; + } + } + + ss << std::endl; + + return ss.str(); +} + +void HelicityModulusMonitor::calculate_helicity_fields() { +#if HAS_CUDA + if (jams::instance().mode() == jams::Mode::GPU) { + interaction_JRR_xx_.multiply_gpu(globals::s, helicity_field_rxrx_, jams::instance().cusparse_handle(), cusparse_stream_.get()); + interaction_JRR_xy_.multiply_gpu(globals::s, helicity_field_rxry_, jams::instance().cusparse_handle(), cusparse_stream_.get()); + interaction_JRR_xz_.multiply_gpu(globals::s, helicity_field_rxrz_, jams::instance().cusparse_handle(), cusparse_stream_.get()); + + interaction_JRR_yy_.multiply_gpu(globals::s, helicity_field_ryry_, jams::instance().cusparse_handle(), cusparse_stream_.get()); + interaction_JRR_yz_.multiply_gpu(globals::s, helicity_field_ryrz_, jams::instance().cusparse_handle(), cusparse_stream_.get()); + interaction_JRR_zz_.multiply_gpu(globals::s, helicity_field_rzrz_, jams::instance().cusparse_handle(), cusparse_stream_.get()); + return; + } +#endif + + interaction_JRR_xx_.multiply(globals::s, helicity_field_rxrx_); + interaction_JRR_xy_.multiply(globals::s, helicity_field_rxry_); + interaction_JRR_xz_.multiply(globals::s, helicity_field_rxrz_); + + interaction_JRR_yy_.multiply(globals::s, helicity_field_ryry_); + interaction_JRR_yz_.multiply(globals::s, helicity_field_ryrz_); + interaction_JRR_zz_.multiply(globals::s, helicity_field_rzrz_); + +} + +void HelicityModulusMonitor::calculate_entropy_fields() { +#if HAS_CUDA + if (jams::instance().mode() == jams::Mode::GPU) { + interaction_JR_x_.multiply_gpu(globals::s, entropy_field_x_, jams::instance().cusparse_handle(), cusparse_stream_.get()); + interaction_JR_y_.multiply_gpu(globals::s, entropy_field_y_, jams::instance().cusparse_handle(), cusparse_stream_.get()); + interaction_JR_z_.multiply_gpu(globals::s, entropy_field_z_, jams::instance().cusparse_handle(), cusparse_stream_.get()); + return; + } +#endif + + interaction_JR_x_.multiply(globals::s, entropy_field_x_); + interaction_JR_y_.multiply(globals::s, entropy_field_y_); + interaction_JR_z_.multiply(globals::s, entropy_field_z_); +} + +Mat3 HelicityModulusMonitor::exchange_total_internal_energy_difference() { + using namespace globals; + Mat3 dU = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; + + for (auto i = 0; i < num_spins; ++i) { + + Vec3 s_i = {s(i,0), s(i,1), s(i,2)}; + + // TO DO: ONLY CALCULATE UPPER TRIANGLE OF THE STIFFNESS TENSOR + + // COULD ADD ACCURACY IF CALCULATE THE MEAN OF THE UPPER AND LOWER + // TRIANGLE + + Vec3 h_i_xx = {helicity_field_rxrx_(i,0), helicity_field_rxrx_(i, 1), helicity_field_rxrx_(i, 2)}; + Vec3 h_i_xy = {helicity_field_rxry_(i,0), helicity_field_rxry_(i, 1), helicity_field_rxry_(i, 2)}; + Vec3 h_i_xz = {helicity_field_rxrz_(i,0), helicity_field_rxrz_(i, 1), helicity_field_rxrz_(i, 2)}; + + Vec3 h_i_yx = h_i_xy; + Vec3 h_i_yy = {helicity_field_ryry_(i,0), helicity_field_ryry_(i, 1), helicity_field_ryry_(i, 2)}; + Vec3 h_i_yz = {helicity_field_ryrz_(i,0), helicity_field_ryrz_(i, 1), helicity_field_ryrz_(i, 2)}; + + Vec3 h_i_zx = h_i_xz; + Vec3 h_i_zy = h_i_yz; + Vec3 h_i_zz = {helicity_field_rzrz_(i,0), helicity_field_rzrz_(i, 1), helicity_field_rzrz_(i, 2)}; + + dU[0][0] += dot(s_i, h_i_xx); + dU[1][1] += dot(s_i, h_i_yy); + dU[2][2] += dot(s_i, h_i_zz); + + dU[0][1] += dot(s_i, h_i_xy); + dU[0][2] += dot(s_i, h_i_xz); + dU[1][2] += dot(s_i, h_i_yz); + + dU[1][0] += dot(s_i, h_i_yx); + dU[2][0] += dot(s_i, h_i_zx); + dU[2][1] += dot(s_i, h_i_zy); + + } + + return dU; +} + +Mat3 HelicityModulusMonitor::exchange_total_entropy() { + using namespace globals; + Mat3 TS = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; + + const Vec3 dot_vec = {1.0, 1.0, 1.0}; + + for (auto i = 0; i < globals::num_spins; ++i) { + + + const Vec3 s_i = {s(i,0), s(i,1), s(i,2)}; + + const Vec3 h_i_x = {entropy_field_x_(i,0), entropy_field_x_(i, 1), entropy_field_x_(i, 2)}; + const Vec3 h_i_y = {entropy_field_y_(i,0), entropy_field_y_(i, 1), entropy_field_y_(i, 2)}; + const Vec3 h_i_z = {entropy_field_z_(i,0), entropy_field_z_(i, 1), entropy_field_z_(i, 2)}; + + const double cross_x = dot(cross(s_i, h_i_x), dot_vec); + const double cross_y = dot(cross(s_i, h_i_y), dot_vec); + const double cross_z = dot(cross(s_i, h_i_z), dot_vec); + + TS[0][0] += cross_x * cross_x; + TS[1][1] += cross_y * cross_y; + TS[2][2] += cross_z * cross_z; + + TS[0][1] += abs(cross_x * cross_y); + TS[1][0] += abs(cross_x * cross_y); + + TS[0][2] += abs(cross_x * cross_z); + TS[2][0] += abs(cross_x * cross_z); + + TS[1][2] += abs(cross_y * cross_z); + TS[2][1] += abs(cross_y * cross_z); + + } + + return TS; +} diff --git a/src/jams/monitors/helicity_modulus.h b/src/jams/monitors/helicity_modulus.h new file mode 100644 index 00000000..32fbe530 --- /dev/null +++ b/src/jams/monitors/helicity_modulus.h @@ -0,0 +1,69 @@ +// +// Created by Sean Stansill [ll14s26s] on 26/09/2022. +// + +#ifndef JAMS_HELICITY_MODULUS_H +#define JAMS_HELICITY_MODULUS_H + +#if HAS_CUDA +#include +#include "jams/cuda/cuda_stream.h" +#include "jams/cuda/cuda_common.h" +#endif + +#include + +#include + +#include "jams/containers/sparse_matrix.h" +#include "jams/core/types.h" +#include "jams/core/solver.h" +#include "jams/core/monitor.h" +#include "jams/core/physics.h" + +class HelicityModulusMonitor : public Monitor { +public: + explicit HelicityModulusMonitor(const libconfig::Setting &settings); + + ~HelicityModulusMonitor() override = default; + + void update(Solver *solver) override; + void post_process() override {}; + +private: + std::ofstream tsv_file; + std::string tsv_header(); + Mat3 exchange_total_internal_energy_difference(); + Mat3 exchange_total_entropy(); + void calculate_helicity_fields(); + void calculate_entropy_fields(); + jams::SparseMatrix interaction_JRR_xx_; + jams::SparseMatrix interaction_JRR_xy_; + jams::SparseMatrix interaction_JRR_xz_; + jams::SparseMatrix interaction_JRR_yy_; + jams::SparseMatrix interaction_JRR_yz_; + jams::SparseMatrix interaction_JRR_zz_; + + jams::SparseMatrix interaction_JR_x_; + jams::SparseMatrix interaction_JR_y_; + jams::SparseMatrix interaction_JR_z_; + + jams::MultiArray helicity_field_rxrx_; + jams::MultiArray helicity_field_rxry_; + jams::MultiArray helicity_field_rxrz_; + + jams::MultiArray helicity_field_ryry_; + jams::MultiArray helicity_field_ryrz_; + jams::MultiArray helicity_field_rzrz_; + + jams::MultiArray entropy_field_x_; + jams::MultiArray entropy_field_y_; + jams::MultiArray entropy_field_z_; + + #if HAS_CUDA + CudaStream cusparse_stream_; // cuda stream to run in + #endif + +}; + +#endif //JAMS_HELICITY_MODULUS_H \ No newline at end of file