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/hamiltonian/exchange_functional.cc b/src/jams/hamiltonian/exchange_functional.cc index e767affc..32511915 100644 --- a/src/jams/hamiltonian/exchange_functional.cc +++ b/src/jams/hamiltonian/exchange_functional.cc @@ -1,9 +1,11 @@ #include #include +#include #include "jams/helpers/output.h" #include "jams/hamiltonian/exchange_functional.h" #include "jams/core/lattice.h" #include "jams/core/globals.h" +#include "exchange_functional.h" using namespace std; @@ -13,6 +15,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,25 +24,134 @@ 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) << "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)" << "\t"; + 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" << "\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; vector nbrs; + // --- 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}); + 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); + 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; + 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); + for (auto n = 0; n < k.size(); ++n) { + k(n) = kvector * n * (kmax / num_k); +// cout << "n = " << n << ", kspace_path_(n) = " << k(n) << endl; + } + + 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); - for (const auto& nbr : nbrs) { + 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->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; + cout << "central coordinate = ( " << rspace_displacement_(central_site_indx) << " )" << endl; + cout << "central coordinate 2 = ( " << rspace_displacement2_(central_site_indx) << " )" << endl; + } + + 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)); + + const auto rij = norm(::lattice->displacement(i, j)); + 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){ + 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 == central_site){ + spectrum_crystal_limit_noave[kk] += tmp; + } + else if(i == central_site_indx){ + spectrum_crystal_limit_noave_central_cite[kk] += tmp; + cout << "NOW!" << endl; + } + counter2++; + } + } + } + // --- for crystal limit spectrum --- } } + 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"; 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)/num_k << "\n"; + // --- for crystal limit spectrum --- + for (auto m = 0; m < spectrum_crystal_limit.size(); m++) { + 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() << "\t"; + 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() << "\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 --- finalize(); } @@ -61,6 +173,18 @@ 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; +} + +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; @@ -74,15 +198,18 @@ 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"])); + 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); } } -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; diff --git a/src/jams/hamiltonian/exchange_functional.h b/src/jams/hamiltonian/exchange_functional.h index 7c5d9c5e..c6a6cb1d 100644 --- a/src/jams/hamiltonian/exchange_functional.h +++ b/src/jams/hamiltonian/exchange_functional.h @@ -4,7 +4,10 @@ #include #include +#include +#include +#include "jams/helpers/random.h" #include "jams/core/hamiltonian.h" #include "jams/core/interactions.h" #include "jams/containers/sparse_matrix.h" @@ -22,7 +25,9 @@ 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); double radius_cutoff_; }; diff --git a/src/jams/hamiltonian/sparse_interaction.cc b/src/jams/hamiltonian/sparse_interaction.cc index ef1bcfd7..38616ce9 100644 --- a/src/jams/hamiltonian/sparse_interaction.cc +++ b/src/jams/hamiltonian/sparse_interaction.cc @@ -19,9 +19,8 @@ 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); +// std::cout << "i= " << i << " j= " << j << " : row = " << 3 * i + m << ", colum = " << 3 * j + m << ", value [J] = " << value/input_unit_conversion_ << std::endl; } } diff --git a/src/jams/helpers/neutrons.cc b/src/jams/helpers/neutrons.cc index da1a121e..cae33778 100644 --- a/src/jams/helpers/neutrons.cc +++ b/src/jams/helpers/neutrons.cc @@ -9,12 +9,22 @@ 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); + +// 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 1dc4a43f..e0f68e03 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); @@ -22,33 +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; 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/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 83c7b7ff..dd80463e 100644 --- a/src/jams/monitors/neutron_scattering_no_lattice.cc +++ b/src/jams/monitors/neutron_scattering_no_lattice.cc @@ -22,22 +22,29 @@ 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; + + 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"]); + evaluate_form_factor_ = true; } if (settings.exists("polarizations")) { 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( @@ -50,6 +57,24 @@ void NeutronScatteringNoLatticeMonitor::update(Solver *solver) { store_kspace_data_on_path(); periodogram_index_++; + if (output_dist_) { + 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(); @@ -58,6 +83,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()) { @@ -72,26 +98,36 @@ void NeutronScatteringNoLatticeMonitor::update(Solver *solver) { } void NeutronScatteringNoLatticeMonitor::configure_kspace_vectors(const libconfig::Setting &settings) { - kvector_ = jams::config_optional(settings, "kvector", kvector_); - - kspace_path_.resize(num_k_); - for (auto i = 0; i < kspace_path_.size(); ++i) { - kspace_path_(i) = kvector_ * i * (kmax_ / num_k_); + 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_+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_);//kmax shoul be in the unit of [2*pi/a(lattice constant)] + cout << "i = " << i << ", kspace_path_(i) = " << kspace_path_(i) << endl; } 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_(i) = lattice->displacement({0,0,0}, lattice->atom_position(i)); + //rspace_displacement_ is the position vector of i-th spin + //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)); } } -jams::MultiArray +//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, 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) { @@ -100,10 +136,18 @@ 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}) { + //(*) 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]); } } @@ -112,14 +156,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) { @@ -195,7 +239,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 +250,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 +287,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 +308,35 @@ 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); - if (norm(r) >= 0.5) continue; + //r should be normalized to get the length of position vector + 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; + //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); +// cout << "do_rspace_windowing_ = true, value = " << do_rspace_windowing_ << endl; } - + //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; } } @@ -351,3 +402,97 @@ jams::MultiArray NeutronScatteringNoLatticeMonitor::average_kspace_ti return average; } + +NeutronScatteringNoLatticeMonitor::~NeutronScatteringNoLatticeMonitor() { + if (output_dist_) { + using namespace std; + using namespace globals; + + const double delta_t = periodogram_props_.sample_time; + + string name4 = seedname + "_radial_dist.tsv"; + ofstream outfile4(name4.c_str()); + outfile4.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"; + + 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(); + + 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"; + } + + outfile4.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::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); + 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){ + 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]++; + if(abs(r_ij[2]) > (delta_r_ * 10.0)){ + count_z(m)++; + } + } + } + for (auto kk = 0; kk < dist.size(); kk++){ + dist[kk] /= num_spins; + count_z(kk) /= num_spins; +// cout << "dist[" << kk << "]/num_spins = " << dist[kk] << endl; +// cout << "dist_z[" << kk << "]/num_spins = " << count_z(kk) << 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"; + } + outfile5.close(); + 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 62e0b081..30acb269 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; @@ -34,7 +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_ = true; + bool do_rspace_windowing_; + bool output_dist_; jams::MultiArray rspace_displacement_; jams::MultiArray kspace_path_; jams::MultiArray kspace_spins_timeseries_; @@ -49,11 +51,21 @@ 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_; + 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(); + std::vector radial_distribution_function_z(); + double distance(const Vec3 &r_ij); + std::vector> time_correlation(unsigned i, unsigned j); }; #endif //JAMS_NEUTRON_SCATTERING_NO_LATTICE_H