Hotfix/exchange functional fixes

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;
+}

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);

src/jams/hamiltonian/exchange_functional.cc

@@ -1,9 +1,11 @@
 #include <functional>
 #include <fstream>
+#include <random>
 #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<double>(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<Atom> nbrs;
+  // --- for crystal limit spectrum ---
+  int num_k = jams::config_required<int>(settings, "num_k");
+  std::vector<complex<double>> spectrum_crystal_limit(num_k+1,{0.0,0.0});
+  std::vector<complex<double>> spectrum_crystal_limit2(num_k+1,{0.0,0.0});
+  std::vector<complex<double>> spectrum_crystal_limit_noave(num_k+1,{0.0,0.0});
+  std::vector<complex<double>> spectrum_crystal_limit_noave_central_cite(num_k+1,{0.0,0.0});
+  double kmax = jams::config_required<double>(settings, "kmax");
+  Vec3 kvector = jams::config_required<Vec3>(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<Vec3, 1> rspace_displacement_;
+  rspace_displacement_.resize(globals::s.size(0));
+  jams::MultiArray<Vec3, 1> rspace_displacement2_;
+  rspace_displacement2_.resize(globals::s.size(0));
+
+  jams::MultiArray<Vec3, 1> 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<double> tmp = { exchange_functional(rij)* (1.0-cos(kTwoPi*kr)),  exchange_functional(rij) * sin(kTwoPi*kr)};
+                  std::complex<double> 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;

src/jams/hamiltonian/exchange_functional.h

@@ -4,7 +4,10 @@
 #include <functional>
 
 #include <libconfig.h++>
+#include <pcg_random.hpp>
+#include <random>
 
+#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_;
 };

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;
   }
 }
 

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;
     }

src/jams/interface/config.h

@@ -9,6 +9,7 @@
 #include <jams/helpers/utils.h>
 #include <jams/core/interactions.h>
 #include "jams/core/types.h"
+#include <iostream>
 
 void config_patch(libconfig::Setting& orig, const libconfig::Setting& patch);
 
@@ -22,33 +23,34 @@ namespace jams {
         if (setting.exists(name)) {
           return config_required<T>(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<>

src/jams/interface/fft.h

@@ -23,6 +23,7 @@
 
 using Complex = std::complex<double>;
 
+//(*) where are these values used?
 namespace jams {
     struct PeriodogramProps {
         int length    = 1000;