From 182bcef09e954a29f765951f2d49927645ba5fc1 Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Thu, 6 Feb 2020 14:36:13 +0000 Subject: [PATCH 01/80] Fixed some bugs & added some comments. --- src/jams/helpers/neutrons.cc | 3 +- .../monitors/neutron_scattering_no_lattice.cc | 37 +++++++++++++------ .../monitors/neutron_scattering_no_lattice.h | 8 ++-- 3 files changed, 31 insertions(+), 17 deletions(-) diff --git a/src/jams/helpers/neutrons.cc b/src/jams/helpers/neutrons.cc index da1a121e..6fe35bf3 100644 --- a/src/jams/helpers/neutrons.cc +++ b/src/jams/helpers/neutrons.cc @@ -9,10 +9,11 @@ namespace jams { // Calculates the approximate neutron form factor at |q| // Approximation and constants from International Tables for Crystallography: Vol. C (pp. 454–461). double form_factor(const Vec3 &q, const double &lattice_parameter, FormFactorG &g, FormFactorJ &j) { + //(*) How is "norm(q)/(4.0*kPi*lattice_parameter)" is related to s=sin(theta)/lambda in the ref.? auto s_sq = pow2(norm(q) / (4.0 * kPi * lattice_parameter)); auto ffq = 0.0; for (auto l : {0, 2, 4, 6}) { - double p = (l == 0) ? 1.0 : s_sq; + double p = (l == 0) ? 1.0 : s_sq; //if l=0 then insert 1.0 to p, otherwise s_sq ffq += g[l] * p * (j[l].A * exp(-j[l].a * s_sq) + j[l].B * exp(-j[l].b * s_sq) + j[l].C * exp(-j[l].c * s_sq) + j[l].D); } diff --git a/src/jams/monitors/neutron_scattering_no_lattice.cc b/src/jams/monitors/neutron_scattering_no_lattice.cc index 83c7b7ff..5628b0ec 100644 --- a/src/jams/monitors/neutron_scattering_no_lattice.cc +++ b/src/jams/monitors/neutron_scattering_no_lattice.cc @@ -32,12 +32,14 @@ NeutronScatteringNoLatticeMonitor::NeutronScatteringNoLatticeMonitor(const libco configure_polarizations(settings["polarizations"]); } + //insert values "length" and "overlap" from cfg if (settings.exists("periodogram")) { configure_periodogram(settings["periodogram"]); } periodogram_props_.sample_time = output_step_freq_ * solver->time_step(); + //zero: set all the elements as 0 zero(static_structure_factor_.resize(kspace_path_.size())); zero(kspace_spins_timeseries_.resize(periodogram_props_.length, kspace_path_.size())); zero(total_unpolarized_neutron_cross_section_.resize( @@ -72,19 +74,26 @@ void NeutronScatteringNoLatticeMonitor::update(Solver *solver) { } void NeutronScatteringNoLatticeMonitor::configure_kspace_vectors(const libconfig::Setting &settings) { - kvector_ = jams::config_optional(settings, "kvector", kvector_); + kmax_ = jams::config_required(settings, "kmax"); + kvector_ = jams::config_required(settings, "kvector"); + num_k_ = jams::config_required(settings, "num_k"); kspace_path_.resize(num_k_); for (auto i = 0; i < kspace_path_.size(); ++i) { - kspace_path_(i) = kvector_ * i * (kmax_ / num_k_); + //kspace_path_ stores the scattering vectors + kspace_path_(i) = kvector_ * i * (kmax_ / (num_k_-1));//kmax shoul be in the unit of [1/a] + cout << "i = " << i << ", kspace_path_(i) = " << kspace_path_(i) << endl; } rspace_displacement_.resize(globals::s.size(0)); for (auto i = 0; i < globals::s.size(0); ++i) { + //rspace_displacement_ is the position vector of i-th spin rspace_displacement_(i) = lattice->displacement({0,0,0}, lattice->atom_position(i)); } } +//This function calculates inelastic magnetic scattering (inelastic magnetic + elastic phonon scattering) +//That's why dirac_delta(f)*s0[i]*s0[j] is subtracted from cross_section(f, k). jams::MultiArray NeutronScatteringNoLatticeMonitor::calculate_unpolarized_cross_section(const jams::MultiArray &spectrum, const jams::MultiArray& elastic_spectrum) { @@ -195,7 +204,7 @@ void NeutronScatteringNoLatticeMonitor::shift_periodogram_overlap() { void NeutronScatteringNoLatticeMonitor::output_static_structure_factor() { const auto num_time_points = kspace_spins_timeseries_.size(0); - + //static_structure_factor=spin-spin correlation at the same time ofstream ofs(seedname + "_static_structure_factor_path_" + to_string(0) + ".tsv"); ofs << "index\t" << "qx\t" << "qy\t" << "qz\t" << "q_A-1\t"; @@ -206,9 +215,10 @@ void NeutronScatteringNoLatticeMonitor::output_static_structure_factor() { for (auto k = 0; k < kspace_path_.size(); ++k) { ofs << fmt::integer << k << "\t"; ofs << fmt::decimal << kspace_path_(k) << "\t"; - ofs << fmt::decimal << kTwoPi * norm(kspace_path_(k)) / (lattice->parameter() * 1e10) << "\t"; - for (auto i : {0,1,2}) { - for (auto j : {0,1,2}) { + //By norm(kspace_path_(k)) we can get the magnitude of each k + ofs << fmt::decimal << kTwoPi * norm(kspace_path_(k)) / (lattice->parameter() * 1e10) << "\t";//[angstrom^(-1)] + for (auto i : {0,1,2}) { //loop over i = 0, 1, 2 (meaning x,y,z components of spins) + for (auto j : {0,1,2}) {//loop over j = 0, 1, 2 (meaning x,y,z components of spins) auto s_a = static_structure_factor_(k)[i] / double(num_time_points); auto s_b = static_structure_factor_(k)[j] / double(num_time_points); auto s_ab = conj(s_a) * s_b; @@ -242,8 +252,8 @@ void NeutronScatteringNoLatticeMonitor::output_neutron_cross_section() { for (auto i = 0; i < (time_points / 2) + 1; ++i) { for (auto j = 0; j < kspace_path_.size(); ++j) { ofs << fmt::integer << j << "\t"; - ofs << fmt::decimal << kspace_path_(j) << "\t"; - ofs << fmt::decimal << kTwoPi * norm(kspace_path_(j)) / (lattice->parameter()*1e10) << "\t"; + ofs << fmt::decimal << kspace_path_(j) << "\t"; //k=k_j [AA^(-1)] + ofs << fmt::decimal << kTwoPi * norm(kspace_path_(j)) / (lattice->parameter()*1e10) << "\t"; //k=2*pi*k_j/a [m^(-1)] ofs << fmt::decimal << (i * freq_delta / 1e12) << "\t"; // THz ofs << fmt::decimal << (i * freq_delta / 1e12) * 4.135668 << "\t"; // meV // cross section output units are Barns Steradian^-1 Joules^-1 unitcell^-1 @@ -263,29 +273,32 @@ void NeutronScatteringNoLatticeMonitor::output_neutron_cross_section() { void NeutronScatteringNoLatticeMonitor::store_kspace_data_on_path() { auto i = periodogram_index_; - + //fill: insert the 3rd variable into 1st-2nd variables + //fill: "kspace_spins_timeseries_(i,0) - kspace_spins_timeseries_(i,0)+kspace_path_.size()" has Vec3cx{0,0} fill(&kspace_spins_timeseries_(i,0), &kspace_spins_timeseries_(i,0) + kspace_path_.size(), Vec3cx{0.0}); for (auto n = 0; n < globals::num_spins; ++n) { Vec3 spin = {globals::s(n,0), globals::s(n,1), globals::s(n,2)}; Vec3 r = rspace_displacement_(n); + //r should be normalized to get the length of position vector if (norm(r) >= 0.5) continue; auto delta_q = kspace_path_(1) - kspace_path_(0); auto window = 1.0; + //default setting of "bool do_rspace_windowing_" was changed to false. if (do_rspace_windowing_) { // blackmann 4 window const double a0 = 0.40217, a1 = 0.49704, a2 = 0.09392, a3 = 0.00183; const double x = (kTwoPi * norm(r)); window = a0 + a1 * cos(x) + a2 * cos(2 * x) + a3 * cos(3 * x); } - + //Fourier transform auto f0 = exp(-kImagTwoPi * dot(delta_q, r)); auto f = Complex{1.0, 0.0}; for (auto k = 0; k < kspace_path_.size(); ++k) { - kspace_spins_timeseries_(i, k) += f * spin * window; - f *= f0; + kspace_spins_timeseries_(i, k) += f * spin * window; + f *= f0; } } diff --git a/src/jams/monitors/neutron_scattering_no_lattice.h b/src/jams/monitors/neutron_scattering_no_lattice.h index 62e0b081..bb93eef0 100644 --- a/src/jams/monitors/neutron_scattering_no_lattice.h +++ b/src/jams/monitors/neutron_scattering_no_lattice.h @@ -34,7 +34,7 @@ class NeutronScatteringNoLatticeMonitor : public Monitor { jams::MultiArray calculate_unpolarized_cross_section(const jams::MultiArray& spectrum, const jams::MultiArray& elastic_spectrum); jams::MultiArray calculate_polarized_cross_sections(const jams::MultiArray& spectrum, const jams::MultiArray& elastic_spectrum, const std::vector& polarizations); - bool do_rspace_windowing_ = true; + bool do_rspace_windowing_ = false; jams::MultiArray rspace_displacement_; jams::MultiArray kspace_path_; jams::MultiArray kspace_spins_timeseries_; @@ -49,9 +49,9 @@ class NeutronScatteringNoLatticeMonitor : public Monitor { int periodogram_index_ = 0; int total_periods_ = 0; - double kmax_ = 50.0; - int num_k_ = 100; - Vec3 kvector_ = {0.0, 0.0, 1.0}; + double kmax_; + int num_k_; + Vec3 kvector_; }; From 1dfb8e9870b836fcdea3c1b2b14db2bb56903fab Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Thu, 6 Feb 2020 22:12:53 +0000 Subject: [PATCH 02/80] Added another comments. --- src/jams/monitors/neutron_scattering_no_lattice.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/jams/monitors/neutron_scattering_no_lattice.cc b/src/jams/monitors/neutron_scattering_no_lattice.cc index 5628b0ec..9a4739cd 100644 --- a/src/jams/monitors/neutron_scattering_no_lattice.cc +++ b/src/jams/monitors/neutron_scattering_no_lattice.cc @@ -79,7 +79,7 @@ void NeutronScatteringNoLatticeMonitor::configure_kspace_vectors(const libconfig num_k_ = jams::config_required(settings, "num_k"); kspace_path_.resize(num_k_); - for (auto i = 0; i < kspace_path_.size(); ++i) { + for (auto i = 0; i < kspace_path_.size(); ++i) {//this line decides the range of kvector (whether -kmax to +kmax or 0 to kmax) //kspace_path_ stores the scattering vectors kspace_path_(i) = kvector_ * i * (kmax_ / (num_k_-1));//kmax shoul be in the unit of [1/a] cout << "i = " << i << ", kspace_path_(i) = " << kspace_path_(i) << endl; From 5514e947244f855f1b77fe7d2f693ca80610655f Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Fri, 7 Feb 2020 15:08:05 +0000 Subject: [PATCH 03/80] Fixed Segmentation fault & now we can specity whether we use the form factors. --- src/jams/monitors/neutron_scattering_no_lattice.cc | 10 +++++++++- src/jams/monitors/neutron_scattering_no_lattice.h | 2 +- 2 files changed, 10 insertions(+), 2 deletions(-) diff --git a/src/jams/monitors/neutron_scattering_no_lattice.cc b/src/jams/monitors/neutron_scattering_no_lattice.cc index 9a4739cd..5e3044f3 100644 --- a/src/jams/monitors/neutron_scattering_no_lattice.cc +++ b/src/jams/monitors/neutron_scattering_no_lattice.cc @@ -26,6 +26,7 @@ NeutronScatteringNoLatticeMonitor::NeutronScatteringNoLatticeMonitor(const libco if (settings.exists("form_factor")) { configure_form_factors(settings["form_factor"]); + evaluate_form_factor_ = true; } if (settings.exists("polarizations")) { @@ -60,6 +61,7 @@ void NeutronScatteringNoLatticeMonitor::update(Solver *solver) { shift_periodogram_overlap(); total_periods_++; + //element_sum: sum up each matrix element element_sum(total_unpolarized_neutron_cross_section_, calculate_unpolarized_cross_section(spectrum, elastic_spectrum)); if (!neutron_polarizations_.empty()) { @@ -109,8 +111,14 @@ NeutronScatteringNoLatticeMonitor::calculate_unpolarized_cross_section(const jam auto Q = unit_vector(kspace_path_(k)); auto s_a = conj(spectrum(f, k)); auto s_b = spectrum(f, k); + double ff = 0.0; - auto ff = pow2(neutron_form_factors_(0, k)); // NOTE: currently only supports one material + if (evaluate_form_factor_) { + ff = pow2(neutron_form_factors_(0, k)); // NOTE: currently only supports one material + } + else{ + ff = 1.0; + } for (auto i : {0, 1, 2}) { for (auto j : {0, 1, 2}) { cross_section(f, k) += ff * (kronecker_delta(i, j) - Q[i] * Q[j]) * ((s_a[i] * s_b[j]) - dirac_delta(f) * s0[i] * s0[j]); diff --git a/src/jams/monitors/neutron_scattering_no_lattice.h b/src/jams/monitors/neutron_scattering_no_lattice.h index bb93eef0..e423da0d 100644 --- a/src/jams/monitors/neutron_scattering_no_lattice.h +++ b/src/jams/monitors/neutron_scattering_no_lattice.h @@ -52,7 +52,7 @@ class NeutronScatteringNoLatticeMonitor : public Monitor { double kmax_; int num_k_; Vec3 kvector_; - + bool evaluate_form_factor_ = false; }; From 36f7a2114071028b9600d4b778e97e6cc59c6d60 Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Wed, 12 Feb 2020 19:41:53 +0000 Subject: [PATCH 04/80] Some more comments are added. --- src/jams/interface/fft.h | 1 + src/jams/monitors/neutron_scattering_no_lattice.cc | 8 ++++++++ 2 files changed, 9 insertions(+) diff --git a/src/jams/interface/fft.h b/src/jams/interface/fft.h index 04985580..e10e9c7a 100644 --- a/src/jams/interface/fft.h +++ b/src/jams/interface/fft.h @@ -23,6 +23,7 @@ using Complex = std::complex; +//(*) where are these values used? namespace jams { struct PeriodogramProps { int length = 1000; diff --git a/src/jams/monitors/neutron_scattering_no_lattice.cc b/src/jams/monitors/neutron_scattering_no_lattice.cc index 5e3044f3..f7762c7c 100644 --- a/src/jams/monitors/neutron_scattering_no_lattice.cc +++ b/src/jams/monitors/neutron_scattering_no_lattice.cc @@ -121,6 +121,8 @@ NeutronScatteringNoLatticeMonitor::calculate_unpolarized_cross_section(const jam } for (auto i : {0, 1, 2}) { for (auto j : {0, 1, 2}) { + //(*) We may have to remove dirac_delta(f) since this term should be subtracted all the time? + //(*) even when f is not equals to 0 (I think elastic terms are always taken into consideration) cross_section(f, k) += ff * (kronecker_delta(i, j) - Q[i] * Q[j]) * ((s_a[i] * s_b[j]) - dirac_delta(f) * s0[i] * s0[j]); } } @@ -212,6 +214,8 @@ void NeutronScatteringNoLatticeMonitor::shift_periodogram_overlap() { void NeutronScatteringNoLatticeMonitor::output_static_structure_factor() { const auto num_time_points = kspace_spins_timeseries_.size(0); + //(*) I think we should take the thermal average using python after running the program. + //(*) In that case, don't we want to save files with different name for each periodgram (to_string(0-total_periods))? //static_structure_factor=spin-spin correlation at the same time ofstream ofs(seedname + "_static_structure_factor_path_" + to_string(0) + ".tsv"); @@ -251,6 +255,10 @@ void NeutronScatteringNoLatticeMonitor::output_neutron_cross_section() { // sample time is here because the fourier transform in time is not an integral // but a discrete sum + //(*) Don't we have to remove (periodogram_props_.sample_time/kTwoPi) cause there shuld be + // another "1/periodogram_props_.sample_time" which comes from Dirac's delta function (Fermi's golden rule)? + //(*) Also we have to multiply kGyromagneticRatio to the prefactor? + //(*) Also we have to remove (0.5 * kNeutronGFactor) since it is already considered in the form factor? auto prefactor = (periodogram_props_.sample_time / double(total_periods_)) * (1.0 / (kTwoPi * kHBar)) * pow2((0.5 * kNeutronGFactor * pow2(kElementaryCharge)) / (kElectronMass * pow2(kSpeedOfLight))); auto barns_unitcell = prefactor / (1e-28); From ba156e24357bd23e28e01489f12afafa7c453f4d Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Fri, 28 Feb 2020 17:38:51 +0000 Subject: [PATCH 05/80] Fixed exchange interaction matrix to be diagonal. --- src/jams/hamiltonian/sparse_interaction.cc | 4 +--- src/jams/monitors/neutron_scattering_no_lattice.cc | 6 ------ 2 files changed, 1 insertion(+), 9 deletions(-) diff --git a/src/jams/hamiltonian/sparse_interaction.cc b/src/jams/hamiltonian/sparse_interaction.cc index ef1bcfd7..82b37c20 100644 --- a/src/jams/hamiltonian/sparse_interaction.cc +++ b/src/jams/hamiltonian/sparse_interaction.cc @@ -19,9 +19,7 @@ void SparseInteractionHamiltonian::insert_interaction_scalar(const int i, const return; } for (auto m = 0; m < 3; ++m) { - for (auto n = 0; n < 3; ++n) { - sparse_matrix_builder_.insert(3 * i + m, 3 * j + n, value); - } + sparse_matrix_builder_.insert(3 * i + m, 3 * j + m, value); } } diff --git a/src/jams/monitors/neutron_scattering_no_lattice.cc b/src/jams/monitors/neutron_scattering_no_lattice.cc index f7762c7c..5ae9ff2f 100644 --- a/src/jams/monitors/neutron_scattering_no_lattice.cc +++ b/src/jams/monitors/neutron_scattering_no_lattice.cc @@ -214,8 +214,6 @@ void NeutronScatteringNoLatticeMonitor::shift_periodogram_overlap() { void NeutronScatteringNoLatticeMonitor::output_static_structure_factor() { const auto num_time_points = kspace_spins_timeseries_.size(0); - //(*) I think we should take the thermal average using python after running the program. - //(*) In that case, don't we want to save files with different name for each periodgram (to_string(0-total_periods))? //static_structure_factor=spin-spin correlation at the same time ofstream ofs(seedname + "_static_structure_factor_path_" + to_string(0) + ".tsv"); @@ -255,10 +253,6 @@ void NeutronScatteringNoLatticeMonitor::output_neutron_cross_section() { // sample time is here because the fourier transform in time is not an integral // but a discrete sum - //(*) Don't we have to remove (periodogram_props_.sample_time/kTwoPi) cause there shuld be - // another "1/periodogram_props_.sample_time" which comes from Dirac's delta function (Fermi's golden rule)? - //(*) Also we have to multiply kGyromagneticRatio to the prefactor? - //(*) Also we have to remove (0.5 * kNeutronGFactor) since it is already considered in the form factor? auto prefactor = (periodogram_props_.sample_time / double(total_periods_)) * (1.0 / (kTwoPi * kHBar)) * pow2((0.5 * kNeutronGFactor * pow2(kElementaryCharge)) / (kElectronMass * pow2(kSpeedOfLight))); auto barns_unitcell = prefactor / (1e-28); From 2852051c3194f11fb958d1216d0017186838f902 Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Wed, 4 Mar 2020 20:59:24 +0000 Subject: [PATCH 06/80] Added step function into exchange_functional for test run. --- src/jams/hamiltonian/exchange_functional.cc | 12 +++++++++++- src/jams/hamiltonian/exchange_functional.h | 1 + src/jams/hamiltonian/sparse_interaction.cc | 1 + 3 files changed, 13 insertions(+), 1 deletion(-) diff --git a/src/jams/hamiltonian/exchange_functional.cc b/src/jams/hamiltonian/exchange_functional.cc index e767affc..b7929139 100644 --- a/src/jams/hamiltonian/exchange_functional.cc +++ b/src/jams/hamiltonian/exchange_functional.cc @@ -61,6 +61,14 @@ double ExchangeFunctionalHamiltonian::functional_kaneyoshi(const double rij, con return J0 * pow2(rij - r0) * exp(-pow2(rij - r0) / (2 * pow2(sigma))); } +double ExchangeFunctionalHamiltonian::functional_step(const double rij, const double J0, const double r_out){ + if (rij < r_out){ + return J0; + } + else + return 0.0; +} + ExchangeFunctionalHamiltonian::ExchangeFunctional ExchangeFunctionalHamiltonian::functional_from_settings(const libconfig::Setting &settings) { using namespace std::placeholders; @@ -75,7 +83,9 @@ ExchangeFunctionalHamiltonian::functional_from_settings(const libconfig::Setting } else if (functional_name == "gaussian") { return bind(functional_gaussian, _1, double(settings["J0"]), double(settings["r0"]), double(settings["sigma"])); } else if (functional_name == "kaneyoshi") { - return bind(functional_kaneyoshi, _1, double(settings["J0"]), double(settings["r0"]), double(settings["sigma"])); + return bind(functional_kaneyoshi, _1, double(settings["J0"]), double(settings["r0"]), double(settings["sigma"])); + } else if (functional_name == "step") { + return bind(functional_step, _1, double(settings["J0"]), double(settings["r_out"])); } else { throw runtime_error("unknown exchange functional: " + functional_name); } diff --git a/src/jams/hamiltonian/exchange_functional.h b/src/jams/hamiltonian/exchange_functional.h index 7c5d9c5e..24ff31fe 100644 --- a/src/jams/hamiltonian/exchange_functional.h +++ b/src/jams/hamiltonian/exchange_functional.h @@ -23,6 +23,7 @@ class ExchangeFunctionalHamiltonian : public SparseInteractionHamiltonian { static double functional_exp(double rij, double J0, double r0, double lengthscale); static double functional_gaussian(double rij, double J0, double r0, double lengthscale); static double functional_kaneyoshi(double rij, double J0, double r0, double lengthscale); + static double functional_step(double rij, double J0, double r_out); double radius_cutoff_; }; diff --git a/src/jams/hamiltonian/sparse_interaction.cc b/src/jams/hamiltonian/sparse_interaction.cc index 82b37c20..38616ce9 100644 --- a/src/jams/hamiltonian/sparse_interaction.cc +++ b/src/jams/hamiltonian/sparse_interaction.cc @@ -20,6 +20,7 @@ void SparseInteractionHamiltonian::insert_interaction_scalar(const int i, const } for (auto m = 0; m < 3; ++m) { sparse_matrix_builder_.insert(3 * i + m, 3 * j + m, value); +// std::cout << "i= " << i << " j= " << j << " : row = " << 3 * i + m << ", colum = " << 3 * j + m << ", value [J] = " << value/input_unit_conversion_ << std::endl; } } From dfec800721ed9de9e2671f4a97d255d1ac340265 Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Fri, 6 Mar 2020 17:55:27 +0000 Subject: [PATCH 07/80] Fixed the size of kspace_path_ array. --- src/jams/monitors/neutron_scattering_no_lattice.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/jams/monitors/neutron_scattering_no_lattice.cc b/src/jams/monitors/neutron_scattering_no_lattice.cc index 5ae9ff2f..4d60d51b 100644 --- a/src/jams/monitors/neutron_scattering_no_lattice.cc +++ b/src/jams/monitors/neutron_scattering_no_lattice.cc @@ -80,10 +80,10 @@ void NeutronScatteringNoLatticeMonitor::configure_kspace_vectors(const libconfig kvector_ = jams::config_required(settings, "kvector"); num_k_ = jams::config_required(settings, "num_k"); - kspace_path_.resize(num_k_); + kspace_path_.resize(num_k_+1); for (auto i = 0; i < kspace_path_.size(); ++i) {//this line decides the range of kvector (whether -kmax to +kmax or 0 to kmax) //kspace_path_ stores the scattering vectors - kspace_path_(i) = kvector_ * i * (kmax_ / (num_k_-1));//kmax shoul be in the unit of [1/a] + kspace_path_(i) = kvector_ * i * (kmax_ / num_k_);//kmax shoul be in the unit of [2*pi/a(lattice constant)] cout << "i = " << i << ", kspace_path_(i) = " << kspace_path_(i) << endl; } From eebde26783be8ed58d95c1cb461864ae32aeb835 Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Sat, 7 Mar 2020 00:15:24 +0000 Subject: [PATCH 08/80] Fixed possible maximum value of distance in r_space and added new function Lattice::get_lattice_dimensions(). --- src/jams/core/lattice.cc | 3 +++ src/jams/core/lattice.h | 2 ++ src/jams/monitors/neutron_scattering_no_lattice.cc | 4 +++- 3 files changed, 8 insertions(+), 1 deletion(-) diff --git a/src/jams/core/lattice.cc b/src/jams/core/lattice.cc index 0b09bb8e..bd965247 100755 --- a/src/jams/core/lattice.cc +++ b/src/jams/core/lattice.cc @@ -1120,3 +1120,6 @@ bool Lattice::has_impurities() const { return !impurity_map_.empty(); } +Vec3i Lattice::get_lattice_dimensions() { + return lattice_dimensions; +} \ No newline at end of file diff --git a/src/jams/core/lattice.h b/src/jams/core/lattice.h index 31982643..f6df0f41 100755 --- a/src/jams/core/lattice.h +++ b/src/jams/core/lattice.h @@ -128,6 +128,8 @@ class Lattice : public Base { const Vec3i &kspace_size() const; + Vec3i get_lattice_dimensions(); + private: void read_materials_from_config(const libconfig::Setting &settings); ImpurityMap read_impurities_from_config(const libconfig::Setting &settings); diff --git a/src/jams/monitors/neutron_scattering_no_lattice.cc b/src/jams/monitors/neutron_scattering_no_lattice.cc index 4d60d51b..953260d0 100644 --- a/src/jams/monitors/neutron_scattering_no_lattice.cc +++ b/src/jams/monitors/neutron_scattering_no_lattice.cc @@ -292,7 +292,9 @@ void NeutronScatteringNoLatticeMonitor::store_kspace_data_on_path() { Vec3 r = rspace_displacement_(n); //r should be normalized to get the length of position vector - if (norm(r) >= 0.5) continue; + Vec3i lattice_dimensions = ::lattice->get_lattice_dimensions(); + //Here we assume lattice_dimensions are all the same (x,y,z) + if (norm(r) >= lattice_dimensions[0]*0.5) continue; auto delta_q = kspace_path_(1) - kspace_path_(0); auto window = 1.0; From 51c49116fb2ba43843ae6f7130d6c60b8d35ff04 Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Wed, 11 Mar 2020 20:00:04 +0000 Subject: [PATCH 09/80] Added multiple (up to three) Gaussian exchange functional. --- src/jams/hamiltonian/exchange_functional.cc | 6 ++++++ src/jams/hamiltonian/exchange_functional.h | 1 + 2 files changed, 7 insertions(+) diff --git a/src/jams/hamiltonian/exchange_functional.cc b/src/jams/hamiltonian/exchange_functional.cc index b7929139..74b51f65 100644 --- a/src/jams/hamiltonian/exchange_functional.cc +++ b/src/jams/hamiltonian/exchange_functional.cc @@ -69,6 +69,10 @@ double ExchangeFunctionalHamiltonian::functional_step(const double rij, const do return 0.0; } +double ExchangeFunctionalHamiltonian::functional_gaussian_multi(const double rij, const double J0, const double r0, const double sigma, const double J0_2, const double r0_2, const double sigma_2, const double J0_3, const double r0_3, const double sigma_3){ + return J0 * exp(-pow2(rij - r0)/(2 * pow2(sigma))) + J0_2 * exp(-pow2(rij - r0_2)/(2 * pow2(sigma_2))) + J0_3 * exp(-pow2(rij - r0_3)/(2 * pow2(sigma_3))); +} + ExchangeFunctionalHamiltonian::ExchangeFunctional ExchangeFunctionalHamiltonian::functional_from_settings(const libconfig::Setting &settings) { using namespace std::placeholders; @@ -82,6 +86,8 @@ ExchangeFunctionalHamiltonian::functional_from_settings(const libconfig::Setting return bind(functional_exp, _1, double(settings["J0"]), double(settings["r0"]), double(settings["sigma"])); } else if (functional_name == "gaussian") { return bind(functional_gaussian, _1, double(settings["J0"]), double(settings["r0"]), double(settings["sigma"])); + } else if (functional_name == "gaussian_multi") { + return bind(functional_gaussian_multi, _1, double(settings["J0"]), double(settings["r0"]), double(settings["sigma"]), double(settings["J0_2"]), double(settings["r0_2"]), double(settings["sigma_2"]), double(settings["J0_3"]), double(settings["r0_3"]), double(settings["sigma_3"])); } else if (functional_name == "kaneyoshi") { return bind(functional_kaneyoshi, _1, double(settings["J0"]), double(settings["r0"]), double(settings["sigma"])); } else if (functional_name == "step") { diff --git a/src/jams/hamiltonian/exchange_functional.h b/src/jams/hamiltonian/exchange_functional.h index 24ff31fe..decbc3e5 100644 --- a/src/jams/hamiltonian/exchange_functional.h +++ b/src/jams/hamiltonian/exchange_functional.h @@ -22,6 +22,7 @@ class ExchangeFunctionalHamiltonian : public SparseInteractionHamiltonian { static double functional_rkky(double rij, double J0, double r0, double k_F); static double functional_exp(double rij, double J0, double r0, double lengthscale); static double functional_gaussian(double rij, double J0, double r0, double lengthscale); + static double functional_gaussian_multi(double rij, double J0, double r0, double lengthscale, double J0_2, double r0_2, double lengthscale_2, double J0_3, double r0_3, double lengthscale_3); static double functional_kaneyoshi(double rij, double J0, double r0, double lengthscale); static double functional_step(double rij, double J0, double r_out); From 201e9127ed15a49ee72ee3acbcf5c2aa4d716192 Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Thu, 12 Mar 2020 22:58:45 +0000 Subject: [PATCH 10/80] Added real-space spin-spin correlation function into the NoLattice Monitor. --- .../monitors/neutron_scattering_no_lattice.cc | 185 +++++++++++++++++- .../monitors/neutron_scattering_no_lattice.h | 12 +- 2 files changed, 192 insertions(+), 5 deletions(-) diff --git a/src/jams/monitors/neutron_scattering_no_lattice.cc b/src/jams/monitors/neutron_scattering_no_lattice.cc index 953260d0..2427d153 100644 --- a/src/jams/monitors/neutron_scattering_no_lattice.cc +++ b/src/jams/monitors/neutron_scattering_no_lattice.cc @@ -38,6 +38,8 @@ NeutronScatteringNoLatticeMonitor::NeutronScatteringNoLatticeMonitor(const libco configure_periodogram(settings["periodogram"]); } + t_diff_r_dep_ = config_required(settings, "t_diff_r_dep_"); + periodogram_props_.sample_time = output_step_freq_ * solver->time_step(); //zero: set all the elements as 0 @@ -53,6 +55,22 @@ void NeutronScatteringNoLatticeMonitor::update(Solver *solver) { store_kspace_data_on_path(); periodogram_index_++; + using namespace globals; + + jams::MultiArray spin_x_time(num_spins); + jams::MultiArray spin_y_time(num_spins); + jams::MultiArray spin_z_time(num_spins); + + for(auto i = 0; i < num_spins; ++i){ + spin_x_time(i) = s(i, 0); + spin_y_time(i) = s(i, 1); + spin_z_time(i) = s(i, 2); + } + + spin_x.push_back(spin_x_time); + spin_y.push_back(spin_y_time); + spin_z.push_back(spin_z_time); + if (is_multiple_of(periodogram_index_, periodogram_props_.length)) { auto elastic_spectrum = average_kspace_timeseries(); @@ -96,13 +114,13 @@ void NeutronScatteringNoLatticeMonitor::configure_kspace_vectors(const libconfig //This function calculates inelastic magnetic scattering (inelastic magnetic + elastic phonon scattering) //That's why dirac_delta(f)*s0[i]*s0[j] is subtracted from cross_section(f, k). -jams::MultiArray +jams::MultiArray, 2> NeutronScatteringNoLatticeMonitor::calculate_unpolarized_cross_section(const jams::MultiArray &spectrum, const jams::MultiArray& elastic_spectrum) { const auto num_freqencies = spectrum.size(0); const auto num_reciprocal_points = kspace_path_.size(); - jams::MultiArray cross_section(num_freqencies, num_reciprocal_points); + jams::MultiArray, 2> cross_section(num_freqencies, num_reciprocal_points); cross_section.zero(); for (auto f = 0; f < num_freqencies; ++f) { @@ -131,14 +149,14 @@ NeutronScatteringNoLatticeMonitor::calculate_unpolarized_cross_section(const jam return cross_section; } -jams::MultiArray +jams::MultiArray, 3> NeutronScatteringNoLatticeMonitor::calculate_polarized_cross_sections(const MultiArray &spectrum, const jams::MultiArray& elastic_spectrum, const vector &polarizations) { const auto num_freqencies = spectrum.size(0); const auto num_reciprocal_points = kspace_path_.size(); - MultiArray convolved(polarizations.size(), num_freqencies, num_reciprocal_points); + MultiArray, 3> convolved(polarizations.size(), num_freqencies, num_reciprocal_points); convolved.zero(); for (auto f = 0; f < num_freqencies; ++f) { @@ -376,3 +394,162 @@ jams::MultiArray NeutronScatteringNoLatticeMonitor::average_kspace_ti return average; } + +NeutronScatteringNoLatticeMonitor::~NeutronScatteringNoLatticeMonitor() { + using namespace std; + using namespace globals; + + const double delta_t = periodogram_props_.sample_time; + + string name3 = seedname + "_S+S-_r_dependence.tsv"; + ofstream outfile3(name3.c_str()); + outfile3.setf(std::ios::right); + + string name4 = seedname + "_radial_dist.tsv"; + ofstream outfile4(name4.c_str()); + outfile4.setf(std::ios::right); + + string name5 = seedname + "_SzSz_r_dependence.tsv"; + ofstream outfile5(name5.c_str()); + outfile5.setf(std::ios::right); + + // header for the r-dependence of the S+S- file + outfile3 << std::setw(20) << "time_diff" << "\t";//(\t=tab) + outfile3 << std::setw(20) << "r_ij" << "\t";//(\t=tab) + outfile3 << std::setw(20) << "Re:S+S-/dist" << "\t"; + outfile3 << std::setw(20) << "Im:S+S-/dist" << "\n"; + + // header for the radial dist file + outfile4 << std::setw(20) << "bin" << "\t"; + outfile4 << std::setw(20) << "dist" << "\n"; + + // header for the r-dependence of the SzSz file + outfile5 << std::setw(20) << "time_diff" << "\t";//(\t=tab) + outfile5 << std::setw(20) << "r_ij" << "\t";//(\t=tab) + outfile5 << std::setw(20) << "Re:SzSz/dist" << "\t"; + outfile5 << std::setw(20) << "Im:SzSz/dist" << "\n"; + + Vec3i lattice_dimensions = ::lattice->get_lattice_dimensions(); + int num_distance = std::ceil(0.5 * lattice_dimensions[0] / delta_r_); + //size of the array: (L/2)/delta_r_ + + vector> r_dependence_SxSy(num_distance, 0.0); + vector> r_dependence_SzSz(num_distance, 0.0); + vector count_r_dep(r_dependence_SxSy.size(),0); + + auto radial_dist = radial_distribution_function(); + + for (unsigned j = 0; j < num_spins; ++j) { + const auto r_j = lattice->atom_position(j); + for (unsigned i = 0; i < num_spins; ++i) { + if(i == j){ + continue; + } + const auto r_i = lattice->atom_position(i); + const auto r_ij = lattice->displacement(r_i, r_j); // using j,i in this order gives r_j - r_i + + unsigned m = std::ceil(this->distance(r_ij)/delta_r_); + if(m < r_dependence_SxSy.size()) { + auto correlation = time_correlation(i, j); + + r_dependence_SxSy[m] += correlation[t_diff_r_dep_]; + auto num_time_samples_z = (unsigned int) (periodogram_props_.length); + for (auto t_ref = 0; t_ref < num_time_samples_z; t_ref++) { + r_dependence_SzSz[m] += spin_z[t_ref](i) * spin_z[t_ref + t_diff_r_dep_](j); + } + count_r_dep[m]++; + } + }//i + }//j + + for (auto m = 0; m < radial_dist.size(); m++) { + outfile4 << std::setw(20) << m << "\t"; + outfile4 << std::setw(20) << radial_dist[m] << "\n"; + } + + for (auto m = 1; m < r_dependence_SxSy.size(); m++) { + outfile3 << std::setw(20) << t_diff_r_dep_ * delta_t << "\t"; + outfile3 << std::setw(20) << m << "\t"; + if (count_r_dep[m] != 0) { + outfile3 << std::setw(20) << r_dependence_SxSy[m].real() / count_r_dep[m] << "\t"; + outfile3 << std::setw(20) << r_dependence_SxSy[m].imag() / count_r_dep[m] << "\n"; + } + else{ + outfile3 << std::setw(20) << 0.0 << "\t"; + outfile3 << std::setw(20) << 0.0 << "\n"; + } + } + + for (auto m = 1; m < r_dependence_SzSz.size(); m++) { + outfile5 << std::setw(20) << t_diff_r_dep_ * delta_t << "\t"; + outfile5 << std::setw(20) << m << "\t"; + if (count_r_dep[m] != 0) { + outfile5 << std::setw(20) << r_dependence_SzSz[m].real() / count_r_dep[m] << "\t"; + outfile5 << std::setw(20) << r_dependence_SzSz[m].imag() / count_r_dep[m] << "\n"; + } + else{ + outfile5 << std::setw(20) << 0.0 << "\t"; + outfile5 << std::setw(20) << 0.0 << "\n"; + } + } + + outfile3.close(); + outfile4.close(); + outfile5.close(); + +} + +double NeutronScatteringNoLatticeMonitor::distance(const Vec3 &r_ij) { + return std::sqrt(r_ij[0] * r_ij[0] + r_ij[1] * r_ij[1] + r_ij[2] * r_ij[2]); +} + +std::vector> NeutronScatteringNoLatticeMonitor::time_correlation(unsigned i, unsigned j) { + const auto num_time_samples = (unsigned int)(periodogram_props_.length); + std::vector> correlation(num_time_samples, 0.0); + + for (auto t_diff = 0; t_diff < correlation.size(); t_diff++) { + for (auto t_ref = 0; t_ref < num_time_samples; t_ref++) { + if ((t_ref + t_diff) < num_time_samples && t_ref < num_time_samples) { + correlation[t_diff] += + spin_x[t_ref](i) * spin_x[t_ref + t_diff](j) + spin_y[t_ref](i) * spin_y[t_ref + t_diff](j) + + kImagOne * (-spin_x[t_ref](i) * spin_y[t_ref + t_diff](j) + spin_y[t_ref](i) * spin_x[t_ref + t_diff](j)); + +// if(t_diff == 0){ +// std::cout << "num_t_samples = " << num_time_samples << std::endl; +// std::cout << "t_diff = " << t_diff << ", t_ref = " << t_ref << std::endl; +// std::cout << "corr (Re) = " << correlation[t_diff].real() << ", corr (Im) = " << correlation[t_diff].imag() << std::endl; +// } + + } + } + } + return correlation; +} + +std::vector NeutronScatteringNoLatticeMonitor::radial_distribution_function() { + using namespace globals; + using namespace std; + Vec3i lattice_dimensions = ::lattice->get_lattice_dimensions(); + //Here we assume lattice_dimensions are all the same (x,y,z) + std::vector dist(std::ceil(0.5 * lattice_dimensions[0] / delta_r_),0.0); + for(unsigned j = 0; j < num_spins; ++j){ + const auto r_j = lattice->atom_position(j); + for(unsigned i = 0; i < num_spins; ++i){ + if(i == j){ + dist[0] = 0; + continue; + } + const auto r_i = lattice->atom_position(i); + const auto r_ij = (lattice->displacement(r_i, r_j)); // using j,i in this order gives r_j - r_i + int m = std::ceil(std::sqrt(r_ij[0] * r_ij[0] + r_ij[1] * r_ij[1] + r_ij[2] * r_ij[2])/delta_r_); + if(m > dist.size()){ + continue; + } + dist[m]++; + } + } + for (auto kk = 0; kk < dist.size(); kk++){ + dist[kk] /= num_spins; + } + return dist; +} \ No newline at end of file diff --git a/src/jams/monitors/neutron_scattering_no_lattice.h b/src/jams/monitors/neutron_scattering_no_lattice.h index e423da0d..67f85f4d 100644 --- a/src/jams/monitors/neutron_scattering_no_lattice.h +++ b/src/jams/monitors/neutron_scattering_no_lattice.h @@ -11,7 +11,8 @@ class NeutronScatteringNoLatticeMonitor : public Monitor { public: explicit NeutronScatteringNoLatticeMonitor(const libconfig::Setting &settings); - ~NeutronScatteringNoLatticeMonitor() override = default; + ~NeutronScatteringNoLatticeMonitor() override; +// ~NeutronScatteringNoLatticeMonitor() override = default; void post_process() override {}; void update(Solver *solver) override; @@ -53,7 +54,16 @@ class NeutronScatteringNoLatticeMonitor : public Monitor { int num_k_; Vec3 kvector_; bool evaluate_form_factor_ = false; + double delta_r_ = 0.001; + int t_diff_r_dep_; + std::vector> spin_x; + std::vector> spin_y; + std::vector> spin_z; + + std::vector radial_distribution_function(); + double distance(const Vec3 &r_ij); + std::vector> time_correlation(unsigned i, unsigned j); }; #endif //JAMS_NEUTRON_SCATTERING_NO_LATTICE_H From 9fec5cbf59b0c4d2e1c16cc0a1ad38fdf2732f65 Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Thu, 23 Apr 2020 16:21:16 +0900 Subject: [PATCH 11/80] Generalized calculating rspace_displacement vectors so that we can impose open boundary --- src/jams/monitors/neutron_scattering_no_lattice.cc | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/jams/monitors/neutron_scattering_no_lattice.cc b/src/jams/monitors/neutron_scattering_no_lattice.cc index 2427d153..accb719a 100644 --- a/src/jams/monitors/neutron_scattering_no_lattice.cc +++ b/src/jams/monitors/neutron_scattering_no_lattice.cc @@ -106,9 +106,12 @@ void NeutronScatteringNoLatticeMonitor::configure_kspace_vectors(const libconfig } rspace_displacement_.resize(globals::s.size(0)); + Vec3i lattice_dimensions = ::lattice->get_lattice_dimensions(); for (auto i = 0; i < globals::s.size(0); ++i) { //rspace_displacement_ is the position vector of i-th spin - rspace_displacement_(i) = lattice->displacement({0,0,0}, lattice->atom_position(i)); + //rspace_displacement_(i) = lattice->displacement({0,0,0}, lattice->atom_position(i)); + //Generalize so that we can impose open boundary + rspace_displacement_(i) = lattice->displacement({lattice_dimensions[0]*0.5,lattice_dimensions[1]*0.5,lattice_dimensions[2]*0.5}, lattice->atom_position(i)); } } From bcec3c12a620776ff8f9bc39b15f1868464644cc Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Thu, 23 Apr 2020 16:43:05 +0900 Subject: [PATCH 12/80] Explicitly output if we do the rspace_windowing and calculating the real-space spin-spin correlations were changed to be optional. --- .../monitors/neutron_scattering_no_lattice.cc | 222 +++++++++--------- .../monitors/neutron_scattering_no_lattice.h | 1 + 2 files changed, 117 insertions(+), 106 deletions(-) diff --git a/src/jams/monitors/neutron_scattering_no_lattice.cc b/src/jams/monitors/neutron_scattering_no_lattice.cc index accb719a..812ac654 100644 --- a/src/jams/monitors/neutron_scattering_no_lattice.cc +++ b/src/jams/monitors/neutron_scattering_no_lattice.cc @@ -23,6 +23,10 @@ NeutronScatteringNoLatticeMonitor::NeutronScatteringNoLatticeMonitor(const libco configure_kspace_vectors(settings); config_optional(settings, "rspace_windowing", do_rspace_windowing_); + cout << "do_rspace_windowing_ = " << do_rspace_windowing_ << endl; + + config_optional(settings, "output_real_space_ss_corr_", output_real_space_ss_corr_); + cout << "output_real_space_ss_corr_ = " << output_real_space_ss_corr_ << endl; if (settings.exists("form_factor")) { configure_form_factors(settings["form_factor"]); @@ -38,7 +42,9 @@ NeutronScatteringNoLatticeMonitor::NeutronScatteringNoLatticeMonitor(const libco configure_periodogram(settings["periodogram"]); } - t_diff_r_dep_ = config_required(settings, "t_diff_r_dep_"); + if (output_real_space_ss_corr_) { + t_diff_r_dep_ = config_required(settings, "t_diff_r_dep_"); + } periodogram_props_.sample_time = output_step_freq_ * solver->time_step(); @@ -55,21 +61,23 @@ void NeutronScatteringNoLatticeMonitor::update(Solver *solver) { store_kspace_data_on_path(); periodogram_index_++; - using namespace globals; + if (output_real_space_ss_corr_) { + using namespace globals; - jams::MultiArray spin_x_time(num_spins); - jams::MultiArray spin_y_time(num_spins); - jams::MultiArray spin_z_time(num_spins); + jams::MultiArray spin_x_time(num_spins); + jams::MultiArray spin_y_time(num_spins); + jams::MultiArray spin_z_time(num_spins); - for(auto i = 0; i < num_spins; ++i){ - spin_x_time(i) = s(i, 0); - spin_y_time(i) = s(i, 1); - spin_z_time(i) = s(i, 2); - } + for(auto i = 0; i < num_spins; ++i){ + spin_x_time(i) = s(i, 0); + spin_y_time(i) = s(i, 1); + spin_z_time(i) = s(i, 2); + } - spin_x.push_back(spin_x_time); - spin_y.push_back(spin_y_time); - spin_z.push_back(spin_z_time); + spin_x.push_back(spin_x_time); + spin_y.push_back(spin_y_time); + spin_z.push_back(spin_z_time); + } if (is_multiple_of(periodogram_index_, periodogram_props_.length)) { @@ -325,6 +333,7 @@ void NeutronScatteringNoLatticeMonitor::store_kspace_data_on_path() { const double a0 = 0.40217, a1 = 0.49704, a2 = 0.09392, a3 = 0.00183; const double x = (kTwoPi * norm(r)); window = a0 + a1 * cos(x) + a2 * cos(2 * x) + a3 * cos(3 * x); + cout << "do_rspace_windowing_ = true, value = " << do_rspace_windowing_ << endl; } //Fourier transform auto f0 = exp(-kImagTwoPi * dot(delta_q, r)); @@ -399,107 +408,108 @@ jams::MultiArray NeutronScatteringNoLatticeMonitor::average_kspace_ti } NeutronScatteringNoLatticeMonitor::~NeutronScatteringNoLatticeMonitor() { - using namespace std; - using namespace globals; - - const double delta_t = periodogram_props_.sample_time; - - string name3 = seedname + "_S+S-_r_dependence.tsv"; - ofstream outfile3(name3.c_str()); - outfile3.setf(std::ios::right); - - string name4 = seedname + "_radial_dist.tsv"; - ofstream outfile4(name4.c_str()); - outfile4.setf(std::ios::right); - - string name5 = seedname + "_SzSz_r_dependence.tsv"; - ofstream outfile5(name5.c_str()); - outfile5.setf(std::ios::right); - - // header for the r-dependence of the S+S- file - outfile3 << std::setw(20) << "time_diff" << "\t";//(\t=tab) - outfile3 << std::setw(20) << "r_ij" << "\t";//(\t=tab) - outfile3 << std::setw(20) << "Re:S+S-/dist" << "\t"; - outfile3 << std::setw(20) << "Im:S+S-/dist" << "\n"; - - // header for the radial dist file - outfile4 << std::setw(20) << "bin" << "\t"; - outfile4 << std::setw(20) << "dist" << "\n"; - - // header for the r-dependence of the SzSz file - outfile5 << std::setw(20) << "time_diff" << "\t";//(\t=tab) - outfile5 << std::setw(20) << "r_ij" << "\t";//(\t=tab) - outfile5 << std::setw(20) << "Re:SzSz/dist" << "\t"; - outfile5 << std::setw(20) << "Im:SzSz/dist" << "\n"; - - Vec3i lattice_dimensions = ::lattice->get_lattice_dimensions(); - int num_distance = std::ceil(0.5 * lattice_dimensions[0] / delta_r_); - //size of the array: (L/2)/delta_r_ - - vector> r_dependence_SxSy(num_distance, 0.0); - vector> r_dependence_SzSz(num_distance, 0.0); - vector count_r_dep(r_dependence_SxSy.size(),0); + if (output_real_space_ss_corr_) { + using namespace std; + using namespace globals; + + const double delta_t = periodogram_props_.sample_time; + + string name3 = seedname + "_S+S-_r_dependence.tsv"; + ofstream outfile3(name3.c_str()); + outfile3.setf(std::ios::right); + + string name4 = seedname + "_radial_dist.tsv"; + ofstream outfile4(name4.c_str()); + outfile4.setf(std::ios::right); + + string name5 = seedname + "_SzSz_r_dependence.tsv"; + ofstream outfile5(name5.c_str()); + outfile5.setf(std::ios::right); + + // header for the r-dependence of the S+S- file + outfile3 << std::setw(20) << "time_diff" << "\t";//(\t=tab) + outfile3 << std::setw(20) << "r_ij" << "\t";//(\t=tab) + outfile3 << std::setw(20) << "Re:S+S-/dist" << "\t"; + outfile3 << std::setw(20) << "Im:S+S-/dist" << "\n"; + + // header for the radial dist file + outfile4 << std::setw(20) << "bin" << "\t"; + outfile4 << std::setw(20) << "dist" << "\n"; + + // header for the r-dependence of the SzSz file + outfile5 << std::setw(20) << "time_diff" << "\t";//(\t=tab) + outfile5 << std::setw(20) << "r_ij" << "\t";//(\t=tab) + outfile5 << std::setw(20) << "Re:SzSz/dist" << "\t"; + outfile5 << std::setw(20) << "Im:SzSz/dist" << "\n"; + + Vec3i lattice_dimensions = ::lattice->get_lattice_dimensions(); + int num_distance = std::ceil(0.5 * lattice_dimensions[0] / delta_r_); + //size of the array: (L/2)/delta_r_ + + vector> r_dependence_SxSy(num_distance, 0.0); + vector> r_dependence_SzSz(num_distance, 0.0); + vector count_r_dep(r_dependence_SxSy.size(),0); + + auto radial_dist = radial_distribution_function(); + + for (unsigned j = 0; j < num_spins; ++j) { + const auto r_j = lattice->atom_position(j); + for (unsigned i = 0; i < num_spins; ++i) { + if(i == j){ + continue; + } + const auto r_i = lattice->atom_position(i); + const auto r_ij = lattice->displacement(r_i, r_j); // using j,i in this order gives r_j - r_i + + unsigned m = std::ceil(this->distance(r_ij)/delta_r_); + if(m < r_dependence_SxSy.size()) { + auto correlation = time_correlation(i, j); + + r_dependence_SxSy[m] += correlation[t_diff_r_dep_]; + auto num_time_samples_z = (unsigned int) (periodogram_props_.length); + for (auto t_ref = 0; t_ref < num_time_samples_z; t_ref++) { + r_dependence_SzSz[m] += spin_z[t_ref](i) * spin_z[t_ref + t_diff_r_dep_](j); + } + count_r_dep[m]++; + } + }//i + }//j - auto radial_dist = radial_distribution_function(); + for (auto m = 0; m < radial_dist.size(); m++) { + outfile4 << std::setw(20) << m << "\t"; + outfile4 << std::setw(20) << radial_dist[m] << "\n"; + } - for (unsigned j = 0; j < num_spins; ++j) { - const auto r_j = lattice->atom_position(j); - for (unsigned i = 0; i < num_spins; ++i) { - if(i == j){ - continue; + for (auto m = 1; m < r_dependence_SxSy.size(); m++) { + outfile3 << std::setw(20) << t_diff_r_dep_ * delta_t << "\t"; + outfile3 << std::setw(20) << m << "\t"; + if (count_r_dep[m] != 0) { + outfile3 << std::setw(20) << r_dependence_SxSy[m].real() / count_r_dep[m] << "\t"; + outfile3 << std::setw(20) << r_dependence_SxSy[m].imag() / count_r_dep[m] << "\n"; } - const auto r_i = lattice->atom_position(i); - const auto r_ij = lattice->displacement(r_i, r_j); // using j,i in this order gives r_j - r_i - - unsigned m = std::ceil(this->distance(r_ij)/delta_r_); - if(m < r_dependence_SxSy.size()) { - auto correlation = time_correlation(i, j); - - r_dependence_SxSy[m] += correlation[t_diff_r_dep_]; - auto num_time_samples_z = (unsigned int) (periodogram_props_.length); - for (auto t_ref = 0; t_ref < num_time_samples_z; t_ref++) { - r_dependence_SzSz[m] += spin_z[t_ref](i) * spin_z[t_ref + t_diff_r_dep_](j); - } - count_r_dep[m]++; + else{ + outfile3 << std::setw(20) << 0.0 << "\t"; + outfile3 << std::setw(20) << 0.0 << "\n"; } - }//i - }//j - - for (auto m = 0; m < radial_dist.size(); m++) { - outfile4 << std::setw(20) << m << "\t"; - outfile4 << std::setw(20) << radial_dist[m] << "\n"; - } - - for (auto m = 1; m < r_dependence_SxSy.size(); m++) { - outfile3 << std::setw(20) << t_diff_r_dep_ * delta_t << "\t"; - outfile3 << std::setw(20) << m << "\t"; - if (count_r_dep[m] != 0) { - outfile3 << std::setw(20) << r_dependence_SxSy[m].real() / count_r_dep[m] << "\t"; - outfile3 << std::setw(20) << r_dependence_SxSy[m].imag() / count_r_dep[m] << "\n"; } - else{ - outfile3 << std::setw(20) << 0.0 << "\t"; - outfile3 << std::setw(20) << 0.0 << "\n"; - } - } - for (auto m = 1; m < r_dependence_SzSz.size(); m++) { - outfile5 << std::setw(20) << t_diff_r_dep_ * delta_t << "\t"; - outfile5 << std::setw(20) << m << "\t"; - if (count_r_dep[m] != 0) { - outfile5 << std::setw(20) << r_dependence_SzSz[m].real() / count_r_dep[m] << "\t"; - outfile5 << std::setw(20) << r_dependence_SzSz[m].imag() / count_r_dep[m] << "\n"; - } - else{ - outfile5 << std::setw(20) << 0.0 << "\t"; - outfile5 << std::setw(20) << 0.0 << "\n"; + for (auto m = 1; m < r_dependence_SzSz.size(); m++) { + outfile5 << std::setw(20) << t_diff_r_dep_ * delta_t << "\t"; + outfile5 << std::setw(20) << m << "\t"; + if (count_r_dep[m] != 0) { + outfile5 << std::setw(20) << r_dependence_SzSz[m].real() / count_r_dep[m] << "\t"; + outfile5 << std::setw(20) << r_dependence_SzSz[m].imag() / count_r_dep[m] << "\n"; + } + else{ + outfile5 << std::setw(20) << 0.0 << "\t"; + outfile5 << std::setw(20) << 0.0 << "\n"; + } } - } - - outfile3.close(); - outfile4.close(); - outfile5.close(); + outfile3.close(); + outfile4.close(); + outfile5.close(); + } } double NeutronScatteringNoLatticeMonitor::distance(const Vec3 &r_ij) { diff --git a/src/jams/monitors/neutron_scattering_no_lattice.h b/src/jams/monitors/neutron_scattering_no_lattice.h index 67f85f4d..ae589210 100644 --- a/src/jams/monitors/neutron_scattering_no_lattice.h +++ b/src/jams/monitors/neutron_scattering_no_lattice.h @@ -36,6 +36,7 @@ class NeutronScatteringNoLatticeMonitor : public Monitor { jams::MultiArray calculate_polarized_cross_sections(const jams::MultiArray& spectrum, const jams::MultiArray& elastic_spectrum, const std::vector& polarizations); bool do_rspace_windowing_ = false; + bool output_real_space_ss_corr_ = false; jams::MultiArray rspace_displacement_; jams::MultiArray kspace_path_; jams::MultiArray kspace_spins_timeseries_; From 0b2131baf5a4e737df361bf13ad6c38e9407bcbb Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Thu, 23 Apr 2020 18:06:00 +0900 Subject: [PATCH 13/80] Explicitly output if we do the rspace_windowing and calculating the real-space spin-spin correlations were changed to be optional 2 --- src/jams/interface/config.h | 1 + 1 file changed, 1 insertion(+) diff --git a/src/jams/interface/config.h b/src/jams/interface/config.h index 1dc4a43f..a8a28137 100644 --- a/src/jams/interface/config.h +++ b/src/jams/interface/config.h @@ -22,6 +22,7 @@ namespace jams { if (setting.exists(name)) { return config_required(setting, name); } else { + std::cout << "config_optional is called but default value will be used. name = " << name << std::endl; return def; } } From 7bd69792aa6ef75d4ae88d6b2ae287022805c631 Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Thu, 23 Apr 2020 21:58:01 +0900 Subject: [PATCH 14/80] Explicitly output if we do the rspace_windowing and calculating the real-space spin-spin correlations were changed to be optional 3 --- src/jams/interface/config.h | 1 + 1 file changed, 1 insertion(+) diff --git a/src/jams/interface/config.h b/src/jams/interface/config.h index a8a28137..ee31c1a1 100644 --- a/src/jams/interface/config.h +++ b/src/jams/interface/config.h @@ -9,6 +9,7 @@ #include #include #include "jams/core/types.h" +#include void config_patch(libconfig::Setting& orig, const libconfig::Setting& patch); From 1268328d6ef33a33d7f8648fa7d7d1e092dd4dda Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Thu, 23 Apr 2020 22:12:52 +0900 Subject: [PATCH 15/80] Explicitly output if we do the rspace_windowing and calculating the real-space spin-spin correlations were changed to be optional 4 --- src/jams/monitors/neutron_scattering_no_lattice.cc | 2 +- src/jams/monitors/neutron_scattering_no_lattice.h | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/jams/monitors/neutron_scattering_no_lattice.cc b/src/jams/monitors/neutron_scattering_no_lattice.cc index 812ac654..2bb123fa 100644 --- a/src/jams/monitors/neutron_scattering_no_lattice.cc +++ b/src/jams/monitors/neutron_scattering_no_lattice.cc @@ -333,7 +333,7 @@ void NeutronScatteringNoLatticeMonitor::store_kspace_data_on_path() { const double a0 = 0.40217, a1 = 0.49704, a2 = 0.09392, a3 = 0.00183; const double x = (kTwoPi * norm(r)); window = a0 + a1 * cos(x) + a2 * cos(2 * x) + a3 * cos(3 * x); - cout << "do_rspace_windowing_ = true, value = " << do_rspace_windowing_ << endl; +// cout << "do_rspace_windowing_ = true, value = " << do_rspace_windowing_ << endl; } //Fourier transform auto f0 = exp(-kImagTwoPi * dot(delta_q, r)); diff --git a/src/jams/monitors/neutron_scattering_no_lattice.h b/src/jams/monitors/neutron_scattering_no_lattice.h index ae589210..4ba11bd7 100644 --- a/src/jams/monitors/neutron_scattering_no_lattice.h +++ b/src/jams/monitors/neutron_scattering_no_lattice.h @@ -35,7 +35,7 @@ class NeutronScatteringNoLatticeMonitor : public Monitor { jams::MultiArray calculate_unpolarized_cross_section(const jams::MultiArray& spectrum, const jams::MultiArray& elastic_spectrum); jams::MultiArray calculate_polarized_cross_sections(const jams::MultiArray& spectrum, const jams::MultiArray& elastic_spectrum, const std::vector& polarizations); - bool do_rspace_windowing_ = false; + bool do_rspace_windowing_ = true; bool output_real_space_ss_corr_ = false; jams::MultiArray rspace_displacement_; jams::MultiArray kspace_path_; From 04838d6f4b86df606e5ebc72feee4bec46aac80e Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Fri, 24 Apr 2020 08:58:30 +0900 Subject: [PATCH 16/80] rspace_window=true & rspace_displacement_(i) is measured compared to origin. --- src/jams/monitors/neutron_scattering_no_lattice.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/jams/monitors/neutron_scattering_no_lattice.cc b/src/jams/monitors/neutron_scattering_no_lattice.cc index 2bb123fa..94f563ab 100644 --- a/src/jams/monitors/neutron_scattering_no_lattice.cc +++ b/src/jams/monitors/neutron_scattering_no_lattice.cc @@ -117,9 +117,9 @@ void NeutronScatteringNoLatticeMonitor::configure_kspace_vectors(const libconfig Vec3i lattice_dimensions = ::lattice->get_lattice_dimensions(); for (auto i = 0; i < globals::s.size(0); ++i) { //rspace_displacement_ is the position vector of i-th spin - //rspace_displacement_(i) = lattice->displacement({0,0,0}, lattice->atom_position(i)); + rspace_displacement_(i) = lattice->displacement({0,0,0}, lattice->atom_position(i)); //Generalize so that we can impose open boundary - rspace_displacement_(i) = lattice->displacement({lattice_dimensions[0]*0.5,lattice_dimensions[1]*0.5,lattice_dimensions[2]*0.5}, lattice->atom_position(i)); + //rspace_displacement_(i) = lattice->displacement({lattice_dimensions[0]*0.5,lattice_dimensions[1]*0.5,lattice_dimensions[2]*0.5}, lattice->atom_position(i)); } } From 34fc3967d0dcf50a62fc1a393f25f1b72a388fd4 Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Fri, 24 Apr 2020 09:05:37 +0900 Subject: [PATCH 17/80] rspace_window=false & rspace_displacement_(i) is measured compared to the centre. --- src/jams/monitors/neutron_scattering_no_lattice.cc | 4 ++-- src/jams/monitors/neutron_scattering_no_lattice.h | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/jams/monitors/neutron_scattering_no_lattice.cc b/src/jams/monitors/neutron_scattering_no_lattice.cc index 94f563ab..2bb123fa 100644 --- a/src/jams/monitors/neutron_scattering_no_lattice.cc +++ b/src/jams/monitors/neutron_scattering_no_lattice.cc @@ -117,9 +117,9 @@ void NeutronScatteringNoLatticeMonitor::configure_kspace_vectors(const libconfig Vec3i lattice_dimensions = ::lattice->get_lattice_dimensions(); for (auto i = 0; i < globals::s.size(0); ++i) { //rspace_displacement_ is the position vector of i-th spin - rspace_displacement_(i) = lattice->displacement({0,0,0}, lattice->atom_position(i)); + //rspace_displacement_(i) = lattice->displacement({0,0,0}, lattice->atom_position(i)); //Generalize so that we can impose open boundary - //rspace_displacement_(i) = lattice->displacement({lattice_dimensions[0]*0.5,lattice_dimensions[1]*0.5,lattice_dimensions[2]*0.5}, lattice->atom_position(i)); + rspace_displacement_(i) = lattice->displacement({lattice_dimensions[0]*0.5,lattice_dimensions[1]*0.5,lattice_dimensions[2]*0.5}, lattice->atom_position(i)); } } diff --git a/src/jams/monitors/neutron_scattering_no_lattice.h b/src/jams/monitors/neutron_scattering_no_lattice.h index 4ba11bd7..ae589210 100644 --- a/src/jams/monitors/neutron_scattering_no_lattice.h +++ b/src/jams/monitors/neutron_scattering_no_lattice.h @@ -35,7 +35,7 @@ class NeutronScatteringNoLatticeMonitor : public Monitor { jams::MultiArray calculate_unpolarized_cross_section(const jams::MultiArray& spectrum, const jams::MultiArray& elastic_spectrum); jams::MultiArray calculate_polarized_cross_sections(const jams::MultiArray& spectrum, const jams::MultiArray& elastic_spectrum, const std::vector& polarizations); - bool do_rspace_windowing_ = true; + bool do_rspace_windowing_ = false; bool output_real_space_ss_corr_ = false; jams::MultiArray rspace_displacement_; jams::MultiArray kspace_path_; From e5ffe9269768e14b6d08ce6a1fc0ad9a515ebb02 Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Fri, 24 Apr 2020 17:43:30 +0900 Subject: [PATCH 18/80] Fixed "do_rspace_windowing_" and "output_real_space_ss_corr_" input. --- src/jams/helpers/neutrons.cc | 9 +++++++++ src/jams/interface/config.h | 12 ++++++------ src/jams/monitors/neutron_scattering_no_lattice.cc | 6 +++--- src/jams/monitors/neutron_scattering_no_lattice.h | 4 ++-- 4 files changed, 20 insertions(+), 11 deletions(-) diff --git a/src/jams/helpers/neutrons.cc b/src/jams/helpers/neutrons.cc index 6fe35bf3..cae33778 100644 --- a/src/jams/helpers/neutrons.cc +++ b/src/jams/helpers/neutrons.cc @@ -16,6 +16,15 @@ namespace jams { double p = (l == 0) ? 1.0 : s_sq; //if l=0 then insert 1.0 to p, otherwise s_sq ffq += g[l] * p * (j[l].A * exp(-j[l].a * s_sq) + j[l].B * exp(-j[l].b * s_sq) + j[l].C * exp(-j[l].c * s_sq) + j[l].D); + +// cout << "j_params[" << l << "].A = " << j[l].A << endl; +// cout << "j_params[" << l << "].a = " << j[l].a << endl; +// cout << "j_params[" << l << "].B = " << j[l].B << endl; +// cout << "j_params[" << l << "].b = " << j[l].b << endl; +// cout << "j_params[" << l << "].C = " << j[l].C << endl; +// cout << "j_params[" << l << "].c = " << j[l].c << endl; +// cout << "j_params[" << l << "].D = " << j[l].D << endl; +// cout << "g_params[" << l << "] = " << g[l] << endl; } return 0.5 * ffq; } diff --git a/src/jams/interface/config.h b/src/jams/interface/config.h index ee31c1a1..e0f68e03 100644 --- a/src/jams/interface/config.h +++ b/src/jams/interface/config.h @@ -23,34 +23,34 @@ namespace jams { if (setting.exists(name)) { return config_required(setting, name); } else { - std::cout << "config_optional is called but default value will be used. name = " << name << std::endl; +// std::cout << "config_optional is called but default value will be used. name = " << name << std::endl; return def; } } template<> inline std::string config_required(const libconfig::Setting &setting, const std::string &name) { - return setting[name].c_str(); + return setting[name].c_str(); } template<> inline bool config_required(const libconfig::Setting &setting, const std::string &name) { - return bool(setting[name]); + return bool(setting[name]); } template<> inline int config_required(const libconfig::Setting &setting, const std::string &name) { - return int(setting[name]); + return int(setting[name]); } template<> inline long config_required(const libconfig::Setting &setting, const std::string &name) { - return long(setting[name]); + return long(setting[name]); } template<> inline unsigned config_required(const libconfig::Setting &setting, const std::string &name) { - return unsigned(setting[name]); + return unsigned(setting[name]); } template<> diff --git a/src/jams/monitors/neutron_scattering_no_lattice.cc b/src/jams/monitors/neutron_scattering_no_lattice.cc index 2bb123fa..25f28141 100644 --- a/src/jams/monitors/neutron_scattering_no_lattice.cc +++ b/src/jams/monitors/neutron_scattering_no_lattice.cc @@ -22,11 +22,11 @@ NeutronScatteringNoLatticeMonitor::NeutronScatteringNoLatticeMonitor(const libco configure_kspace_vectors(settings); - config_optional(settings, "rspace_windowing", do_rspace_windowing_); + do_rspace_windowing_ = config_optional(settings, "rspace_windowing", false); cout << "do_rspace_windowing_ = " << do_rspace_windowing_ << endl; - config_optional(settings, "output_real_space_ss_corr_", output_real_space_ss_corr_); - cout << "output_real_space_ss_corr_ = " << output_real_space_ss_corr_ << endl; + output_real_space_ss_corr_ = config_optional(settings, "output_real_space_ss_corr_", false); + cout << "output_real_space_ss_corr_ = " << output_real_space_ss_corr_ << endl; if (settings.exists("form_factor")) { configure_form_factors(settings["form_factor"]); diff --git a/src/jams/monitors/neutron_scattering_no_lattice.h b/src/jams/monitors/neutron_scattering_no_lattice.h index ae589210..3fdf2854 100644 --- a/src/jams/monitors/neutron_scattering_no_lattice.h +++ b/src/jams/monitors/neutron_scattering_no_lattice.h @@ -35,8 +35,8 @@ class NeutronScatteringNoLatticeMonitor : public Monitor { jams::MultiArray calculate_unpolarized_cross_section(const jams::MultiArray& spectrum, const jams::MultiArray& elastic_spectrum); jams::MultiArray calculate_polarized_cross_sections(const jams::MultiArray& spectrum, const jams::MultiArray& elastic_spectrum, const std::vector& polarizations); - bool do_rspace_windowing_ = false; - bool output_real_space_ss_corr_ = false; + bool do_rspace_windowing_; + bool output_real_space_ss_corr_; jams::MultiArray rspace_displacement_; jams::MultiArray kspace_path_; jams::MultiArray kspace_spins_timeseries_; From f28247435b7e5f5f523c4e68ad0d27568a0ebde2 Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Sun, 7 Jun 2020 22:18:39 +0900 Subject: [PATCH 19/80] Deleted unused time correlation of spins. --- .../monitors/neutron_scattering_no_lattice.cc | 135 ++++++------------ .../monitors/neutron_scattering_no_lattice.h | 5 +- 2 files changed, 44 insertions(+), 96 deletions(-) diff --git a/src/jams/monitors/neutron_scattering_no_lattice.cc b/src/jams/monitors/neutron_scattering_no_lattice.cc index 25f28141..2c963982 100644 --- a/src/jams/monitors/neutron_scattering_no_lattice.cc +++ b/src/jams/monitors/neutron_scattering_no_lattice.cc @@ -25,8 +25,8 @@ NeutronScatteringNoLatticeMonitor::NeutronScatteringNoLatticeMonitor(const libco do_rspace_windowing_ = config_optional(settings, "rspace_windowing", false); cout << "do_rspace_windowing_ = " << do_rspace_windowing_ << endl; - output_real_space_ss_corr_ = config_optional(settings, "output_real_space_ss_corr_", false); - cout << "output_real_space_ss_corr_ = " << output_real_space_ss_corr_ << endl; + output_dist_ = config_optional(settings, "output_dist", false); + cout << "output_dist_ = " << output_dist_ << endl; if (settings.exists("form_factor")) { configure_form_factors(settings["form_factor"]); @@ -42,10 +42,6 @@ NeutronScatteringNoLatticeMonitor::NeutronScatteringNoLatticeMonitor(const libco configure_periodogram(settings["periodogram"]); } - if (output_real_space_ss_corr_) { - t_diff_r_dep_ = config_required(settings, "t_diff_r_dep_"); - } - periodogram_props_.sample_time = output_step_freq_ * solver->time_step(); //zero: set all the elements as 0 @@ -61,7 +57,7 @@ void NeutronScatteringNoLatticeMonitor::update(Solver *solver) { store_kspace_data_on_path(); periodogram_index_++; - if (output_real_space_ss_corr_) { + if (output_dist_) { using namespace globals; jams::MultiArray spin_x_time(num_spins); @@ -408,39 +404,27 @@ jams::MultiArray NeutronScatteringNoLatticeMonitor::average_kspace_ti } NeutronScatteringNoLatticeMonitor::~NeutronScatteringNoLatticeMonitor() { - if (output_real_space_ss_corr_) { + if (output_dist_) { using namespace std; using namespace globals; const double delta_t = periodogram_props_.sample_time; - string name3 = seedname + "_S+S-_r_dependence.tsv"; - ofstream outfile3(name3.c_str()); - outfile3.setf(std::ios::right); - string name4 = seedname + "_radial_dist.tsv"; ofstream outfile4(name4.c_str()); outfile4.setf(std::ios::right); - string name5 = seedname + "_SzSz_r_dependence.tsv"; + string name5 = seedname + "_radial_dist_z.tsv"; ofstream outfile5(name5.c_str()); outfile5.setf(std::ios::right); - // header for the r-dependence of the S+S- file - outfile3 << std::setw(20) << "time_diff" << "\t";//(\t=tab) - outfile3 << std::setw(20) << "r_ij" << "\t";//(\t=tab) - outfile3 << std::setw(20) << "Re:S+S-/dist" << "\t"; - outfile3 << std::setw(20) << "Im:S+S-/dist" << "\n"; - // header for the radial dist file outfile4 << std::setw(20) << "bin" << "\t"; outfile4 << std::setw(20) << "dist" << "\n"; - // header for the r-dependence of the SzSz file - outfile5 << std::setw(20) << "time_diff" << "\t";//(\t=tab) - outfile5 << std::setw(20) << "r_ij" << "\t";//(\t=tab) - outfile5 << std::setw(20) << "Re:SzSz/dist" << "\t"; - outfile5 << std::setw(20) << "Im:SzSz/dist" << "\n"; + // header for the radial dist (z) file + outfile5 << std::setw(20) << "bin" << "\t"; + outfile5 << std::setw(20) << "dist" << "\n"; Vec3i lattice_dimensions = ::lattice->get_lattice_dimensions(); int num_distance = std::ceil(0.5 * lattice_dimensions[0] / delta_r_); @@ -451,62 +435,18 @@ NeutronScatteringNoLatticeMonitor::~NeutronScatteringNoLatticeMonitor() { vector count_r_dep(r_dependence_SxSy.size(),0); auto radial_dist = radial_distribution_function(); - - for (unsigned j = 0; j < num_spins; ++j) { - const auto r_j = lattice->atom_position(j); - for (unsigned i = 0; i < num_spins; ++i) { - if(i == j){ - continue; - } - const auto r_i = lattice->atom_position(i); - const auto r_ij = lattice->displacement(r_i, r_j); // using j,i in this order gives r_j - r_i - - unsigned m = std::ceil(this->distance(r_ij)/delta_r_); - if(m < r_dependence_SxSy.size()) { - auto correlation = time_correlation(i, j); - - r_dependence_SxSy[m] += correlation[t_diff_r_dep_]; - auto num_time_samples_z = (unsigned int) (periodogram_props_.length); - for (auto t_ref = 0; t_ref < num_time_samples_z; t_ref++) { - r_dependence_SzSz[m] += spin_z[t_ref](i) * spin_z[t_ref + t_diff_r_dep_](j); - } - count_r_dep[m]++; - } - }//i - }//j + auto radial_dist_z = radial_distribution_function_z(); for (auto m = 0; m < radial_dist.size(); m++) { outfile4 << std::setw(20) << m << "\t"; outfile4 << std::setw(20) << radial_dist[m] << "\n"; } - for (auto m = 1; m < r_dependence_SxSy.size(); m++) { - outfile3 << std::setw(20) << t_diff_r_dep_ * delta_t << "\t"; - outfile3 << std::setw(20) << m << "\t"; - if (count_r_dep[m] != 0) { - outfile3 << std::setw(20) << r_dependence_SxSy[m].real() / count_r_dep[m] << "\t"; - outfile3 << std::setw(20) << r_dependence_SxSy[m].imag() / count_r_dep[m] << "\n"; - } - else{ - outfile3 << std::setw(20) << 0.0 << "\t"; - outfile3 << std::setw(20) << 0.0 << "\n"; - } - } - - for (auto m = 1; m < r_dependence_SzSz.size(); m++) { - outfile5 << std::setw(20) << t_diff_r_dep_ * delta_t << "\t"; + for (auto m = 0; m < radial_dist.size(); m++) { outfile5 << std::setw(20) << m << "\t"; - if (count_r_dep[m] != 0) { - outfile5 << std::setw(20) << r_dependence_SzSz[m].real() / count_r_dep[m] << "\t"; - outfile5 << std::setw(20) << r_dependence_SzSz[m].imag() / count_r_dep[m] << "\n"; - } - else{ - outfile5 << std::setw(20) << 0.0 << "\t"; - outfile5 << std::setw(20) << 0.0 << "\n"; - } + outfile5 << std::setw(20) << radial_dist[m] << "\n"; } - outfile3.close(); outfile4.close(); outfile5.close(); } @@ -516,29 +456,6 @@ double NeutronScatteringNoLatticeMonitor::distance(const Vec3 &r_ij) { return std::sqrt(r_ij[0] * r_ij[0] + r_ij[1] * r_ij[1] + r_ij[2] * r_ij[2]); } -std::vector> NeutronScatteringNoLatticeMonitor::time_correlation(unsigned i, unsigned j) { - const auto num_time_samples = (unsigned int)(periodogram_props_.length); - std::vector> correlation(num_time_samples, 0.0); - - for (auto t_diff = 0; t_diff < correlation.size(); t_diff++) { - for (auto t_ref = 0; t_ref < num_time_samples; t_ref++) { - if ((t_ref + t_diff) < num_time_samples && t_ref < num_time_samples) { - correlation[t_diff] += - spin_x[t_ref](i) * spin_x[t_ref + t_diff](j) + spin_y[t_ref](i) * spin_y[t_ref + t_diff](j) - + kImagOne * (-spin_x[t_ref](i) * spin_y[t_ref + t_diff](j) + spin_y[t_ref](i) * spin_x[t_ref + t_diff](j)); - -// if(t_diff == 0){ -// std::cout << "num_t_samples = " << num_time_samples << std::endl; -// std::cout << "t_diff = " << t_diff << ", t_ref = " << t_ref << std::endl; -// std::cout << "corr (Re) = " << correlation[t_diff].real() << ", corr (Im) = " << correlation[t_diff].imag() << std::endl; -// } - - } - } - } - return correlation; -} - std::vector NeutronScatteringNoLatticeMonitor::radial_distribution_function() { using namespace globals; using namespace std; @@ -563,6 +480,36 @@ std::vector NeutronScatteringNoLatticeMonitor::radial_distribution_funct } for (auto kk = 0; kk < dist.size(); kk++){ dist[kk] /= num_spins; + cout << "dist[" << kk << "] = " << dist[kk] << ", dist/num_spins = " << dist[kk]/num_spins << endl; } return dist; +} + +std::vector NeutronScatteringNoLatticeMonitor::radial_distribution_function_z() { + using namespace globals; + using namespace std; + Vec3i lattice_dimensions = ::lattice->get_lattice_dimensions(); + //Here we assume lattice_dimensions are all the same (x,y,z) + std::vector dist_z(std::ceil(0.5 * lattice_dimensions[0] / delta_r_),0.0); + for(unsigned j = 0; j < num_spins; ++j){ + const auto r_j = lattice->atom_position(j); + for(unsigned i = 0; i < num_spins; ++i){ + if(i == j){ + dist_z[0] = 0; + continue; + } + const auto r_i = lattice->atom_position(i); + const auto r_ij = (lattice->displacement(r_i, r_j)); // using j,i in this order gives r_j - r_i + int m = std::ceil(r_ij[2]/delta_r_); + if(m > dist_z.size()){ + continue; + } + dist_z[m]++; + } + } + for (auto kk = 0; kk < dist_z.size(); kk++){ + dist_z[kk] /= num_spins; + cout << "dist_z[" << kk << "] = " << dist_z[kk] << ", dist_z/num_spins = " << dist_z[kk]/num_spins << endl; + } + return dist_z; } \ No newline at end of file diff --git a/src/jams/monitors/neutron_scattering_no_lattice.h b/src/jams/monitors/neutron_scattering_no_lattice.h index 3fdf2854..be8d1615 100644 --- a/src/jams/monitors/neutron_scattering_no_lattice.h +++ b/src/jams/monitors/neutron_scattering_no_lattice.h @@ -36,7 +36,7 @@ class NeutronScatteringNoLatticeMonitor : public Monitor { jams::MultiArray calculate_polarized_cross_sections(const jams::MultiArray& spectrum, const jams::MultiArray& elastic_spectrum, const std::vector& polarizations); bool do_rspace_windowing_; - bool output_real_space_ss_corr_; + bool output_dist_; jams::MultiArray rspace_displacement_; jams::MultiArray kspace_path_; jams::MultiArray kspace_spins_timeseries_; @@ -55,7 +55,7 @@ class NeutronScatteringNoLatticeMonitor : public Monitor { int num_k_; Vec3 kvector_; bool evaluate_form_factor_ = false; - double delta_r_ = 0.001; + double delta_r_ = 0.01; int t_diff_r_dep_; std::vector> spin_x; @@ -63,6 +63,7 @@ class NeutronScatteringNoLatticeMonitor : public Monitor { std::vector> spin_z; std::vector radial_distribution_function(); + std::vector radial_distribution_function_z(); double distance(const Vec3 &r_ij); std::vector> time_correlation(unsigned i, unsigned j); }; From 645ea0a3bf279d63b3c10b1c9afeb7b023e1bc26 Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Sun, 7 Jun 2020 22:19:16 +0900 Subject: [PATCH 20/80] Added function calculating crystal limit spectrum. --- src/jams/hamiltonian/exchange_functional.cc | 40 +++++++++++++++++++++ 1 file changed, 40 insertions(+) diff --git a/src/jams/hamiltonian/exchange_functional.cc b/src/jams/hamiltonian/exchange_functional.cc index 74b51f65..b95f49a3 100644 --- a/src/jams/hamiltonian/exchange_functional.cc +++ b/src/jams/hamiltonian/exchange_functional.cc @@ -13,6 +13,7 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se radius_cutoff_ = jams::config_required(settings, "r_cutoff"); cout << " cutoff radius: " << jams::fmt::decimal << radius_cutoff_ << "\n"; + cout << " max cutoff radius: " << lattice->max_interaction_radius() << "\n"; if (radius_cutoff_ > lattice->max_interaction_radius()) { throw std::runtime_error("cutoff radius is larger than the maximum radius which avoids self interaction"); @@ -21,8 +22,29 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se ofstream of(seedname + "_exchange_functional.tsv"); output_exchange_functional(of, exchange_functional, radius_cutoff_); + string name3 = seedname + "_spectrum_crystal_limit.tsv"; + ofstream outfile3(name3.c_str()); + outfile3.setf(std::ios::right); + + // header for crystal-limit spectrum file + outfile3 << std::setw(20) << "k" << "\t";//(\t=tab) + outfile3 << std::setw(20) << "Re:E(k)" << "\t"; + outfile3 << std::setw(20) << "Im:E(k)" << "\n"; + auto counter = 0; vector nbrs; + // --- for crystal limit spectrum --- + int num_k = jams::config_required(settings, "num_k"); + std::vector> spectrum_crystal_limit(num_k,{0.0,0.0}); + double kmax = jams::config_required(settings, "kmax"); + Vec3 kvector = jams::config_required(settings, "kvector"); + jams::MultiArray k; + k.resize(num_k+1); + for (auto n = 0; n < k.size(); ++n) { + k(n) = kvector * n * (kmax / num_k); +// cout << "n = " << n << ", kspace_path_(n) = " << k(n) << endl; + } + // --- for crystal limit spectrum --- for (auto i = 0; i < globals::num_spins; ++i) { nbrs.clear(); ::lattice->atom_neighbours(i, radius_cutoff_, nbrs); @@ -33,13 +55,31 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se continue; } const auto rij = norm(lattice->displacement(i, j)); + // --- for crystal limit spectrum --- + const auto rij_vec = lattice->displacement(i, j); this->insert_interaction_scalar(i, j, input_unit_conversion_ * exchange_functional(rij)); counter++; + for (auto kk = 0; kk < spectrum_crystal_limit.size(); kk++){ + double kr = std::inner_product(k(kk).begin(), k(kk).end(), rij_vec.begin(), 0.0); + std::complex tmp = {input_unit_conversion_ * exchange_functional(rij)* (1.0-cos(kr)), input_unit_conversion_ * exchange_functional(rij) * sin(kr)}; + spectrum_crystal_limit[kk] += tmp; + } + // --- for crystal limit spectrum --- } } + for (auto kk = 0; kk < spectrum_crystal_limit.size(); kk++) { + spectrum_crystal_limit[kk] /= globals::num_spins; + } cout << " total interactions " << jams::fmt::integer << counter << "\n"; cout << " average interactions per spin " << jams::fmt::decimal << counter / double(globals::num_spins) << "\n"; + // --- for crystal limit spectrum --- + for (auto m = 0; m < spectrum_crystal_limit.size(); m++) { + cout << " spectrum_crystal_limit (" << m << ") = " << spectrum_crystal_limit[m] << "\n"; + outfile3 << std::setw(20) << k(m) << "\t"; + outfile3 << std::setw(20) << spectrum_crystal_limit[m] << "\n"; + } + // --- for crystal limit spectrum --- finalize(); } From 6f71f3df87bba84b719a9aaea313d1f6c50a2326 Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Mon, 8 Jun 2020 12:24:27 +0900 Subject: [PATCH 21/80] [Fixed] Added function calculating crystal limit spectrum. --- src/jams/hamiltonian/exchange_functional.cc | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/jams/hamiltonian/exchange_functional.cc b/src/jams/hamiltonian/exchange_functional.cc index b95f49a3..e9c343a6 100644 --- a/src/jams/hamiltonian/exchange_functional.cc +++ b/src/jams/hamiltonian/exchange_functional.cc @@ -35,7 +35,7 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se vector nbrs; // --- for crystal limit spectrum --- int num_k = jams::config_required(settings, "num_k"); - std::vector> spectrum_crystal_limit(num_k,{0.0,0.0}); + std::vector> spectrum_crystal_limit(num_k+1,{0.0,0.0}); double kmax = jams::config_required(settings, "kmax"); Vec3 kvector = jams::config_required(settings, "kvector"); jams::MultiArray k; @@ -75,9 +75,12 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se cout << " average interactions per spin " << jams::fmt::decimal << counter / double(globals::num_spins) << "\n"; // --- for crystal limit spectrum --- for (auto m = 0; m < spectrum_crystal_limit.size(); m++) { - cout << " spectrum_crystal_limit (" << m << ") = " << spectrum_crystal_limit[m] << "\n"; +// cout << " spectrum_crystal_limit (" << m << ") = " << spectrum_crystal_limit[m] << "\n"; +// cout << " real (" << m << ") = " << spectrum_crystal_limit[m].real() << "\n"; +// cout << " imag (" << m << ") = " << spectrum_crystal_limit[m].imag() << "\n"; outfile3 << std::setw(20) << k(m) << "\t"; - outfile3 << std::setw(20) << spectrum_crystal_limit[m] << "\n"; + outfile3 << std::setw(20) << spectrum_crystal_limit[m].real() << "\n"; + outfile3 << std::setw(20) << spectrum_crystal_limit[m].imag() << "\n"; } // --- for crystal limit spectrum --- From 3479d2fb5b3e2886667d82dc959be223e0978244 Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Mon, 8 Jun 2020 12:30:24 +0900 Subject: [PATCH 22/80] [Fixed-2] Added function calculating crystal limit spectrum. --- src/jams/hamiltonian/exchange_functional.cc | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/src/jams/hamiltonian/exchange_functional.cc b/src/jams/hamiltonian/exchange_functional.cc index e9c343a6..6ea0f1cc 100644 --- a/src/jams/hamiltonian/exchange_functional.cc +++ b/src/jams/hamiltonian/exchange_functional.cc @@ -27,7 +27,9 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se outfile3.setf(std::ios::right); // header for crystal-limit spectrum file - outfile3 << std::setw(20) << "k" << "\t";//(\t=tab) + outfile3 << std::setw(20) << "kx" << "\t";//(\t=tab) + outfile3 << std::setw(20) << "ky" << "\t";//(\t=tab) + outfile3 << std::setw(20) << "kz" << "\t";//(\t=tab) outfile3 << std::setw(20) << "Re:E(k)" << "\t"; outfile3 << std::setw(20) << "Im:E(k)" << "\n"; @@ -78,8 +80,10 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se // cout << " spectrum_crystal_limit (" << m << ") = " << spectrum_crystal_limit[m] << "\n"; // cout << " real (" << m << ") = " << spectrum_crystal_limit[m].real() << "\n"; // cout << " imag (" << m << ") = " << spectrum_crystal_limit[m].imag() << "\n"; - outfile3 << std::setw(20) << k(m) << "\t"; - outfile3 << std::setw(20) << spectrum_crystal_limit[m].real() << "\n"; + outfile3 << std::setw(20) << k(m)[0] << "\t"; + outfile3 << std::setw(20) << k(m)[1] << "\t"; + outfile3 << std::setw(20) << k(m)[2] << "\t"; + outfile3 << std::setw(20) << spectrum_crystal_limit[m].real() << "\t"; outfile3 << std::setw(20) << spectrum_crystal_limit[m].imag() << "\n"; } // --- for crystal limit spectrum --- From c04d0933d641ff6a82ba3f0245d598654f16d67c Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Mon, 8 Jun 2020 13:10:05 +0900 Subject: [PATCH 23/80] [Fixed-3] Added function calculating crystal limit spectrum. --- src/jams/hamiltonian/exchange_functional.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/jams/hamiltonian/exchange_functional.cc b/src/jams/hamiltonian/exchange_functional.cc index 6ea0f1cc..602c767c 100644 --- a/src/jams/hamiltonian/exchange_functional.cc +++ b/src/jams/hamiltonian/exchange_functional.cc @@ -63,7 +63,7 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se counter++; for (auto kk = 0; kk < spectrum_crystal_limit.size(); kk++){ double kr = std::inner_product(k(kk).begin(), k(kk).end(), rij_vec.begin(), 0.0); - std::complex tmp = {input_unit_conversion_ * exchange_functional(rij)* (1.0-cos(kr)), input_unit_conversion_ * exchange_functional(rij) * sin(kr)}; + std::complex tmp = { exchange_functional(rij)* (1.0-cos(kr)), exchange_functional(rij) * sin(kr)}; spectrum_crystal_limit[kk] += tmp; } // --- for crystal limit spectrum --- From 25d4639ca96b71b61f445851a39db4aec35fec1b Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Mon, 8 Jun 2020 13:43:10 +0900 Subject: [PATCH 24/80] [Fixed-4] Added function calculating crystal limit spectrum. --- src/jams/hamiltonian/exchange_functional.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/jams/hamiltonian/exchange_functional.cc b/src/jams/hamiltonian/exchange_functional.cc index 602c767c..e558420a 100644 --- a/src/jams/hamiltonian/exchange_functional.cc +++ b/src/jams/hamiltonian/exchange_functional.cc @@ -63,7 +63,7 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se counter++; for (auto kk = 0; kk < spectrum_crystal_limit.size(); kk++){ double kr = std::inner_product(k(kk).begin(), k(kk).end(), rij_vec.begin(), 0.0); - std::complex tmp = { exchange_functional(rij)* (1.0-cos(kr)), exchange_functional(rij) * sin(kr)}; + std::complex tmp = { exchange_functional(rij)* (1.0-cos(kTwoPi*kr)), exchange_functional(rij) * sin(kTwoPi*kr)}; spectrum_crystal_limit[kk] += tmp; } // --- for crystal limit spectrum --- From 829a05b61bbbf4e56a3df8958a57316ef7b64dad Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Mon, 8 Jun 2020 13:54:47 +0900 Subject: [PATCH 25/80] [Fixed-5] Added function calculating crystal limit spectrum. --- src/jams/hamiltonian/exchange_functional.cc | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/jams/hamiltonian/exchange_functional.cc b/src/jams/hamiltonian/exchange_functional.cc index e558420a..1df17579 100644 --- a/src/jams/hamiltonian/exchange_functional.cc +++ b/src/jams/hamiltonian/exchange_functional.cc @@ -34,6 +34,7 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se outfile3 << std::setw(20) << "Im:E(k)" << "\n"; auto counter = 0; + auto counter2 = 0; vector nbrs; // --- for crystal limit spectrum --- int num_k = jams::config_required(settings, "num_k"); @@ -65,6 +66,9 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se double kr = std::inner_product(k(kk).begin(), k(kk).end(), rij_vec.begin(), 0.0); std::complex tmp = { exchange_functional(rij)* (1.0-cos(kTwoPi*kr)), exchange_functional(rij) * sin(kTwoPi*kr)}; spectrum_crystal_limit[kk] += tmp; + if(kr > 0.0){ + counter2++; + } } // --- for crystal limit spectrum --- } @@ -75,6 +79,7 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se cout << " total interactions " << jams::fmt::integer << counter << "\n"; cout << " average interactions per spin " << jams::fmt::decimal << counter / double(globals::num_spins) << "\n"; + cout << " average interactions per spin (kr > 0) " << jams::fmt::decimal << counter2 / double(globals::num_spins)/spectrum_crystal_limit.size() << "\n"; // --- for crystal limit spectrum --- for (auto m = 0; m < spectrum_crystal_limit.size(); m++) { // cout << " spectrum_crystal_limit (" << m << ") = " << spectrum_crystal_limit[m] << "\n"; From 5b0769db0b69468243db4bf245069243f1551474 Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Mon, 8 Jun 2020 14:04:34 +0900 Subject: [PATCH 26/80] [Fixed-6] Added function calculating crystal limit spectrum. --- src/jams/hamiltonian/exchange_functional.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/jams/hamiltonian/exchange_functional.cc b/src/jams/hamiltonian/exchange_functional.cc index 1df17579..35e12571 100644 --- a/src/jams/hamiltonian/exchange_functional.cc +++ b/src/jams/hamiltonian/exchange_functional.cc @@ -66,7 +66,7 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se double kr = std::inner_product(k(kk).begin(), k(kk).end(), rij_vec.begin(), 0.0); std::complex tmp = { exchange_functional(rij)* (1.0-cos(kTwoPi*kr)), exchange_functional(rij) * sin(kTwoPi*kr)}; spectrum_crystal_limit[kk] += tmp; - if(kr > 0.0){ + if(kr != 0.0){ counter2++; } } @@ -79,7 +79,7 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se cout << " total interactions " << jams::fmt::integer << counter << "\n"; cout << " average interactions per spin " << jams::fmt::decimal << counter / double(globals::num_spins) << "\n"; - cout << " average interactions per spin (kr > 0) " << jams::fmt::decimal << counter2 / double(globals::num_spins)/spectrum_crystal_limit.size() << "\n"; + cout << " average interactions per spin (kr != 0) " << jams::fmt::decimal << counter2 / double(globals::num_spins)/num_k << "\n"; // --- for crystal limit spectrum --- for (auto m = 0; m < spectrum_crystal_limit.size(); m++) { // cout << " spectrum_crystal_limit (" << m << ") = " << spectrum_crystal_limit[m] << "\n"; From 455514ef464a71c4ed9ff3a88f064f32b882e8a3 Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Mon, 8 Jun 2020 14:10:06 +0900 Subject: [PATCH 27/80] [Fixed-6] Added function calculating crystal limit spectrum. --- src/jams/hamiltonian/exchange_functional.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/jams/hamiltonian/exchange_functional.cc b/src/jams/hamiltonian/exchange_functional.cc index 35e12571..bf397572 100644 --- a/src/jams/hamiltonian/exchange_functional.cc +++ b/src/jams/hamiltonian/exchange_functional.cc @@ -64,9 +64,9 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se counter++; for (auto kk = 0; kk < spectrum_crystal_limit.size(); kk++){ double kr = std::inner_product(k(kk).begin(), k(kk).end(), rij_vec.begin(), 0.0); - std::complex tmp = { exchange_functional(rij)* (1.0-cos(kTwoPi*kr)), exchange_functional(rij) * sin(kTwoPi*kr)}; - spectrum_crystal_limit[kk] += tmp; if(kr != 0.0){ + std::complex tmp = { exchange_functional(rij)* (1.0-cos(kTwoPi*kr)), exchange_functional(rij) * sin(kTwoPi*kr)}; + spectrum_crystal_limit[kk] += tmp; counter2++; } } From 127cb8df300e7c2902cf750ad838ce12d046f739 Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Mon, 8 Jun 2020 14:43:36 +0900 Subject: [PATCH 28/80] Fixed radial dist functions. --- .../monitors/neutron_scattering_no_lattice.cc | 18 ++++++++++-------- .../monitors/neutron_scattering_no_lattice.h | 2 +- 2 files changed, 11 insertions(+), 9 deletions(-) diff --git a/src/jams/monitors/neutron_scattering_no_lattice.cc b/src/jams/monitors/neutron_scattering_no_lattice.cc index 2c963982..50147beb 100644 --- a/src/jams/monitors/neutron_scattering_no_lattice.cc +++ b/src/jams/monitors/neutron_scattering_no_lattice.cc @@ -420,31 +420,33 @@ NeutronScatteringNoLatticeMonitor::~NeutronScatteringNoLatticeMonitor() { // header for the radial dist file outfile4 << std::setw(20) << "bin" << "\t"; + outfile4 << std::setw(20) << "r/L" << "\t"; outfile4 << std::setw(20) << "dist" << "\n"; // header for the radial dist (z) file outfile5 << std::setw(20) << "bin" << "\t"; + outfile5 << std::setw(20) << "r/L" << "\t"; outfile5 << std::setw(20) << "dist" << "\n"; Vec3i lattice_dimensions = ::lattice->get_lattice_dimensions(); int num_distance = std::ceil(0.5 * lattice_dimensions[0] / delta_r_); //size of the array: (L/2)/delta_r_ - vector> r_dependence_SxSy(num_distance, 0.0); - vector> r_dependence_SzSz(num_distance, 0.0); - vector count_r_dep(r_dependence_SxSy.size(),0); - auto radial_dist = radial_distribution_function(); auto radial_dist_z = radial_distribution_function_z(); for (auto m = 0; m < radial_dist.size(); m++) { outfile4 << std::setw(20) << m << "\t"; + //Here we assume lattice_dimensions are all the same (x,y,z) + outfile4 << std::setw(20) << m * delta_r_ << "\t"; outfile4 << std::setw(20) << radial_dist[m] << "\n"; } - for (auto m = 0; m < radial_dist.size(); m++) { + for (auto m = 0; m < radial_dist_z.size(); m++) { outfile5 << std::setw(20) << m << "\t"; - outfile5 << std::setw(20) << radial_dist[m] << "\n"; + //Here we assume lattice_dimensions are all the same (x,y,z) + outfile5 << std::setw(20) << m * delta_r_ << "\t"; + outfile5 << std::setw(20) << radial_dist_z[m] << "\n"; } outfile4.close(); @@ -480,7 +482,7 @@ std::vector NeutronScatteringNoLatticeMonitor::radial_distribution_funct } for (auto kk = 0; kk < dist.size(); kk++){ dist[kk] /= num_spins; - cout << "dist[" << kk << "] = " << dist[kk] << ", dist/num_spins = " << dist[kk]/num_spins << endl; +// cout << "dist[" << kk << "] = " << dist[kk] << ", dist/num_spins = " << dist[kk]/num_spins << endl; } return dist; } @@ -509,7 +511,7 @@ std::vector NeutronScatteringNoLatticeMonitor::radial_distribution_funct } for (auto kk = 0; kk < dist_z.size(); kk++){ dist_z[kk] /= num_spins; - cout << "dist_z[" << kk << "] = " << dist_z[kk] << ", dist_z/num_spins = " << dist_z[kk]/num_spins << endl; +// cout << "dist_z[" << kk << "] = " << dist_z[kk] << ", dist_z/num_spins = " << dist_z[kk]/num_spins << endl; } return dist_z; } \ No newline at end of file diff --git a/src/jams/monitors/neutron_scattering_no_lattice.h b/src/jams/monitors/neutron_scattering_no_lattice.h index be8d1615..30acb269 100644 --- a/src/jams/monitors/neutron_scattering_no_lattice.h +++ b/src/jams/monitors/neutron_scattering_no_lattice.h @@ -55,7 +55,7 @@ class NeutronScatteringNoLatticeMonitor : public Monitor { int num_k_; Vec3 kvector_; bool evaluate_form_factor_ = false; - double delta_r_ = 0.01; + double delta_r_ = 0.001; int t_diff_r_dep_; std::vector> spin_x; From b1a9623f46a5c9114de1f725980b38c6d1d8c04d Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Mon, 8 Jun 2020 15:43:39 +0900 Subject: [PATCH 29/80] Fixed radial dist functions 2 --- src/jams/monitors/neutron_scattering_no_lattice.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/jams/monitors/neutron_scattering_no_lattice.cc b/src/jams/monitors/neutron_scattering_no_lattice.cc index 50147beb..da1087e5 100644 --- a/src/jams/monitors/neutron_scattering_no_lattice.cc +++ b/src/jams/monitors/neutron_scattering_no_lattice.cc @@ -501,7 +501,7 @@ std::vector NeutronScatteringNoLatticeMonitor::radial_distribution_funct continue; } const auto r_i = lattice->atom_position(i); - const auto r_ij = (lattice->displacement(r_i, r_j)); // using j,i in this order gives r_j - r_i + const auto r_ij = (lattice->displacement(r_j, r_i)); int m = std::ceil(r_ij[2]/delta_r_); if(m > dist_z.size()){ continue; From 5b83f9c5f8920ce8434ccbc962bd990bea66ccca Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Mon, 8 Jun 2020 15:46:05 +0900 Subject: [PATCH 30/80] Fixed radial dist functions 3 --- src/jams/monitors/neutron_scattering_no_lattice.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/jams/monitors/neutron_scattering_no_lattice.cc b/src/jams/monitors/neutron_scattering_no_lattice.cc index da1087e5..08624730 100644 --- a/src/jams/monitors/neutron_scattering_no_lattice.cc +++ b/src/jams/monitors/neutron_scattering_no_lattice.cc @@ -502,7 +502,7 @@ std::vector NeutronScatteringNoLatticeMonitor::radial_distribution_funct } const auto r_i = lattice->atom_position(i); const auto r_ij = (lattice->displacement(r_j, r_i)); - int m = std::ceil(r_ij[2]/delta_r_); + int m = std::ceil(abs(r_ij[2])/delta_r_); if(m > dist_z.size()){ continue; } From 02c51b8551b0d094b12458f1970cee9cbc5b8397 Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Mon, 8 Jun 2020 16:19:46 +0900 Subject: [PATCH 31/80] Fixed radial dist functions 4 --- .../monitors/neutron_scattering_no_lattice.cc | 77 +++++++------------ 1 file changed, 29 insertions(+), 48 deletions(-) diff --git a/src/jams/monitors/neutron_scattering_no_lattice.cc b/src/jams/monitors/neutron_scattering_no_lattice.cc index 08624730..d4ae988f 100644 --- a/src/jams/monitors/neutron_scattering_no_lattice.cc +++ b/src/jams/monitors/neutron_scattering_no_lattice.cc @@ -414,26 +414,16 @@ NeutronScatteringNoLatticeMonitor::~NeutronScatteringNoLatticeMonitor() { ofstream outfile4(name4.c_str()); outfile4.setf(std::ios::right); - string name5 = seedname + "_radial_dist_z.tsv"; - ofstream outfile5(name5.c_str()); - outfile5.setf(std::ios::right); - // header for the radial dist file outfile4 << std::setw(20) << "bin" << "\t"; outfile4 << std::setw(20) << "r/L" << "\t"; outfile4 << std::setw(20) << "dist" << "\n"; - // header for the radial dist (z) file - outfile5 << std::setw(20) << "bin" << "\t"; - outfile5 << std::setw(20) << "r/L" << "\t"; - outfile5 << std::setw(20) << "dist" << "\n"; - Vec3i lattice_dimensions = ::lattice->get_lattice_dimensions(); int num_distance = std::ceil(0.5 * lattice_dimensions[0] / delta_r_); //size of the array: (L/2)/delta_r_ auto radial_dist = radial_distribution_function(); - auto radial_dist_z = radial_distribution_function_z(); for (auto m = 0; m < radial_dist.size(); m++) { outfile4 << std::setw(20) << m << "\t"; @@ -442,15 +432,7 @@ NeutronScatteringNoLatticeMonitor::~NeutronScatteringNoLatticeMonitor() { outfile4 << std::setw(20) << radial_dist[m] << "\n"; } - for (auto m = 0; m < radial_dist_z.size(); m++) { - outfile5 << std::setw(20) << m << "\t"; - //Here we assume lattice_dimensions are all the same (x,y,z) - outfile5 << std::setw(20) << m * delta_r_ << "\t"; - outfile5 << std::setw(20) << radial_dist_z[m] << "\n"; - } - outfile4.close(); - outfile5.close(); } } @@ -464,6 +446,11 @@ std::vector NeutronScatteringNoLatticeMonitor::radial_distribution_funct Vec3i lattice_dimensions = ::lattice->get_lattice_dimensions(); //Here we assume lattice_dimensions are all the same (x,y,z) std::vector dist(std::ceil(0.5 * lattice_dimensions[0] / delta_r_),0.0); + jams::MultiArray count_z; + count_z.resize(dist.size()); + for(auto mm=0; mm < count_z.size(); mm++){ + count_z(mm) = 0.0; + } for(unsigned j = 0; j < num_spins; ++j){ const auto r_j = lattice->atom_position(j); for(unsigned i = 0; i < num_spins; ++i){ @@ -474,44 +461,38 @@ std::vector NeutronScatteringNoLatticeMonitor::radial_distribution_funct const auto r_i = lattice->atom_position(i); const auto r_ij = (lattice->displacement(r_i, r_j)); // using j,i in this order gives r_j - r_i int m = std::ceil(std::sqrt(r_ij[0] * r_ij[0] + r_ij[1] * r_ij[1] + r_ij[2] * r_ij[2])/delta_r_); - if(m > dist.size()){ + if(m >= dist.size()){ continue; } dist[m]++; + if(abs(r_ij[2]) > 0.0){ + count_z(m)++; + } } } for (auto kk = 0; kk < dist.size(); kk++){ dist[kk] /= num_spins; -// cout << "dist[" << kk << "] = " << dist[kk] << ", dist/num_spins = " << dist[kk]/num_spins << endl; + count_z(kk) /= num_spins; +// cout << "dist[" << kk << "]/num_spins = " << dist[kk] << endl; +// cout << "dist_z[" << kk << "]/num_spins = " << count_z(kk) << endl; } - return dist; -} -std::vector NeutronScatteringNoLatticeMonitor::radial_distribution_function_z() { - using namespace globals; - using namespace std; - Vec3i lattice_dimensions = ::lattice->get_lattice_dimensions(); - //Here we assume lattice_dimensions are all the same (x,y,z) - std::vector dist_z(std::ceil(0.5 * lattice_dimensions[0] / delta_r_),0.0); - for(unsigned j = 0; j < num_spins; ++j){ - const auto r_j = lattice->atom_position(j); - for(unsigned i = 0; i < num_spins; ++i){ - if(i == j){ - dist_z[0] = 0; - continue; - } - const auto r_i = lattice->atom_position(i); - const auto r_ij = (lattice->displacement(r_j, r_i)); - int m = std::ceil(abs(r_ij[2])/delta_r_); - if(m > dist_z.size()){ - continue; - } - dist_z[m]++; - } - } - for (auto kk = 0; kk < dist_z.size(); kk++){ - dist_z[kk] /= num_spins; -// cout << "dist_z[" << kk << "] = " << dist_z[kk] << ", dist_z/num_spins = " << dist_z[kk]/num_spins << endl; + string name5 = seedname + "_radial_dist_z.tsv"; + ofstream outfile5(name5.c_str()); + outfile5.setf(std::ios::right); + + // header for the radial dist (z) file + outfile5 << std::setw(20) << "bin" << "\t"; + outfile5 << std::setw(20) << "r/L" << "\t"; + outfile5 << std::setw(20) << "z/L" << "\t"; + outfile5 << std::setw(20) << "dist" << "\n"; + + for (auto m = 0; m < dist.size(); m++) { + outfile5 << std::setw(20) << m << "\t"; + //Here we assume lattice_dimensions are all the same (x,y,z) + outfile5 << std::setw(20) << m * delta_r_ << "\t"; + outfile5 << std::setw(20) << count_z(m) << "\n"; } - return dist_z; + outfile5.close(); + return dist; } \ No newline at end of file From 982c42382775fa795d2e444648ea54431ad4017d Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Mon, 8 Jun 2020 16:25:52 +0900 Subject: [PATCH 32/80] Fixed radial dist functions 5 --- src/jams/monitors/neutron_scattering_no_lattice.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/jams/monitors/neutron_scattering_no_lattice.cc b/src/jams/monitors/neutron_scattering_no_lattice.cc index d4ae988f..dd80463e 100644 --- a/src/jams/monitors/neutron_scattering_no_lattice.cc +++ b/src/jams/monitors/neutron_scattering_no_lattice.cc @@ -465,7 +465,7 @@ std::vector NeutronScatteringNoLatticeMonitor::radial_distribution_funct continue; } dist[m]++; - if(abs(r_ij[2]) > 0.0){ + if(abs(r_ij[2]) > (delta_r_ * 10.0)){ count_z(m)++; } } From d41f99a33084ba53b53cdda23a6599cbcf1e087f Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Mon, 8 Jun 2020 23:00:43 +0900 Subject: [PATCH 33/80] Added function calculating crystal limit spectrum ver2. --- src/jams/hamiltonian/exchange_functional.cc | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/jams/hamiltonian/exchange_functional.cc b/src/jams/hamiltonian/exchange_functional.cc index bf397572..dbc3a156 100644 --- a/src/jams/hamiltonian/exchange_functional.cc +++ b/src/jams/hamiltonian/exchange_functional.cc @@ -39,6 +39,7 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se // --- for crystal limit spectrum --- int num_k = jams::config_required(settings, "num_k"); std::vector> spectrum_crystal_limit(num_k+1,{0.0,0.0}); + std::vector> spectrum_crystal_limit2(num_k+1,{0.0,0.0}); double kmax = jams::config_required(settings, "kmax"); Vec3 kvector = jams::config_required(settings, "kvector"); jams::MultiArray k; @@ -66,7 +67,9 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se double kr = std::inner_product(k(kk).begin(), k(kk).end(), rij_vec.begin(), 0.0); if(kr != 0.0){ std::complex tmp = { exchange_functional(rij)* (1.0-cos(kTwoPi*kr)), exchange_functional(rij) * sin(kTwoPi*kr)}; + std::complex tmp2 = { exchange_functional(rij)* (1.0-cos(kTwoPi*kr)) /(4*kPi*rij*rij), exchange_functional(rij) * sin(kTwoPi*kr)/(4*kPi*rij*rij)}; spectrum_crystal_limit[kk] += tmp; + spectrum_crystal_limit2[kk] += tmp2; counter2++; } } @@ -75,6 +78,7 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se } for (auto kk = 0; kk < spectrum_crystal_limit.size(); kk++) { spectrum_crystal_limit[kk] /= globals::num_spins; + spectrum_crystal_limit2[kk] /= globals::num_spins; } cout << " total interactions " << jams::fmt::integer << counter << "\n"; @@ -90,6 +94,8 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se outfile3 << std::setw(20) << k(m)[2] << "\t"; outfile3 << std::setw(20) << spectrum_crystal_limit[m].real() << "\t"; outfile3 << std::setw(20) << spectrum_crystal_limit[m].imag() << "\n"; + outfile3 << std::setw(20) << spectrum_crystal_limit2[m].real() << "\t"; + outfile3 << std::setw(20) << spectrum_crystal_limit2[m].imag() << "\n"; } // --- for crystal limit spectrum --- From a2a10cc8bf93832831bd76583867a0edd85111cd Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Tue, 9 Jun 2020 16:59:42 +0900 Subject: [PATCH 34/80] [Fixed] Added function calculating crystal limit spectrum ver2. --- src/jams/hamiltonian/exchange_functional.cc | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/jams/hamiltonian/exchange_functional.cc b/src/jams/hamiltonian/exchange_functional.cc index dbc3a156..0ad82124 100644 --- a/src/jams/hamiltonian/exchange_functional.cc +++ b/src/jams/hamiltonian/exchange_functional.cc @@ -31,7 +31,9 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se outfile3 << std::setw(20) << "ky" << "\t";//(\t=tab) outfile3 << std::setw(20) << "kz" << "\t";//(\t=tab) outfile3 << std::setw(20) << "Re:E(k)" << "\t"; - outfile3 << std::setw(20) << "Im:E(k)" << "\n"; + outfile3 << std::setw(20) << "Im:E(k)" << "\t"; + outfile3 << std::setw(20) << "Re:E(k)-2" << "\t"; + outfile3 << std::setw(20) << "Im:E(k)-2" << "\n"; auto counter = 0; auto counter2 = 0; @@ -93,7 +95,7 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se outfile3 << std::setw(20) << k(m)[1] << "\t"; outfile3 << std::setw(20) << k(m)[2] << "\t"; outfile3 << std::setw(20) << spectrum_crystal_limit[m].real() << "\t"; - outfile3 << std::setw(20) << spectrum_crystal_limit[m].imag() << "\n"; + outfile3 << std::setw(20) << spectrum_crystal_limit[m].imag() << "\t"; outfile3 << std::setw(20) << spectrum_crystal_limit2[m].real() << "\t"; outfile3 << std::setw(20) << spectrum_crystal_limit2[m].imag() << "\n"; } From a407e663507206d5539484d90c70adb2a53587c2 Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Thu, 11 Jun 2020 16:41:44 +0900 Subject: [PATCH 35/80] Added 3rd version of crystal limit. --- src/jams/hamiltonian/exchange_functional.cc | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/jams/hamiltonian/exchange_functional.cc b/src/jams/hamiltonian/exchange_functional.cc index 0ad82124..b04f2333 100644 --- a/src/jams/hamiltonian/exchange_functional.cc +++ b/src/jams/hamiltonian/exchange_functional.cc @@ -34,6 +34,8 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se outfile3 << std::setw(20) << "Im:E(k)" << "\t"; outfile3 << std::setw(20) << "Re:E(k)-2" << "\t"; outfile3 << std::setw(20) << "Im:E(k)-2" << "\n"; + outfile3 << std::setw(20) << "Re:E(k)-3" << "\n"; + outfile3 << std::setw(20) << "Im:E(k)-3" << "\n"; auto counter = 0; auto counter2 = 0; @@ -42,6 +44,7 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se int num_k = jams::config_required(settings, "num_k"); std::vector> spectrum_crystal_limit(num_k+1,{0.0,0.0}); std::vector> spectrum_crystal_limit2(num_k+1,{0.0,0.0}); + std::vector> spectrum_crystal_limit_noave(num_k+1,{0.0,0.0}); double kmax = jams::config_required(settings, "kmax"); Vec3 kvector = jams::config_required(settings, "kvector"); jams::MultiArray k; @@ -72,6 +75,9 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se std::complex tmp2 = { exchange_functional(rij)* (1.0-cos(kTwoPi*kr)) /(4*kPi*rij*rij), exchange_functional(rij) * sin(kTwoPi*kr)/(4*kPi*rij*rij)}; spectrum_crystal_limit[kk] += tmp; spectrum_crystal_limit2[kk] += tmp2; + if(i == 0){ + spectrum_crystal_limit_noave[kk] += tmp; + } counter2++; } } @@ -98,6 +104,8 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se outfile3 << std::setw(20) << spectrum_crystal_limit[m].imag() << "\t"; outfile3 << std::setw(20) << spectrum_crystal_limit2[m].real() << "\t"; outfile3 << std::setw(20) << spectrum_crystal_limit2[m].imag() << "\n"; + outfile3 << std::setw(20) << spectrum_crystal_limit_noave[m].real() << "\n"; + outfile3 << std::setw(20) << spectrum_crystal_limit_noave[m].imag() << "\n"; } // --- for crystal limit spectrum --- From 3360f854e0f3c3b35eadf143eed5550aa2d51fa4 Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Thu, 11 Jun 2020 18:24:18 +0900 Subject: [PATCH 36/80] We can change r_ij range in the crystal limit. --- src/jams/hamiltonian/exchange_functional.cc | 31 +++++++++++++-------- 1 file changed, 19 insertions(+), 12 deletions(-) diff --git a/src/jams/hamiltonian/exchange_functional.cc b/src/jams/hamiltonian/exchange_functional.cc index b04f2333..d2fe269e 100644 --- a/src/jams/hamiltonian/exchange_functional.cc +++ b/src/jams/hamiltonian/exchange_functional.cc @@ -47,6 +47,11 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se std::vector> spectrum_crystal_limit_noave(num_k+1,{0.0,0.0}); double kmax = jams::config_required(settings, "kmax"); Vec3 kvector = jams::config_required(settings, "kvector"); + double rij_min = jams::config_optional(settings, "rij_min", 0.0); + double rij_max = jams::config_optional(settings, "rij_max", radius_cutoff_); + cout << "crystal limit: rij_min = " << rij_min << endl; + cout << "crystal limit: rij_max = " << rij_max << endl; + jams::MultiArray k; k.resize(num_k+1); for (auto n = 0; n < k.size(); ++n) { @@ -64,21 +69,23 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se continue; } const auto rij = norm(lattice->displacement(i, j)); - // --- for crystal limit spectrum --- - const auto rij_vec = lattice->displacement(i, j); this->insert_interaction_scalar(i, j, input_unit_conversion_ * exchange_functional(rij)); counter++; - for (auto kk = 0; kk < spectrum_crystal_limit.size(); kk++){ - double kr = std::inner_product(k(kk).begin(), k(kk).end(), rij_vec.begin(), 0.0); - if(kr != 0.0){ - std::complex tmp = { exchange_functional(rij)* (1.0-cos(kTwoPi*kr)), exchange_functional(rij) * sin(kTwoPi*kr)}; - std::complex tmp2 = { exchange_functional(rij)* (1.0-cos(kTwoPi*kr)) /(4*kPi*rij*rij), exchange_functional(rij) * sin(kTwoPi*kr)/(4*kPi*rij*rij)}; - spectrum_crystal_limit[kk] += tmp; - spectrum_crystal_limit2[kk] += tmp2; - if(i == 0){ - spectrum_crystal_limit_noave[kk] += tmp; + // --- for crystal limit spectrum --- + if(rij < rij_max && rij > rij_min){ + const auto rij_vec = lattice->displacement(i, j); + for (auto kk = 0; kk < spectrum_crystal_limit.size(); kk++){ + double kr = std::inner_product(k(kk).begin(), k(kk).end(), rij_vec.begin(), 0.0); + if(kr != 0.0){ + std::complex tmp = { exchange_functional(rij)* (1.0-cos(kTwoPi*kr)), exchange_functional(rij) * sin(kTwoPi*kr)}; + std::complex tmp2 = { exchange_functional(rij)* (1.0-cos(kTwoPi*kr)) /(4*kPi*rij*rij), exchange_functional(rij) * sin(kTwoPi*kr)/(4*kPi*rij*rij)}; + spectrum_crystal_limit[kk] += tmp; + spectrum_crystal_limit2[kk] += tmp2; + if(i == 0){ + spectrum_crystal_limit_noave[kk] += tmp; + } + counter2++; } - counter2++; } } // --- for crystal limit spectrum --- From 2a30a445edd32d2a90a990f24ed3f5938a9e5671 Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Thu, 11 Jun 2020 18:44:01 +0900 Subject: [PATCH 37/80] [Fixed] We can change r_ij range in the crystal limit. --- src/jams/hamiltonian/exchange_functional.cc | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/jams/hamiltonian/exchange_functional.cc b/src/jams/hamiltonian/exchange_functional.cc index d2fe269e..4d08d034 100644 --- a/src/jams/hamiltonian/exchange_functional.cc +++ b/src/jams/hamiltonian/exchange_functional.cc @@ -33,8 +33,8 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se outfile3 << std::setw(20) << "Re:E(k)" << "\t"; outfile3 << std::setw(20) << "Im:E(k)" << "\t"; outfile3 << std::setw(20) << "Re:E(k)-2" << "\t"; - outfile3 << std::setw(20) << "Im:E(k)-2" << "\n"; - outfile3 << std::setw(20) << "Re:E(k)-3" << "\n"; + outfile3 << std::setw(20) << "Im:E(k)-2" << "\t"; + outfile3 << std::setw(20) << "Re:E(k)-3" << "\t"; outfile3 << std::setw(20) << "Im:E(k)-3" << "\n"; auto counter = 0; @@ -110,8 +110,8 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se outfile3 << std::setw(20) << spectrum_crystal_limit[m].real() << "\t"; outfile3 << std::setw(20) << spectrum_crystal_limit[m].imag() << "\t"; outfile3 << std::setw(20) << spectrum_crystal_limit2[m].real() << "\t"; - outfile3 << std::setw(20) << spectrum_crystal_limit2[m].imag() << "\n"; - outfile3 << std::setw(20) << spectrum_crystal_limit_noave[m].real() << "\n"; + outfile3 << std::setw(20) << spectrum_crystal_limit2[m].imag() << "\t"; + outfile3 << std::setw(20) << spectrum_crystal_limit_noave[m].real() << "\t"; outfile3 << std::setw(20) << spectrum_crystal_limit_noave[m].imag() << "\n"; } // --- for crystal limit spectrum --- From 1a28fa596ae50f491c1a8835c68cf928ea2663d6 Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Thu, 11 Jun 2020 19:12:50 +0900 Subject: [PATCH 38/80] [Fixed-2] We can change r_ij range in the crystal limit. --- src/jams/hamiltonian/exchange_functional.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/jams/hamiltonian/exchange_functional.cc b/src/jams/hamiltonian/exchange_functional.cc index 4d08d034..079aa188 100644 --- a/src/jams/hamiltonian/exchange_functional.cc +++ b/src/jams/hamiltonian/exchange_functional.cc @@ -81,7 +81,7 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se std::complex tmp2 = { exchange_functional(rij)* (1.0-cos(kTwoPi*kr)) /(4*kPi*rij*rij), exchange_functional(rij) * sin(kTwoPi*kr)/(4*kPi*rij*rij)}; spectrum_crystal_limit[kk] += tmp; spectrum_crystal_limit2[kk] += tmp2; - if(i == 0){ + if(i == 25000){ spectrum_crystal_limit_noave[kk] += tmp; } counter2++; From 472125aecf220d7b86609c77a71c2a66a8e04232 Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Thu, 11 Jun 2020 23:25:31 +0900 Subject: [PATCH 39/80] Central site can be specified via cfg file. --- src/jams/hamiltonian/exchange_functional.cc | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/jams/hamiltonian/exchange_functional.cc b/src/jams/hamiltonian/exchange_functional.cc index 079aa188..d607c931 100644 --- a/src/jams/hamiltonian/exchange_functional.cc +++ b/src/jams/hamiltonian/exchange_functional.cc @@ -49,8 +49,10 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se Vec3 kvector = jams::config_required(settings, "kvector"); double rij_min = jams::config_optional(settings, "rij_min", 0.0); double rij_max = jams::config_optional(settings, "rij_max", radius_cutoff_); + int central_site = jams::config_optional(settings, "central_site", 0); cout << "crystal limit: rij_min = " << rij_min << endl; cout << "crystal limit: rij_max = " << rij_max << endl; + cout << "central_site: i = " << central_site << endl; jams::MultiArray k; k.resize(num_k+1); @@ -81,7 +83,7 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se std::complex tmp2 = { exchange_functional(rij)* (1.0-cos(kTwoPi*kr)) /(4*kPi*rij*rij), exchange_functional(rij) * sin(kTwoPi*kr)/(4*kPi*rij*rij)}; spectrum_crystal_limit[kk] += tmp; spectrum_crystal_limit2[kk] += tmp2; - if(i == 25000){ + if(i == central_site){ spectrum_crystal_limit_noave[kk] += tmp; } counter2++; From 6eb3185a51d6a1f07e5f8edff1f5eb0fbae7945b Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Mon, 15 Jun 2020 18:08:14 +0900 Subject: [PATCH 40/80] Added crystal limit calculation where central site = center of the system. --- src/jams/hamiltonian/exchange_functional.cc | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/src/jams/hamiltonian/exchange_functional.cc b/src/jams/hamiltonian/exchange_functional.cc index d607c931..76d16394 100644 --- a/src/jams/hamiltonian/exchange_functional.cc +++ b/src/jams/hamiltonian/exchange_functional.cc @@ -45,6 +45,7 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se std::vector> spectrum_crystal_limit(num_k+1,{0.0,0.0}); std::vector> spectrum_crystal_limit2(num_k+1,{0.0,0.0}); std::vector> spectrum_crystal_limit_noave(num_k+1,{0.0,0.0}); + std::vector> spectrum_crystal_limit_noave_central_cite(num_k+1,{0.0,0.0}); double kmax = jams::config_required(settings, "kmax"); Vec3 kvector = jams::config_required(settings, "kvector"); double rij_min = jams::config_optional(settings, "rij_min", 0.0); @@ -65,6 +66,12 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se nbrs.clear(); ::lattice->atom_neighbours(i, radius_cutoff_, nbrs); + rspace_displacement_(i) = lattice->displacement({lattice_dimensions[0]*0.5,lattice_dimensions[1]*0.5,lattice_dimensions[2]*0.5}, lattice->atom_position(i)); + if( norm(rspace_displacement_(i)) < 0.001023 ){ + int central_cite_indx = i; + cout << "central_site index: i = " << central_site_indx << endl; + } + for (const auto& nbr : nbrs) { const auto j = nbr.id; if (i == j) { @@ -86,6 +93,9 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se if(i == central_site){ spectrum_crystal_limit_noave[kk] += tmp; } + else if(i == central_site_indx){ + spectrum_crystal_limit_noave_central_cite[kk] += tmp; + } counter2++; } } @@ -114,7 +124,9 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se outfile3 << std::setw(20) << spectrum_crystal_limit2[m].real() << "\t"; outfile3 << std::setw(20) << spectrum_crystal_limit2[m].imag() << "\t"; outfile3 << std::setw(20) << spectrum_crystal_limit_noave[m].real() << "\t"; - outfile3 << std::setw(20) << spectrum_crystal_limit_noave[m].imag() << "\n"; + outfile3 << std::setw(20) << spectrum_crystal_limit_noave[m].imag() << "\t"; + outfile3 << std::setw(20) << spectrum_crystal_limit_noave_central_cite[m].real() << "\t"; + outfile3 << std::setw(20) << spectrum_crystal_limit_noave_central_cite[m].imag() << "\n"; } // --- for crystal limit spectrum --- From d5bb71febab1fa04efebf52c3151e36e6d46b178 Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Mon, 15 Jun 2020 18:12:01 +0900 Subject: [PATCH 41/80] [Fixed] Added crystal limit calculation where central site = center of the system. --- src/jams/hamiltonian/exchange_functional.cc | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/jams/hamiltonian/exchange_functional.cc b/src/jams/hamiltonian/exchange_functional.cc index 76d16394..997338ce 100644 --- a/src/jams/hamiltonian/exchange_functional.cc +++ b/src/jams/hamiltonian/exchange_functional.cc @@ -54,6 +54,7 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se cout << "crystal limit: rij_min = " << rij_min << endl; cout << "crystal limit: rij_max = " << rij_max << endl; cout << "central_site: i = " << central_site << endl; + int central_site_indx; jams::MultiArray k; k.resize(num_k+1); @@ -68,7 +69,7 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se rspace_displacement_(i) = lattice->displacement({lattice_dimensions[0]*0.5,lattice_dimensions[1]*0.5,lattice_dimensions[2]*0.5}, lattice->atom_position(i)); if( norm(rspace_displacement_(i)) < 0.001023 ){ - int central_cite_indx = i; + central_site_indx = i; cout << "central_site index: i = " << central_site_indx << endl; } From 7fa6b7a16cf70a1b8ce92d403354105ac7394570 Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Mon, 15 Jun 2020 18:14:35 +0900 Subject: [PATCH 42/80] [Fixed-2] Added crystal limit calculation where central site = center of the system. --- src/jams/hamiltonian/exchange_functional.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/jams/hamiltonian/exchange_functional.cc b/src/jams/hamiltonian/exchange_functional.cc index 997338ce..0835a9bd 100644 --- a/src/jams/hamiltonian/exchange_functional.cc +++ b/src/jams/hamiltonian/exchange_functional.cc @@ -83,7 +83,7 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se counter++; // --- for crystal limit spectrum --- if(rij < rij_max && rij > rij_min){ - const auto rij_vec = lattice->displacement(i, j); + const auto rij_vec = ::lattice->displacement(i, j); for (auto kk = 0; kk < spectrum_crystal_limit.size(); kk++){ double kr = std::inner_product(k(kk).begin(), k(kk).end(), rij_vec.begin(), 0.0); if(kr != 0.0){ From 98b3b68762886ae36e12118b45ea3947fa153bd4 Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Mon, 15 Jun 2020 18:28:21 +0900 Subject: [PATCH 43/80] [Fixed-3] Added crystal limit calculation where central site = center of the system. --- src/jams/hamiltonian/exchange_functional.cc | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/jams/hamiltonian/exchange_functional.cc b/src/jams/hamiltonian/exchange_functional.cc index 0835a9bd..3005e6c7 100644 --- a/src/jams/hamiltonian/exchange_functional.cc +++ b/src/jams/hamiltonian/exchange_functional.cc @@ -54,7 +54,9 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se cout << "crystal limit: rij_min = " << rij_min << endl; cout << "crystal limit: rij_max = " << rij_max << endl; cout << "central_site: i = " << central_site << endl; - int central_site_indx; + int central_site_indx = 0; + jams::MultiArray rspace_displacement_; + rspace_displacement_.resize(globals::s.size(0)); jams::MultiArray k; k.resize(num_k+1); @@ -67,7 +69,7 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se nbrs.clear(); ::lattice->atom_neighbours(i, radius_cutoff_, nbrs); - rspace_displacement_(i) = lattice->displacement({lattice_dimensions[0]*0.5,lattice_dimensions[1]*0.5,lattice_dimensions[2]*0.5}, lattice->atom_position(i)); + rspace_displacement_(i) = ::lattice->displacement({0.5,0.5,0.5}, lattice->atom_position(i)); if( norm(rspace_displacement_(i)) < 0.001023 ){ central_site_indx = i; cout << "central_site index: i = " << central_site_indx << endl; @@ -78,7 +80,7 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se if (i == j) { continue; } - const auto rij = norm(lattice->displacement(i, j)); + const auto rij = norm(::lattice->displacement(i, j)); this->insert_interaction_scalar(i, j, input_unit_conversion_ * exchange_functional(rij)); counter++; // --- for crystal limit spectrum --- From 7b313ca09b077247546d2ba2d6605315fda0b7fc Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Mon, 15 Jun 2020 23:41:34 +0900 Subject: [PATCH 44/80] [Fixed-4] Added crystal limit calculation where central site = center of the system. --- src/jams/hamiltonian/exchange_functional.cc | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/src/jams/hamiltonian/exchange_functional.cc b/src/jams/hamiltonian/exchange_functional.cc index 3005e6c7..c14d1f1e 100644 --- a/src/jams/hamiltonian/exchange_functional.cc +++ b/src/jams/hamiltonian/exchange_functional.cc @@ -35,7 +35,9 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se outfile3 << std::setw(20) << "Re:E(k)-2" << "\t"; outfile3 << std::setw(20) << "Im:E(k)-2" << "\t"; outfile3 << std::setw(20) << "Re:E(k)-3" << "\t"; - outfile3 << std::setw(20) << "Im:E(k)-3" << "\n"; + outfile3 << std::setw(20) << "Im:E(k)-3" << "\t"; + outfile3 << std::setw(20) << "Re:E(k)-4" << "\t"; + outfile3 << std::setw(20) << "Im:E(k)-4" << "\n"; auto counter = 0; auto counter2 = 0; @@ -69,8 +71,9 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se nbrs.clear(); ::lattice->atom_neighbours(i, radius_cutoff_, nbrs); - rspace_displacement_(i) = ::lattice->displacement({0.5,0.5,0.5}, lattice->atom_position(i)); - if( norm(rspace_displacement_(i)) < 0.001023 ){ + Vec3i lattice_dimensions = ::lattice->get_lattice_dimensions(); + rspace_displacement_(i) = ::lattice->displacement({lattice_dimensions[0]*0.5,lattice_dimensions[1]*0.5,lattice_dimensions[2]*0.5}, lattice->atom_position(i)); + if( norm(rspace_displacement_(i)) < lattice_dimensions[0]*0.001023 ){ central_site_indx = i; cout << "central_site index: i = " << central_site_indx << endl; } @@ -98,6 +101,7 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se } else if(i == central_site_indx){ spectrum_crystal_limit_noave_central_cite[kk] += tmp; + cout << "NOW!" << endl; } counter2++; } From c2e057c792ad81e7d0a1dc77c6afa1cab12bde5a Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Tue, 16 Jun 2020 08:59:44 +0900 Subject: [PATCH 45/80] [Fixed-5] Added crystal limit calculation where central site = center of the system. --- src/jams/hamiltonian/exchange_functional.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/jams/hamiltonian/exchange_functional.cc b/src/jams/hamiltonian/exchange_functional.cc index c14d1f1e..bdf85d76 100644 --- a/src/jams/hamiltonian/exchange_functional.cc +++ b/src/jams/hamiltonian/exchange_functional.cc @@ -73,7 +73,7 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se Vec3i lattice_dimensions = ::lattice->get_lattice_dimensions(); rspace_displacement_(i) = ::lattice->displacement({lattice_dimensions[0]*0.5,lattice_dimensions[1]*0.5,lattice_dimensions[2]*0.5}, lattice->atom_position(i)); - if( norm(rspace_displacement_(i)) < lattice_dimensions[0]*0.001023 ){ + if( norm(rspace_displacement_(i)) < lattice_dimensions[0]*0.002047 ){ central_site_indx = i; cout << "central_site index: i = " << central_site_indx << endl; } From 75ea9cedd6b51c4faf42cec0537da7fac9f99466 Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Tue, 16 Jun 2020 10:39:46 +0900 Subject: [PATCH 46/80] [Fixed-6] Added crystal limit calculation where central site = center of the system. --- src/jams/hamiltonian/exchange_functional.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/jams/hamiltonian/exchange_functional.cc b/src/jams/hamiltonian/exchange_functional.cc index bdf85d76..75832f02 100644 --- a/src/jams/hamiltonian/exchange_functional.cc +++ b/src/jams/hamiltonian/exchange_functional.cc @@ -73,7 +73,7 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se Vec3i lattice_dimensions = ::lattice->get_lattice_dimensions(); rspace_displacement_(i) = ::lattice->displacement({lattice_dimensions[0]*0.5,lattice_dimensions[1]*0.5,lattice_dimensions[2]*0.5}, lattice->atom_position(i)); - if( norm(rspace_displacement_(i)) < lattice_dimensions[0]*0.002047 ){ + if( norm(rspace_displacement_(i)) < lattice_dimensions[0]*0.005117 ){ //0.002047 central_site_indx = i; cout << "central_site index: i = " << central_site_indx << endl; } From 647037ed431fa492108af9e8e76bc418fec776cc Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Tue, 16 Jun 2020 10:43:28 +0900 Subject: [PATCH 47/80] [Fixed-7] Added crystal limit calculation where central site = center of the system. --- src/jams/hamiltonian/exchange_functional.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/jams/hamiltonian/exchange_functional.cc b/src/jams/hamiltonian/exchange_functional.cc index 75832f02..98b31014 100644 --- a/src/jams/hamiltonian/exchange_functional.cc +++ b/src/jams/hamiltonian/exchange_functional.cc @@ -73,7 +73,7 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se Vec3i lattice_dimensions = ::lattice->get_lattice_dimensions(); rspace_displacement_(i) = ::lattice->displacement({lattice_dimensions[0]*0.5,lattice_dimensions[1]*0.5,lattice_dimensions[2]*0.5}, lattice->atom_position(i)); - if( norm(rspace_displacement_(i)) < lattice_dimensions[0]*0.005117 ){ //0.002047 + if( norm(rspace_displacement_(i)) < lattice_dimensions[0]*0.1 ){ //0.002047 central_site_indx = i; cout << "central_site index: i = " << central_site_indx << endl; } From 23547bc3c6dc25092c200ae4c2db06f7d5dd8776 Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Tue, 16 Jun 2020 10:47:18 +0900 Subject: [PATCH 48/80] [Fixed-8] Added crystal limit calculation where central site = center of the system. --- src/jams/hamiltonian/exchange_functional.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/jams/hamiltonian/exchange_functional.cc b/src/jams/hamiltonian/exchange_functional.cc index 98b31014..a443dabc 100644 --- a/src/jams/hamiltonian/exchange_functional.cc +++ b/src/jams/hamiltonian/exchange_functional.cc @@ -73,7 +73,7 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se Vec3i lattice_dimensions = ::lattice->get_lattice_dimensions(); rspace_displacement_(i) = ::lattice->displacement({lattice_dimensions[0]*0.5,lattice_dimensions[1]*0.5,lattice_dimensions[2]*0.5}, lattice->atom_position(i)); - if( norm(rspace_displacement_(i)) < lattice_dimensions[0]*0.1 ){ //0.002047 + if( norm(rspace_displacement_(i)) < lattice_dimensions[0]*0.01 ){ //0.002047 central_site_indx = i; cout << "central_site index: i = " << central_site_indx << endl; } From 14f7b17a239120363e19e04c44a0342bbc9c8800 Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Tue, 16 Jun 2020 10:55:30 +0900 Subject: [PATCH 49/80] [Fixed-9] Added crystal limit calculation where central site = center of the system. --- src/jams/hamiltonian/exchange_functional.cc | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/jams/hamiltonian/exchange_functional.cc b/src/jams/hamiltonian/exchange_functional.cc index a443dabc..9db70c92 100644 --- a/src/jams/hamiltonian/exchange_functional.cc +++ b/src/jams/hamiltonian/exchange_functional.cc @@ -59,6 +59,8 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se int central_site_indx = 0; jams::MultiArray rspace_displacement_; rspace_displacement_.resize(globals::s.size(0)); + jams::MultiArray rspace_displacement2_; + rspace_displacement2_.resize(globals::s.size(0)); jams::MultiArray k; k.resize(num_k+1); @@ -73,9 +75,11 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se Vec3i lattice_dimensions = ::lattice->get_lattice_dimensions(); rspace_displacement_(i) = ::lattice->displacement({lattice_dimensions[0]*0.5,lattice_dimensions[1]*0.5,lattice_dimensions[2]*0.5}, lattice->atom_position(i)); + rspace_displacement2_(i) = ::lattice->displacement({0,0,0}, lattice->atom_position(i)); if( norm(rspace_displacement_(i)) < lattice_dimensions[0]*0.01 ){ //0.002047 central_site_indx = i; cout << "central_site index: i = " << central_site_indx << endl; + cout << "central coordinate = ( " << rspace_displacement2_(central_site_indx) << " )" << endl; } for (const auto& nbr : nbrs) { From b3a339401bfd0a9e08c47b5eff56a68c3d774285 Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Tue, 16 Jun 2020 20:45:29 +0900 Subject: [PATCH 50/80] [Fixed-10] Added crystal limit calculation where central site = center of the system. --- src/jams/hamiltonian/exchange_functional.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/jams/hamiltonian/exchange_functional.cc b/src/jams/hamiltonian/exchange_functional.cc index 9db70c92..a2902bc1 100644 --- a/src/jams/hamiltonian/exchange_functional.cc +++ b/src/jams/hamiltonian/exchange_functional.cc @@ -76,7 +76,7 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se Vec3i lattice_dimensions = ::lattice->get_lattice_dimensions(); rspace_displacement_(i) = ::lattice->displacement({lattice_dimensions[0]*0.5,lattice_dimensions[1]*0.5,lattice_dimensions[2]*0.5}, lattice->atom_position(i)); rspace_displacement2_(i) = ::lattice->displacement({0,0,0}, lattice->atom_position(i)); - if( norm(rspace_displacement_(i)) < lattice_dimensions[0]*0.01 ){ //0.002047 + if( norm(rspace_displacement_(i)) < lattice_dimensions[0]*0.05 ){ //0.002047 central_site_indx = i; cout << "central_site index: i = " << central_site_indx << endl; cout << "central coordinate = ( " << rspace_displacement2_(central_site_indx) << " )" << endl; From 88add6811c455804d481aae06e98db2e0886a1e9 Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Tue, 16 Jun 2020 20:57:34 +0900 Subject: [PATCH 51/80] [Fixed-11] Added crystal limit calculation where central site = center of the system. --- src/jams/hamiltonian/exchange_functional.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/jams/hamiltonian/exchange_functional.cc b/src/jams/hamiltonian/exchange_functional.cc index a2902bc1..8560d822 100644 --- a/src/jams/hamiltonian/exchange_functional.cc +++ b/src/jams/hamiltonian/exchange_functional.cc @@ -76,7 +76,7 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se Vec3i lattice_dimensions = ::lattice->get_lattice_dimensions(); rspace_displacement_(i) = ::lattice->displacement({lattice_dimensions[0]*0.5,lattice_dimensions[1]*0.5,lattice_dimensions[2]*0.5}, lattice->atom_position(i)); rspace_displacement2_(i) = ::lattice->displacement({0,0,0}, lattice->atom_position(i)); - if( norm(rspace_displacement_(i)) < lattice_dimensions[0]*0.05 ){ //0.002047 + if( norm(rspace_displacement_(i)) < lattice_dimensions[0]*0.03 ){ //0.002047 central_site_indx = i; cout << "central_site index: i = " << central_site_indx << endl; cout << "central coordinate = ( " << rspace_displacement2_(central_site_indx) << " )" << endl; From efd249eaba0e4b6fe4cf6fe27e24c39c4f375976 Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Tue, 16 Jun 2020 21:00:28 +0900 Subject: [PATCH 52/80] [Fixed-12] Added crystal limit calculation where central site = center of the system. --- src/jams/hamiltonian/exchange_functional.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/jams/hamiltonian/exchange_functional.cc b/src/jams/hamiltonian/exchange_functional.cc index 8560d822..9a0951b2 100644 --- a/src/jams/hamiltonian/exchange_functional.cc +++ b/src/jams/hamiltonian/exchange_functional.cc @@ -76,7 +76,7 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se Vec3i lattice_dimensions = ::lattice->get_lattice_dimensions(); rspace_displacement_(i) = ::lattice->displacement({lattice_dimensions[0]*0.5,lattice_dimensions[1]*0.5,lattice_dimensions[2]*0.5}, lattice->atom_position(i)); rspace_displacement2_(i) = ::lattice->displacement({0,0,0}, lattice->atom_position(i)); - if( norm(rspace_displacement_(i)) < lattice_dimensions[0]*0.03 ){ //0.002047 + if( norm(rspace_displacement_(i)) < lattice_dimensions[0]*0.02 ){ //0.002047 central_site_indx = i; cout << "central_site index: i = " << central_site_indx << endl; cout << "central coordinate = ( " << rspace_displacement2_(central_site_indx) << " )" << endl; From 695befb3e20346d2bdb18917265a2d8ccf476ae8 Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Tue, 16 Jun 2020 21:13:49 +0900 Subject: [PATCH 53/80] [Fixed-13] Added crystal limit calculation where central site = center of the system. --- src/jams/hamiltonian/exchange_functional.cc | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/jams/hamiltonian/exchange_functional.cc b/src/jams/hamiltonian/exchange_functional.cc index 9a0951b2..d5313b69 100644 --- a/src/jams/hamiltonian/exchange_functional.cc +++ b/src/jams/hamiltonian/exchange_functional.cc @@ -75,11 +75,12 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se Vec3i lattice_dimensions = ::lattice->get_lattice_dimensions(); rspace_displacement_(i) = ::lattice->displacement({lattice_dimensions[0]*0.5,lattice_dimensions[1]*0.5,lattice_dimensions[2]*0.5}, lattice->atom_position(i)); - rspace_displacement2_(i) = ::lattice->displacement({0,0,0}, lattice->atom_position(i)); + rspace_displacement2_(i) = ::lattice->displacement(lattice->atom_position(i),{0,0,0}); if( norm(rspace_displacement_(i)) < lattice_dimensions[0]*0.02 ){ //0.002047 central_site_indx = i; cout << "central_site index: i = " << central_site_indx << endl; - cout << "central coordinate = ( " << rspace_displacement2_(central_site_indx) << " )" << endl; + cout << "central coordinate = ( " << rspace_displacement_(central_site_indx) << " )" << endl; + cout << "central coordinate 2 = ( " << rspace_displacement2_(central_site_indx) << " )" << endl; } for (const auto& nbr : nbrs) { From 4c67a2f6ceaef9405bedb3d020e62618f45da2ad Mon Sep 17 00:00:00 2001 From: Mai Kameda Date: Tue, 16 Jun 2020 21:20:06 +0900 Subject: [PATCH 54/80] [Fixed-14] Added crystal limit calculation where central site = center of the system. --- src/jams/hamiltonian/exchange_functional.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/jams/hamiltonian/exchange_functional.cc b/src/jams/hamiltonian/exchange_functional.cc index d5313b69..3bea8945 100644 --- a/src/jams/hamiltonian/exchange_functional.cc +++ b/src/jams/hamiltonian/exchange_functional.cc @@ -75,7 +75,7 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se Vec3i lattice_dimensions = ::lattice->get_lattice_dimensions(); rspace_displacement_(i) = ::lattice->displacement({lattice_dimensions[0]*0.5,lattice_dimensions[1]*0.5,lattice_dimensions[2]*0.5}, lattice->atom_position(i)); - rspace_displacement2_(i) = ::lattice->displacement(lattice->atom_position(i),{0,0,0}); + rspace_displacement2_(i) = ::lattice->atom_position(i); if( norm(rspace_displacement_(i)) < lattice_dimensions[0]*0.02 ){ //0.002047 central_site_indx = i; cout << "central_site index: i = " << central_site_indx << endl; From e6aa1d896b208113b79ac7245b0fa39b71df814a Mon Sep 17 00:00:00 2001 From: mkameda <39291377+mkameda@users.noreply.github.com> Date: Wed, 24 Nov 2021 08:55:45 +0900 Subject: [PATCH 55/80] Added random potential --- src/jams/hamiltonian/exchange_functional.cc | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/src/jams/hamiltonian/exchange_functional.cc b/src/jams/hamiltonian/exchange_functional.cc index e767affc..cbc45ef1 100644 --- a/src/jams/hamiltonian/exchange_functional.cc +++ b/src/jams/hamiltonian/exchange_functional.cc @@ -61,6 +61,16 @@ double ExchangeFunctionalHamiltonian::functional_kaneyoshi(const double rij, con return J0 * pow2(rij - r0) * exp(-pow2(rij - r0) / (2 * pow2(sigma))); } +double ExchangeFunctionalHamiltonian::functional_random(const double rij, const double J0, const double rc, const double width){ + std::uniform_real_distribution<> rand_potential(0.0, width * J0); //0.0 < width < 1.0 //rc is cutoff for the random potential + std::mt19937 rand_src(12345); //seed=12345 + if(rij < rc) { + return rand_potential(rand_src); + } else{ + return 0.0; + } +} + ExchangeFunctionalHamiltonian::ExchangeFunctional ExchangeFunctionalHamiltonian::functional_from_settings(const libconfig::Setting &settings) { using namespace std::placeholders; @@ -76,6 +86,8 @@ ExchangeFunctionalHamiltonian::functional_from_settings(const libconfig::Setting return bind(functional_gaussian, _1, double(settings["J0"]), double(settings["r0"]), double(settings["sigma"])); } else if (functional_name == "kaneyoshi") { return bind(functional_kaneyoshi, _1, double(settings["J0"]), double(settings["r0"]), double(settings["sigma"])); + } else if (functional_name == "random") { + return bind(functional_random, _1, double(settings["J0"]), double(settings["rc"]), double(settings["width"])); } else { throw runtime_error("unknown exchange functional: " + functional_name); } From dad06600677e62ad1b7e98616ce94735b02526b1 Mon Sep 17 00:00:00 2001 From: mkameda <39291377+mkameda@users.noreply.github.com> Date: Wed, 24 Nov 2021 10:29:53 +0900 Subject: [PATCH 56/80] [Fixed] Added random potential --- src/jams/hamiltonian/exchange_functional.cc | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/src/jams/hamiltonian/exchange_functional.cc b/src/jams/hamiltonian/exchange_functional.cc index 3bea8945..4dcb547a 100644 --- a/src/jams/hamiltonian/exchange_functional.cc +++ b/src/jams/hamiltonian/exchange_functional.cc @@ -174,6 +174,16 @@ double ExchangeFunctionalHamiltonian::functional_gaussian_multi(const double rij return J0 * exp(-pow2(rij - r0)/(2 * pow2(sigma))) + J0_2 * exp(-pow2(rij - r0_2)/(2 * pow2(sigma_2))) + J0_3 * exp(-pow2(rij - r0_3)/(2 * pow2(sigma_3))); } +double ExchangeFunctionalHamiltonian::functional_random(const double rij, const double J0, const double rc, const double width){ + std::uniform_real_distribution<> rand_potential(0.0, width * J0); //0.0 < width < 1.0 //rc is cutoff for the random potential + std::mt19937 rand_src(12345); //seed=12345 + if(rij < rc) { + return rand_potential(rand_src); + } else{ + return 0.0; + } +} + ExchangeFunctionalHamiltonian::ExchangeFunctional ExchangeFunctionalHamiltonian::functional_from_settings(const libconfig::Setting &settings) { using namespace std::placeholders; @@ -193,6 +203,8 @@ ExchangeFunctionalHamiltonian::functional_from_settings(const libconfig::Setting return bind(functional_kaneyoshi, _1, double(settings["J0"]), double(settings["r0"]), double(settings["sigma"])); } else if (functional_name == "step") { return bind(functional_step, _1, double(settings["J0"]), double(settings["r_out"])); + } else if (functional_name == "random") { + return bind(functional_random, _1, double(settings["J0"]), double(settings["rc"]), double(settings["width"])); } else { throw runtime_error("unknown exchange functional: " + functional_name); } From 47f58a593891645231b18e6709ad801d2ac7435a Mon Sep 17 00:00:00 2001 From: mkameda <39291377+mkameda@users.noreply.github.com> Date: Wed, 24 Nov 2021 11:31:57 +0900 Subject: [PATCH 57/80] [Fixed-2] Added random potential --- src/jams/hamiltonian/exchange_functional.h | 1 + 1 file changed, 1 insertion(+) diff --git a/src/jams/hamiltonian/exchange_functional.h b/src/jams/hamiltonian/exchange_functional.h index decbc3e5..dedf6c48 100644 --- a/src/jams/hamiltonian/exchange_functional.h +++ b/src/jams/hamiltonian/exchange_functional.h @@ -25,6 +25,7 @@ class ExchangeFunctionalHamiltonian : public SparseInteractionHamiltonian { static double functional_gaussian_multi(double rij, double J0, double r0, double lengthscale, double J0_2, double r0_2, double lengthscale_2, double J0_3, double r0_3, double lengthscale_3); static double functional_kaneyoshi(double rij, double J0, double r0, double lengthscale); static double functional_step(double rij, double J0, double r_out); + static double functional_random(const double rij, const double J0, const double rc, const double width); double radius_cutoff_; }; From 9c614434abec47369e49fa20612751e097aa7ed7 Mon Sep 17 00:00:00 2001 From: mkameda <39291377+mkameda@users.noreply.github.com> Date: Wed, 24 Nov 2021 11:34:58 +0900 Subject: [PATCH 58/80] [Fixed-3] Added random potential --- src/jams/hamiltonian/exchange_functional.cc | 1 + 1 file changed, 1 insertion(+) diff --git a/src/jams/hamiltonian/exchange_functional.cc b/src/jams/hamiltonian/exchange_functional.cc index 4dcb547a..0821e028 100644 --- a/src/jams/hamiltonian/exchange_functional.cc +++ b/src/jams/hamiltonian/exchange_functional.cc @@ -1,5 +1,6 @@ #include #include +#include #include "jams/helpers/output.h" #include "jams/hamiltonian/exchange_functional.h" #include "jams/core/lattice.h" From c946314ecd48df23bdb054e617c1b8dab6448daa Mon Sep 17 00:00:00 2001 From: mkameda <39291377+mkameda@users.noreply.github.com> Date: Wed, 24 Nov 2021 11:55:19 +0900 Subject: [PATCH 59/80] [Fixed-4] Added random potential --- src/jams/hamiltonian/exchange_functional.cc | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/jams/hamiltonian/exchange_functional.cc b/src/jams/hamiltonian/exchange_functional.cc index 0821e028..c332d635 100644 --- a/src/jams/hamiltonian/exchange_functional.cc +++ b/src/jams/hamiltonian/exchange_functional.cc @@ -176,10 +176,12 @@ double ExchangeFunctionalHamiltonian::functional_gaussian_multi(const double rij } double ExchangeFunctionalHamiltonian::functional_random(const double rij, const double J0, const double rc, const double width){ - std::uniform_real_distribution<> rand_potential(0.0, width * J0); //0.0 < width < 1.0 //rc is cutoff for the random potential + std::uniform_real_distribution<> rand_potential(-width * J0, width * J0); //-1.0 < width < 1.0 std::mt19937 rand_src(12345); //seed=12345 - if(rij < rc) { - return rand_potential(rand_src); + double rand = rand_potential(rand_src); + if(rij < r_out) { + std::cout << "random potential at r = " << rij << " is " << rand << << std::endl; + return J0 + rand; } else{ return 0.0; } From d2e27965bf848b2fb81fc697c5d47d0564ecda0a Mon Sep 17 00:00:00 2001 From: mkameda <39291377+mkameda@users.noreply.github.com> Date: Wed, 24 Nov 2021 11:55:49 +0900 Subject: [PATCH 60/80] [Fixed-5] Added random potential --- src/jams/hamiltonian/exchange_functional.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/jams/hamiltonian/exchange_functional.cc b/src/jams/hamiltonian/exchange_functional.cc index c332d635..8db1bf64 100644 --- a/src/jams/hamiltonian/exchange_functional.cc +++ b/src/jams/hamiltonian/exchange_functional.cc @@ -175,7 +175,7 @@ double ExchangeFunctionalHamiltonian::functional_gaussian_multi(const double rij return J0 * exp(-pow2(rij - r0)/(2 * pow2(sigma))) + J0_2 * exp(-pow2(rij - r0_2)/(2 * pow2(sigma_2))) + J0_3 * exp(-pow2(rij - r0_3)/(2 * pow2(sigma_3))); } -double ExchangeFunctionalHamiltonian::functional_random(const double rij, const double J0, const double rc, const double width){ +double ExchangeFunctionalHamiltonian::functional_random(const double rij, const double J0, const double r_out, const double width){ std::uniform_real_distribution<> rand_potential(-width * J0, width * J0); //-1.0 < width < 1.0 std::mt19937 rand_src(12345); //seed=12345 double rand = rand_potential(rand_src); @@ -207,7 +207,7 @@ ExchangeFunctionalHamiltonian::functional_from_settings(const libconfig::Setting } else if (functional_name == "step") { return bind(functional_step, _1, double(settings["J0"]), double(settings["r_out"])); } else if (functional_name == "random") { - return bind(functional_random, _1, double(settings["J0"]), double(settings["rc"]), double(settings["width"])); + return bind(functional_random, _1, double(settings["J0"]), double(settings["r_out"]), double(settings["width"])); } else { throw runtime_error("unknown exchange functional: " + functional_name); } From 294017d7aa240b29cc34eedc6a4dbb6a45c68c6e Mon Sep 17 00:00:00 2001 From: mkameda <39291377+mkameda@users.noreply.github.com> Date: Wed, 24 Nov 2021 11:57:39 +0900 Subject: [PATCH 61/80] [Fixed-6] Added random potential --- src/jams/hamiltonian/exchange_functional.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/jams/hamiltonian/exchange_functional.h b/src/jams/hamiltonian/exchange_functional.h index dedf6c48..ae531a03 100644 --- a/src/jams/hamiltonian/exchange_functional.h +++ b/src/jams/hamiltonian/exchange_functional.h @@ -25,7 +25,7 @@ class ExchangeFunctionalHamiltonian : public SparseInteractionHamiltonian { static double functional_gaussian_multi(double rij, double J0, double r0, double lengthscale, double J0_2, double r0_2, double lengthscale_2, double J0_3, double r0_3, double lengthscale_3); static double functional_kaneyoshi(double rij, double J0, double r0, double lengthscale); static double functional_step(double rij, double J0, double r_out); - static double functional_random(const double rij, const double J0, const double rc, const double width); + static double functional_random(const double rij, const double J0, const double r_out, const double width); double radius_cutoff_; }; From fe86feff870e6d2a036f87e2267ab3ce0bead0be Mon Sep 17 00:00:00 2001 From: mkameda <39291377+mkameda@users.noreply.github.com> Date: Wed, 24 Nov 2021 11:58:48 +0900 Subject: [PATCH 62/80] [Fixed-7] Added random potential --- src/jams/hamiltonian/exchange_functional.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/jams/hamiltonian/exchange_functional.cc b/src/jams/hamiltonian/exchange_functional.cc index 8db1bf64..3595f429 100644 --- a/src/jams/hamiltonian/exchange_functional.cc +++ b/src/jams/hamiltonian/exchange_functional.cc @@ -180,7 +180,7 @@ double ExchangeFunctionalHamiltonian::functional_random(const double rij, const std::mt19937 rand_src(12345); //seed=12345 double rand = rand_potential(rand_src); if(rij < r_out) { - std::cout << "random potential at r = " << rij << " is " << rand << << std::endl; + std::cout << "random potential at r = " << rij << " is " << rand << std::endl; return J0 + rand; } else{ return 0.0; From 1f8fc7140c0bf7d0ea3e50862e3cb5e2c89b018f Mon Sep 17 00:00:00 2001 From: mkameda <39291377+mkameda@users.noreply.github.com> Date: Wed, 24 Nov 2021 12:20:59 +0900 Subject: [PATCH 63/80] [Fixed-8] Added random potential --- src/jams/hamiltonian/exchange_functional.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/jams/hamiltonian/exchange_functional.cc b/src/jams/hamiltonian/exchange_functional.cc index 3595f429..b9a09cf4 100644 --- a/src/jams/hamiltonian/exchange_functional.cc +++ b/src/jams/hamiltonian/exchange_functional.cc @@ -180,7 +180,7 @@ double ExchangeFunctionalHamiltonian::functional_random(const double rij, const std::mt19937 rand_src(12345); //seed=12345 double rand = rand_potential(rand_src); if(rij < r_out) { - std::cout << "random potential at r = " << rij << " is " << rand << std::endl; + std::cout << "random potential at r = " << jams::fmt::decimal << rij << " is " << jams::fmt::sci << rand << std::endl; return J0 + rand; } else{ return 0.0; From bed062112596cb326e8ef818ba3a65739948e86e Mon Sep 17 00:00:00 2001 From: mkameda <39291377+mkameda@users.noreply.github.com> Date: Wed, 24 Nov 2021 12:31:51 +0900 Subject: [PATCH 64/80] [Fixed-9] Added random potential --- src/jams/hamiltonian/exchange_functional.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/jams/hamiltonian/exchange_functional.cc b/src/jams/hamiltonian/exchange_functional.cc index b9a09cf4..e7cc7c5a 100644 --- a/src/jams/hamiltonian/exchange_functional.cc +++ b/src/jams/hamiltonian/exchange_functional.cc @@ -180,7 +180,7 @@ double ExchangeFunctionalHamiltonian::functional_random(const double rij, const std::mt19937 rand_src(12345); //seed=12345 double rand = rand_potential(rand_src); if(rij < r_out) { - std::cout << "random potential at r = " << jams::fmt::decimal << rij << " is " << jams::fmt::sci << rand << std::endl; +// std::cout << "random potential at r = " << jams::fmt::decimal << rij << " is " << jams::fmt::sci << rand << std::endl; return J0 + rand; } else{ return 0.0; From 5ff4f3703d0ac83444eee13cf010730939f83a40 Mon Sep 17 00:00:00 2001 From: mkameda <39291377+mkameda@users.noreply.github.com> Date: Wed, 24 Nov 2021 14:54:56 +0900 Subject: [PATCH 65/80] [Fixed-10] Added random potential --- src/jams/hamiltonian/exchange_functional.cc | 5 ++--- src/jams/hamiltonian/exchange_functional.h | 2 ++ 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/src/jams/hamiltonian/exchange_functional.cc b/src/jams/hamiltonian/exchange_functional.cc index e7cc7c5a..f5062825 100644 --- a/src/jams/hamiltonian/exchange_functional.cc +++ b/src/jams/hamiltonian/exchange_functional.cc @@ -177,10 +177,9 @@ double ExchangeFunctionalHamiltonian::functional_gaussian_multi(const double rij double ExchangeFunctionalHamiltonian::functional_random(const double rij, const double J0, const double r_out, const double width){ std::uniform_real_distribution<> rand_potential(-width * J0, width * J0); //-1.0 < width < 1.0 - std::mt19937 rand_src(12345); //seed=12345 - double rand = rand_potential(rand_src); + double rand = rand_potential(rand_src_()); if(rij < r_out) { -// std::cout << "random potential at r = " << jams::fmt::decimal << rij << " is " << jams::fmt::sci << rand << std::endl; + std::cout << "random potential at r = " << jams::fmt::decimal << rij << " is " << jams::fmt::sci << rand << std::endl; return J0 + rand; } else{ return 0.0; diff --git a/src/jams/hamiltonian/exchange_functional.h b/src/jams/hamiltonian/exchange_functional.h index ae531a03..abc7cf9b 100644 --- a/src/jams/hamiltonian/exchange_functional.h +++ b/src/jams/hamiltonian/exchange_functional.h @@ -27,6 +27,8 @@ class ExchangeFunctionalHamiltonian : public SparseInteractionHamiltonian { static double functional_step(double rij, double J0, double r_out); static double functional_random(const double rij, const double J0, const double r_out, const double width); + std::mt19937 rand_src_(12345); //seed=12345 + double radius_cutoff_; }; From a16f37279234c6632b1319fc65071ef0fd24f772 Mon Sep 17 00:00:00 2001 From: mkameda <39291377+mkameda@users.noreply.github.com> Date: Wed, 24 Nov 2021 14:58:54 +0900 Subject: [PATCH 66/80] [Fixed-11] Added random potential --- src/jams/hamiltonian/exchange_functional.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/jams/hamiltonian/exchange_functional.cc b/src/jams/hamiltonian/exchange_functional.cc index f5062825..c8d7be16 100644 --- a/src/jams/hamiltonian/exchange_functional.cc +++ b/src/jams/hamiltonian/exchange_functional.cc @@ -177,7 +177,7 @@ double ExchangeFunctionalHamiltonian::functional_gaussian_multi(const double rij double ExchangeFunctionalHamiltonian::functional_random(const double rij, const double J0, const double r_out, const double width){ std::uniform_real_distribution<> rand_potential(-width * J0, width * J0); //-1.0 < width < 1.0 - double rand = rand_potential(rand_src_()); + double rand = rand_potential(rand_src_); if(rij < r_out) { std::cout << "random potential at r = " << jams::fmt::decimal << rij << " is " << jams::fmt::sci << rand << std::endl; return J0 + rand; From eff777ff7e91e4a2c89e30620a08f6cd51e50139 Mon Sep 17 00:00:00 2001 From: mkameda <39291377+mkameda@users.noreply.github.com> Date: Wed, 24 Nov 2021 15:02:51 +0900 Subject: [PATCH 67/80] [Fixed-12] Added random potential --- src/jams/hamiltonian/exchange_functional.cc | 4 ++-- src/jams/hamiltonian/exchange_functional.h | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/jams/hamiltonian/exchange_functional.cc b/src/jams/hamiltonian/exchange_functional.cc index c8d7be16..c57b3c6f 100644 --- a/src/jams/hamiltonian/exchange_functional.cc +++ b/src/jams/hamiltonian/exchange_functional.cc @@ -177,9 +177,9 @@ double ExchangeFunctionalHamiltonian::functional_gaussian_multi(const double rij double ExchangeFunctionalHamiltonian::functional_random(const double rij, const double J0, const double r_out, const double width){ std::uniform_real_distribution<> rand_potential(-width * J0, width * J0); //-1.0 < width < 1.0 - double rand = rand_potential(rand_src_); + double rand = rand_potential(random_generator_); if(rij < r_out) { - std::cout << "random potential at r = " << jams::fmt::decimal << rij << " is " << jams::fmt::sci << rand << std::endl; + std::cout << "random potential at r = " << rij << " is " << jams::fmt::sci << rand << std::endl; return J0 + rand; } else{ return 0.0; diff --git a/src/jams/hamiltonian/exchange_functional.h b/src/jams/hamiltonian/exchange_functional.h index abc7cf9b..9621cd9e 100644 --- a/src/jams/hamiltonian/exchange_functional.h +++ b/src/jams/hamiltonian/exchange_functional.h @@ -27,7 +27,7 @@ class ExchangeFunctionalHamiltonian : public SparseInteractionHamiltonian { static double functional_step(double rij, double J0, double r_out); static double functional_random(const double rij, const double J0, const double r_out, const double width); - std::mt19937 rand_src_(12345); //seed=12345 + pcg32_k1024 random_generator_ = pcg_extras::seed_seq_from(jams::random_generator()); double radius_cutoff_; }; From d3ab6e2a4c9d4ccd8aeb5c66d681c575b8abbfdb Mon Sep 17 00:00:00 2001 From: mkameda <39291377+mkameda@users.noreply.github.com> Date: Wed, 24 Nov 2021 15:15:20 +0900 Subject: [PATCH 68/80] [Fixed-14] Added random potential --- src/jams/hamiltonian/exchange_functional.cc | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/jams/hamiltonian/exchange_functional.cc b/src/jams/hamiltonian/exchange_functional.cc index c57b3c6f..9d4ae328 100644 --- a/src/jams/hamiltonian/exchange_functional.cc +++ b/src/jams/hamiltonian/exchange_functional.cc @@ -176,11 +176,12 @@ double ExchangeFunctionalHamiltonian::functional_gaussian_multi(const double rij } double ExchangeFunctionalHamiltonian::functional_random(const double rij, const double J0, const double r_out, const double width){ - std::uniform_real_distribution<> rand_potential(-width * J0, width * J0); //-1.0 < width < 1.0 +// std::mt19937 rand_src(12345); //seed=12345 + std::uniform_real_distribution<> rand_potential(-width, width); double rand = rand_potential(random_generator_); if(rij < r_out) { std::cout << "random potential at r = " << rij << " is " << jams::fmt::sci << rand << std::endl; - return J0 + rand; + return J0 * (1 + rand); } else{ return 0.0; } From 4ef29474d90fe4cd4b33832a4b02f37edd074c16 Mon Sep 17 00:00:00 2001 From: mkameda <39291377+mkameda@users.noreply.github.com> Date: Wed, 24 Nov 2021 15:34:29 +0900 Subject: [PATCH 69/80] [Fixed-15] Added random potential --- src/jams/hamiltonian/exchange_functional.cc | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/jams/hamiltonian/exchange_functional.cc b/src/jams/hamiltonian/exchange_functional.cc index 9d4ae328..332dcc0d 100644 --- a/src/jams/hamiltonian/exchange_functional.cc +++ b/src/jams/hamiltonian/exchange_functional.cc @@ -176,9 +176,9 @@ double ExchangeFunctionalHamiltonian::functional_gaussian_multi(const double rij } double ExchangeFunctionalHamiltonian::functional_random(const double rij, const double J0, const double r_out, const double width){ -// std::mt19937 rand_src(12345); //seed=12345 + std::mt19937 rand_src(12345); //seed=12345 std::uniform_real_distribution<> rand_potential(-width, width); - double rand = rand_potential(random_generator_); + double rand = rand_potential(rand_src); if(rij < r_out) { std::cout << "random potential at r = " << rij << " is " << jams::fmt::sci << rand << std::endl; return J0 * (1 + rand); @@ -213,8 +213,7 @@ ExchangeFunctionalHamiltonian::functional_from_settings(const libconfig::Setting } } -void -ExchangeFunctionalHamiltonian::output_exchange_functional(std::ostream &os, +void ExchangeFunctionalHamiltonian::output_exchange_functional(std::ostream &os, const ExchangeFunctionalHamiltonian::ExchangeFunctional& functional, double r_cutoff, double delta_r) { os << "radius exchange" << "\n"; double r = 0.0; From c0a8426f0d8209edda1bc333edc337fdb3f4af20 Mon Sep 17 00:00:00 2001 From: mkameda <39291377+mkameda@users.noreply.github.com> Date: Wed, 24 Nov 2021 15:38:15 +0900 Subject: [PATCH 70/80] [Fixed-16] Added random potential --- src/jams/hamiltonian/exchange_functional.h | 1 + 1 file changed, 1 insertion(+) diff --git a/src/jams/hamiltonian/exchange_functional.h b/src/jams/hamiltonian/exchange_functional.h index 9621cd9e..637448fd 100644 --- a/src/jams/hamiltonian/exchange_functional.h +++ b/src/jams/hamiltonian/exchange_functional.h @@ -4,6 +4,7 @@ #include #include +#include #include "jams/core/hamiltonian.h" #include "jams/core/interactions.h" From bc4f4bbce7bcde2ba49e06a780c960d99b12297a Mon Sep 17 00:00:00 2001 From: mkameda <39291377+mkameda@users.noreply.github.com> Date: Wed, 24 Nov 2021 15:40:08 +0900 Subject: [PATCH 71/80] [Fixed-17] Added random potential --- src/jams/hamiltonian/exchange_functional.h | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/jams/hamiltonian/exchange_functional.h b/src/jams/hamiltonian/exchange_functional.h index 637448fd..c15b125c 100644 --- a/src/jams/hamiltonian/exchange_functional.h +++ b/src/jams/hamiltonian/exchange_functional.h @@ -5,7 +5,11 @@ #include #include +#include +#include +#include "jams/helpers/random.h" +#include "jblib/containers/array.h" #include "jams/core/hamiltonian.h" #include "jams/core/interactions.h" #include "jams/containers/sparse_matrix.h" From af5f5ab07af9c345c427b825c4d56ee8a54dd9b8 Mon Sep 17 00:00:00 2001 From: mkameda <39291377+mkameda@users.noreply.github.com> Date: Wed, 24 Nov 2021 15:42:40 +0900 Subject: [PATCH 72/80] [Fixed-18] Added random potential --- src/jams/hamiltonian/exchange_functional.h | 1 - 1 file changed, 1 deletion(-) diff --git a/src/jams/hamiltonian/exchange_functional.h b/src/jams/hamiltonian/exchange_functional.h index c15b125c..c94da583 100644 --- a/src/jams/hamiltonian/exchange_functional.h +++ b/src/jams/hamiltonian/exchange_functional.h @@ -4,7 +4,6 @@ #include #include -#include #include #include From 5e85d491077c9382cc15af31793e305dbb01c14d Mon Sep 17 00:00:00 2001 From: mkameda <39291377+mkameda@users.noreply.github.com> Date: Wed, 24 Nov 2021 15:44:21 +0900 Subject: [PATCH 73/80] [Fixed-19] Added random potential --- src/jams/hamiltonian/exchange_functional.h | 1 - 1 file changed, 1 deletion(-) diff --git a/src/jams/hamiltonian/exchange_functional.h b/src/jams/hamiltonian/exchange_functional.h index c94da583..80809f9f 100644 --- a/src/jams/hamiltonian/exchange_functional.h +++ b/src/jams/hamiltonian/exchange_functional.h @@ -8,7 +8,6 @@ #include #include "jams/helpers/random.h" -#include "jblib/containers/array.h" #include "jams/core/hamiltonian.h" #include "jams/core/interactions.h" #include "jams/containers/sparse_matrix.h" From a6feb09d27b92b8ee5f0a505308c1609b1ff735c Mon Sep 17 00:00:00 2001 From: mkameda <39291377+mkameda@users.noreply.github.com> Date: Wed, 24 Nov 2021 15:45:38 +0900 Subject: [PATCH 74/80] [Fixed-20] Added random potential --- src/jams/hamiltonian/exchange_functional.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/jams/hamiltonian/exchange_functional.h b/src/jams/hamiltonian/exchange_functional.h index 80809f9f..dd9e42bd 100644 --- a/src/jams/hamiltonian/exchange_functional.h +++ b/src/jams/hamiltonian/exchange_functional.h @@ -4,7 +4,7 @@ #include #include -#include +#include #include #include "jams/helpers/random.h" From 3e9018c6f595fa9e06da4237ab12754c5e302053 Mon Sep 17 00:00:00 2001 From: mkameda <39291377+mkameda@users.noreply.github.com> Date: Wed, 24 Nov 2021 15:50:17 +0900 Subject: [PATCH 75/80] [Fixed-21] Added random potential --- src/jams/hamiltonian/exchange_functional.cc | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/jams/hamiltonian/exchange_functional.cc b/src/jams/hamiltonian/exchange_functional.cc index 332dcc0d..87431f4a 100644 --- a/src/jams/hamiltonian/exchange_functional.cc +++ b/src/jams/hamiltonian/exchange_functional.cc @@ -5,6 +5,7 @@ #include "jams/hamiltonian/exchange_functional.h" #include "jams/core/lattice.h" #include "jams/core/globals.h" +#include "exchange_functional.h" using namespace std; @@ -176,9 +177,9 @@ double ExchangeFunctionalHamiltonian::functional_gaussian_multi(const double rij } double ExchangeFunctionalHamiltonian::functional_random(const double rij, const double J0, const double r_out, const double width){ - std::mt19937 rand_src(12345); //seed=12345 +// std::mt19937 rand_src(12345); //seed=12345 std::uniform_real_distribution<> rand_potential(-width, width); - double rand = rand_potential(rand_src); + double rand = rand_potential(random_generator_); if(rij < r_out) { std::cout << "random potential at r = " << rij << " is " << jams::fmt::sci << rand << std::endl; return J0 * (1 + rand); From 7bba8af28b9dc54235dab19a93159c987a632648 Mon Sep 17 00:00:00 2001 From: mkameda <39291377+mkameda@users.noreply.github.com> Date: Wed, 24 Nov 2021 15:53:36 +0900 Subject: [PATCH 76/80] [Fixed-22] Added random potential --- src/jams/hamiltonian/exchange_functional.cc | 2 +- src/jams/hamiltonian/exchange_functional.h | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/jams/hamiltonian/exchange_functional.cc b/src/jams/hamiltonian/exchange_functional.cc index 87431f4a..6025e766 100644 --- a/src/jams/hamiltonian/exchange_functional.cc +++ b/src/jams/hamiltonian/exchange_functional.cc @@ -179,7 +179,7 @@ double ExchangeFunctionalHamiltonian::functional_gaussian_multi(const double rij double ExchangeFunctionalHamiltonian::functional_random(const double rij, const double J0, const double r_out, const double width){ // std::mt19937 rand_src(12345); //seed=12345 std::uniform_real_distribution<> rand_potential(-width, width); - double rand = rand_potential(random_generator_); + double rand = rand_potential(rand_src_); if(rij < r_out) { std::cout << "random potential at r = " << rij << " is " << jams::fmt::sci << rand << std::endl; return J0 * (1 + rand); diff --git a/src/jams/hamiltonian/exchange_functional.h b/src/jams/hamiltonian/exchange_functional.h index dd9e42bd..0d1fabc1 100644 --- a/src/jams/hamiltonian/exchange_functional.h +++ b/src/jams/hamiltonian/exchange_functional.h @@ -31,7 +31,7 @@ class ExchangeFunctionalHamiltonian : public SparseInteractionHamiltonian { static double functional_random(const double rij, const double J0, const double r_out, const double width); pcg32_k1024 random_generator_ = pcg_extras::seed_seq_from(jams::random_generator()); - + std::mt19937 rand_src_(12345); double radius_cutoff_; }; From 185597638886ed5d06f016ece3bcd3a79cb2f378 Mon Sep 17 00:00:00 2001 From: mkameda <39291377+mkameda@users.noreply.github.com> Date: Wed, 24 Nov 2021 15:56:26 +0900 Subject: [PATCH 77/80] [Fixed-23] Added random potential --- src/jams/hamiltonian/exchange_functional.h | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/jams/hamiltonian/exchange_functional.h b/src/jams/hamiltonian/exchange_functional.h index 0d1fabc1..4ab3bcc5 100644 --- a/src/jams/hamiltonian/exchange_functional.h +++ b/src/jams/hamiltonian/exchange_functional.h @@ -31,7 +31,8 @@ class ExchangeFunctionalHamiltonian : public SparseInteractionHamiltonian { static double functional_random(const double rij, const double J0, const double r_out, const double width); pcg32_k1024 random_generator_ = pcg_extras::seed_seq_from(jams::random_generator()); - std::mt19937 rand_src_(12345); + int seed = 12345; + std::mt19937 rand_src_(seed); double radius_cutoff_; }; From 649b47f2403a98952b9e10384e621bae794c24fa Mon Sep 17 00:00:00 2001 From: mkameda <39291377+mkameda@users.noreply.github.com> Date: Wed, 24 Nov 2021 16:11:11 +0900 Subject: [PATCH 78/80] [Fixed-24] Added random potential --- src/jams/hamiltonian/exchange_functional.cc | 4 ++-- src/jams/hamiltonian/exchange_functional.h | 2 -- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/src/jams/hamiltonian/exchange_functional.cc b/src/jams/hamiltonian/exchange_functional.cc index 6025e766..53fa5bf9 100644 --- a/src/jams/hamiltonian/exchange_functional.cc +++ b/src/jams/hamiltonian/exchange_functional.cc @@ -177,9 +177,9 @@ double ExchangeFunctionalHamiltonian::functional_gaussian_multi(const double rij } double ExchangeFunctionalHamiltonian::functional_random(const double rij, const double J0, const double r_out, const double width){ -// std::mt19937 rand_src(12345); //seed=12345 + std::mt19937 rand_src(rij); //seed=rij std::uniform_real_distribution<> rand_potential(-width, width); - double rand = rand_potential(rand_src_); + double rand = rand_potential(rand_src); if(rij < r_out) { std::cout << "random potential at r = " << rij << " is " << jams::fmt::sci << rand << std::endl; return J0 * (1 + rand); diff --git a/src/jams/hamiltonian/exchange_functional.h b/src/jams/hamiltonian/exchange_functional.h index 4ab3bcc5..0a0a348e 100644 --- a/src/jams/hamiltonian/exchange_functional.h +++ b/src/jams/hamiltonian/exchange_functional.h @@ -31,8 +31,6 @@ class ExchangeFunctionalHamiltonian : public SparseInteractionHamiltonian { static double functional_random(const double rij, const double J0, const double r_out, const double width); pcg32_k1024 random_generator_ = pcg_extras::seed_seq_from(jams::random_generator()); - int seed = 12345; - std::mt19937 rand_src_(seed); double radius_cutoff_; }; From 4243b8cabe7abd38141da4f9e7317d53b471646f Mon Sep 17 00:00:00 2001 From: mkameda <39291377+mkameda@users.noreply.github.com> Date: Wed, 24 Nov 2021 16:45:33 +0900 Subject: [PATCH 79/80] [Fixed-25] Added random potential --- src/jams/hamiltonian/exchange_functional.cc | 37 ++++++++++----------- src/jams/hamiltonian/exchange_functional.h | 2 -- 2 files changed, 17 insertions(+), 22 deletions(-) diff --git a/src/jams/hamiltonian/exchange_functional.cc b/src/jams/hamiltonian/exchange_functional.cc index 53fa5bf9..1b2ad253 100644 --- a/src/jams/hamiltonian/exchange_functional.cc +++ b/src/jams/hamiltonian/exchange_functional.cc @@ -70,7 +70,11 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se k(n) = kvector * n * (kmax / num_k); // cout << "n = " << n << ", kspace_path_(n) = " << k(n) << endl; } - // --- for crystal limit spectrum --- + + if(settings.exists("random")){ + std::mt19937 rand_src(12345); //seed=12345 + } + for (auto i = 0; i < globals::num_spins; ++i) { nbrs.clear(); ::lattice->atom_neighbours(i, radius_cutoff_, nbrs); @@ -85,13 +89,23 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se cout << "central coordinate 2 = ( " << rspace_displacement2_(central_site_indx) << " )" << endl; } - for (const auto& nbr : nbrs) { + if(settings.exists("random")){ + std::uniform_real_distribution<> rand_potential(-width, width); + double rand = rand_potential(rand_src); + } + + for (const auto& nbr : nbrs) { const auto j = nbr.id; if (i == j) { continue; } + const auto rij = norm(::lattice->displacement(i, j)); - this->insert_interaction_scalar(i, j, input_unit_conversion_ * exchange_functional(rij)); + if(settings.exists("random")){ + this->insert_interaction_scalar(i, j, input_unit_conversion_ * (exchange_functional(rij) * (1 + rand))); + } else{ + this->insert_interaction_scalar(i, j, input_unit_conversion_ * exchange_functional(rij)); + } counter++; // --- for crystal limit spectrum --- if(rij < rij_max && rij > rij_min){ @@ -127,9 +141,6 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se cout << " average interactions per spin (kr != 0) " << jams::fmt::decimal << counter2 / double(globals::num_spins)/num_k << "\n"; // --- for crystal limit spectrum --- for (auto m = 0; m < spectrum_crystal_limit.size(); m++) { -// cout << " spectrum_crystal_limit (" << m << ") = " << spectrum_crystal_limit[m] << "\n"; -// cout << " real (" << m << ") = " << spectrum_crystal_limit[m].real() << "\n"; -// cout << " imag (" << m << ") = " << spectrum_crystal_limit[m].imag() << "\n"; outfile3 << std::setw(20) << k(m)[0] << "\t"; outfile3 << std::setw(20) << k(m)[1] << "\t"; outfile3 << std::setw(20) << k(m)[2] << "\t"; @@ -176,18 +187,6 @@ double ExchangeFunctionalHamiltonian::functional_gaussian_multi(const double rij return J0 * exp(-pow2(rij - r0)/(2 * pow2(sigma))) + J0_2 * exp(-pow2(rij - r0_2)/(2 * pow2(sigma_2))) + J0_3 * exp(-pow2(rij - r0_3)/(2 * pow2(sigma_3))); } -double ExchangeFunctionalHamiltonian::functional_random(const double rij, const double J0, const double r_out, const double width){ - std::mt19937 rand_src(rij); //seed=rij - std::uniform_real_distribution<> rand_potential(-width, width); - double rand = rand_potential(rand_src); - if(rij < r_out) { - std::cout << "random potential at r = " << rij << " is " << jams::fmt::sci << rand << std::endl; - return J0 * (1 + rand); - } else{ - return 0.0; - } -} - ExchangeFunctionalHamiltonian::ExchangeFunctional ExchangeFunctionalHamiltonian::functional_from_settings(const libconfig::Setting &settings) { using namespace std::placeholders; @@ -207,8 +206,6 @@ ExchangeFunctionalHamiltonian::functional_from_settings(const libconfig::Setting return bind(functional_kaneyoshi, _1, double(settings["J0"]), double(settings["r0"]), double(settings["sigma"])); } else if (functional_name == "step") { return bind(functional_step, _1, double(settings["J0"]), double(settings["r_out"])); - } else if (functional_name == "random") { - return bind(functional_random, _1, double(settings["J0"]), double(settings["r_out"]), double(settings["width"])); } else { throw runtime_error("unknown exchange functional: " + functional_name); } diff --git a/src/jams/hamiltonian/exchange_functional.h b/src/jams/hamiltonian/exchange_functional.h index 0a0a348e..c6a6cb1d 100644 --- a/src/jams/hamiltonian/exchange_functional.h +++ b/src/jams/hamiltonian/exchange_functional.h @@ -28,9 +28,7 @@ class ExchangeFunctionalHamiltonian : public SparseInteractionHamiltonian { static double functional_gaussian_multi(double rij, double J0, double r0, double lengthscale, double J0_2, double r0_2, double lengthscale_2, double J0_3, double r0_3, double lengthscale_3); static double functional_kaneyoshi(double rij, double J0, double r0, double lengthscale); static double functional_step(double rij, double J0, double r_out); - static double functional_random(const double rij, const double J0, const double r_out, const double width); - pcg32_k1024 random_generator_ = pcg_extras::seed_seq_from(jams::random_generator()); double radius_cutoff_; }; From 0e6ac070172b02a27683dacc29cfa7cd8c0e8cb7 Mon Sep 17 00:00:00 2001 From: mkameda <39291377+mkameda@users.noreply.github.com> Date: Wed, 24 Nov 2021 16:48:19 +0900 Subject: [PATCH 80/80] [Fixed-26] Added random potential --- src/jams/hamiltonian/exchange_functional.cc | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/jams/hamiltonian/exchange_functional.cc b/src/jams/hamiltonian/exchange_functional.cc index 1b2ad253..32511915 100644 --- a/src/jams/hamiltonian/exchange_functional.cc +++ b/src/jams/hamiltonian/exchange_functional.cc @@ -71,9 +71,7 @@ ExchangeFunctionalHamiltonian::ExchangeFunctionalHamiltonian(const libconfig::Se // cout << "n = " << n << ", kspace_path_(n) = " << k(n) << endl; } - if(settings.exists("random")){ - std::mt19937 rand_src(12345); //seed=12345 - } + std::mt19937 rand_src(12345); //seed=12345 for (auto i = 0; i < globals::num_spins; ++i) { nbrs.clear();