Merge pull request #20 from maxscheurer/thole_damping

Thole damping

cppe/core/cppe_state.cc

@@ -50,7 +50,7 @@ void CppeState::calculate_static_energies_and_fields() {
   NuclearFields nfields(m_mol, m_potentials);
   m_nuc_fields = nfields.compute();
   // Multipole fields
-  MultipoleFields mul_fields(m_potentials);
+  MultipoleFields mul_fields(m_potentials, m_options);
   m_multipole_fields = mul_fields.compute();
 }
 

cppe/core/electric_fields.cc

@@ -6,10 +6,7 @@
 
 namespace libcppe {
 
-Eigen::VectorXd NuclearFields::compute(bool damp_core) {
-  if (damp_core) {
-    throw std::runtime_error("damping not implemented");
-  }
+Eigen::VectorXd NuclearFields::compute() {
   std::vector<Eigen::MatrixXi> Tk_coeffs = Tk_coefficients(5);
   Eigen::VectorXd nuc_fields             = Eigen::VectorXd::Zero(3 * m_n_polsites);
 #pragma omp parallel for firstprivate(Tk_coeffs)
@@ -29,10 +26,7 @@ Eigen::VectorXd NuclearFields::compute(bool damp_core) {
   return nuc_fields;
 }
 
-Eigen::VectorXd MultipoleFields::compute(bool damp) {
-  if (damp) {
-    throw std::runtime_error("damping not implemented");
-  }
+Eigen::VectorXd MultipoleFields::compute() {
   std::vector<Eigen::MatrixXi> Tk_coeffs = Tk_coefficients(5);
   Eigen::VectorXd mult_fields            = Eigen::VectorXd::Zero(3 * m_n_polsites);
 // Field at site of potential1 caused by all other sites (also non-polarizable
@@ -41,21 +35,40 @@ Eigen::VectorXd MultipoleFields::compute(bool damp) {
   for (size_t i = 0; i < m_n_polsites; i++) {
     size_t site_counter   = 3 * i;
     Potential& potential1 = m_polsites[i];
+    double alpha_i_isotropic =
+          potential1.get_polarizabilities()[0].get_isotropic_value();  // for damping
     for (size_t j = 0; j < m_potentials.size(); j++) {
       Potential& potential2 =
             m_potentials[j];  // all other multipoles create el. field at site i
       if (potential1.index == potential2.index) continue;
       if (potential1.excludes_site(potential2.index)) continue;
       Eigen::Vector3d diff =
             potential1.get_site_position() - potential2.get_site_position();
+
+      // Thole damping: potential2 needs to be polarizable
+      bool damp_enabled        = m_options.damp_multipole;
+      double alpha_j_isotropic = 0.0;
+      if (!potential2.is_polarizable()) {
+        damp_enabled = false;
+      } else {
+        alpha_j_isotropic =
+              potential2.get_polarizabilities()[0].get_isotropic_value();  // for damping
+      }
+
       // std::cout << "-- created by site " << potential2.index << std::endl;
       for (auto& mul : potential2.get_multipoles()) {
         // if (std::all_of(mul.get_values().begin(), mul.get_values().end(),
         //                 [](double v) { return std::abs(v) == 0.0; })) {
         //   continue;
         // }
-        Eigen::VectorXd Fi =
-              multipole_derivative(mul.m_k, 1, diff, mul.get_values_vec(), Tk_coeffs);
+        Eigen::VectorXd Fi;
+        if (damp_enabled) {
+          Fi = multipole_derivative(mul.m_k, 1, diff, mul.get_values_vec(), Tk_coeffs,
+                                    m_options.damping_factor_multipole, alpha_i_isotropic,
+                                    alpha_j_isotropic);
+        } else {
+          Fi = multipole_derivative(mul.m_k, 1, diff, mul.get_values_vec(), Tk_coeffs);
+        }
         mult_fields(site_counter) += Fi(0);
         mult_fields(site_counter + 1) += Fi(1);
         mult_fields(site_counter + 2) += Fi(2);
@@ -119,7 +132,15 @@ void InducedMoments::compute(const Eigen::VectorXd& total_fields,
           continue;
         }
         Eigen::Vector3d diff = pot2.get_site_position() - pot1.get_site_position();
-        Eigen::VectorXd T2   = Tk_tensor(2, diff, Tk_coeffs);
+        Eigen::VectorXd T2;
+        if (m_options.damp_induced) {
+          Polarizability& alpha_i = pot1.get_polarizabilities()[0];
+          Polarizability& alpha_j = pot2.get_polarizabilities()[0];
+          T2 = Tk_tensor(2, diff, Tk_coeffs, m_options.damping_factor_induced,
+                         alpha_i.get_isotropic_value(), alpha_j.get_isotropic_value());
+        } else {
+          T2 = Tk_tensor(2, diff, Tk_coeffs);
+        }
         Ftmp += smat_vec(T2, induced_moments.segment(m, 3), 1.0);
       }
       // keep value to calculate residual

cppe/core/electric_fields.hh

@@ -22,7 +22,7 @@ class ElectricFields {
     m_n_polsites = m_polsites.size();
   };
   ~ElectricFields(){};
-  virtual Eigen::VectorXd compute(bool damp = false) = 0;
+  virtual Eigen::VectorXd compute() = 0;
 };
 
 class NuclearFields : public ElectricFields {
@@ -32,13 +32,17 @@ class NuclearFields : public ElectricFields {
  public:
   NuclearFields(Molecule mol, std::vector<Potential> potentials)
         : ElectricFields(potentials), m_mol(mol){};
-  Eigen::VectorXd compute(bool damp = false);
+  Eigen::VectorXd compute();
 };
 
 class MultipoleFields : public ElectricFields {
  public:
-  MultipoleFields(std::vector<Potential> potentials) : ElectricFields(potentials){};
-  Eigen::VectorXd compute(bool damp = false);
+  MultipoleFields(std::vector<Potential> potentials, PeOptions options)
+        : ElectricFields(potentials), m_options(options){};
+  Eigen::VectorXd compute();
+
+ private:
+  PeOptions m_options;
 };
 
 class InducedMoments {

cppe/core/math.cc

@@ -17,17 +17,17 @@ Eigen::Vector3d smat_vec(const Eigen::VectorXd& mat, const Eigen::Vector3d& vec,
   return result;
 }
 
-// TODO: add option for damping
 Eigen::VectorXd Tk_tensor(int k, const Eigen::Vector3d& Rij,
-                          std::vector<Eigen::MatrixXi>& Tk_coeffs) {
+                          std::vector<Eigen::MatrixXi>& Tk_coeffs, double damping_factor,
+                          double alpha_i, double alpha_j) {
   int x, y, z;
   Eigen::VectorXd Tk(multipole_components(k));
   int idx = 0;
   for (x = k; x > -1; x--) {
     for (y = k; y > -1; y--) {
       for (z = k; z > -1; z--) {
         if (x + y + z != k) continue;
-        Tk[idx] = T(Rij, x, y, z, Tk_coeffs);
+        Tk[idx] = T(Rij, x, y, z, Tk_coeffs, damping_factor, alpha_i, alpha_j);
         idx++;
       }
     }
@@ -39,7 +39,9 @@ Eigen::VectorXd Tk_tensor(int k, const Eigen::Vector3d& Rij,
 // only 1st derivative supported
 Eigen::VectorXd multipole_derivative(int k, int l, const Eigen::Vector3d& Rji,
                                      Eigen::VectorXd Mkj,
-                                     std::vector<Eigen::MatrixXi>& Tk_coeffs) {
+                                     std::vector<Eigen::MatrixXi>& Tk_coeffs,
+                                     double damping_factor, double alpha_i,
+                                     double alpha_j) {
   if (l > 1) throw std::runtime_error("Only 1st derivatives supported for multipoles");
   Eigen::VectorXd Fi = Eigen::VectorXd::Zero(3);
 
@@ -54,7 +56,7 @@ Eigen::VectorXd multipole_derivative(int k, int l, const Eigen::Vector3d& Rji,
   double symfac;
   // std::cout << "mul k = " << k << std::endl;
   // std::cout << "l = " << l << std::endl;
-  Eigen::VectorXd Tk = Tk_tensor(k + l, Rji, Tk_coeffs);
+  Eigen::VectorXd Tk = Tk_tensor(k + l, Rji, Tk_coeffs, damping_factor, alpha_i, alpha_j);
   for (x = k + l; x > -1; x--) {
     for (y = k + l; y > -1; y--) {
       for (z = k + l; z > -1; z--) {
@@ -97,16 +99,57 @@ int xyz2idx(int x, int y, int z) {
 }
 
 double T(const Eigen::Vector3d& Rij, int x, int y, int z,
-         std::vector<Eigen::MatrixXi>& Cijn) {
+         std::vector<Eigen::MatrixXi>& Cijn, double damping_factor, double alpha_i,
+         double alpha_j) {
   double t = 0.0;
   double R = Rij.norm();
   double Cx, Cy, Cz;
+  int k = x + y + z;
+
+  // Thole Damping
+  std::vector<double> scr_facs;
+  if (damping_factor != 0.0) {
+    if (alpha_i <= 0.0 || alpha_j <= 0.0) {
+      throw std::runtime_error(
+            "Thole damping only valid for non-zero"
+            " isotropic polarizabilities.");
+    }
+    if (k > 3) {
+      throw std::runtime_error("Thole damping only implemented up to third order.");
+    }
+    // Molecular Simulation, 32:6, 471-484, DOI: 10.1080/08927020600631270
+    // v = factor * u, with u = R / (alpha_i * alpha_j)**(1/6)
+    double v = damping_factor * R / std::pow(alpha_i * alpha_j, 1.0 / 6.0);
+    scr_facs = thole_screening_factors(v, k);
+  }
+
   for (size_t l = 0; l <= x; l++) {
     Cx = Cijn[0](x, l) * pow((Rij(0) / R), l);
     for (size_t m = 0; m <= y; m++) {
       Cy = Cx * Cijn[l + x](y, m) * pow((Rij(1) / R), m);
       for (size_t n = 0; n <= z; n++) {
-        Cz = Cy * Cijn[l + x + m + y](z, n) * pow((Rij(2) / R), n);
+        Cz     = Cy * Cijn[l + x + m + y](z, n) * pow((Rij(2) / R), n);
+        int kk = l + m + n;
+        // Thole damping
+        if (scr_facs.size()) {
+          if (kk == 0) {
+            if (k == 0) {
+              Cz *= scr_facs[0];
+            } else if (k == 2) {
+              Cz *= scr_facs[1];
+            }
+          } else if (kk == 1) {
+            if (k == 1) {
+              Cz *= scr_facs[1];
+            } else if (k == 3) {
+              Cz *= scr_facs[2];
+            }
+          } else if (kk == 2) {
+            if (k == 2) Cz *= scr_facs[2];
+          } else if (kk == 3) {
+            if (k == 3) Cz *= scr_facs[3];
+          }
+        }
         t += Cz;
       }
     }
@@ -115,6 +158,29 @@ double T(const Eigen::Vector3d& Rij, int x, int y, int z,
   return t;
 }
 
+std::vector<double> thole_screening_factors(double v, int k) {
+  std::vector<double> ret;
+  // screening factors for potential, field, field gradient, field Hessian
+  double f_v, f_E, f_T, f_D = 1.0;
+  if (k >= 0) {
+    f_v = 1.0 - (0.5 * v + 1.0) * std::exp(-v);
+    ret.push_back(f_v);
+  }
+  if (k >= 1) {
+    f_E = f_v - (0.5 * std::pow(v, 2) + 0.5 * v) * std::exp(-v);
+    ret.push_back(f_E);
+  }
+  if (k >= 2) {
+    f_T = f_E - 1.0 / 6.0 * std::pow(v, 3) * std::exp(-v);
+    ret.push_back(f_T);
+  }
+  if (k >= 3) {
+    f_D = f_T - 1.0 / 30.0 * std::pow(v, 4) * std::exp(-v);
+    ret.push_back(f_D);
+  }
+  return ret;
+}
+
 std::vector<Eigen::MatrixXi> Tk_coefficients(int max_order) {
   int maxi = 2 * max_order + 3;
   std::vector<Eigen::MatrixXi> Cijn;
@@ -219,4 +285,4 @@ std::vector<double> prefactors_nuclei(unsigned k) {
 
 int multipole_components(int k) { return (k + 1) * (k + 2) / 2; }
 
-}  // namespace libcppe
+}  // namespace libcppe

cppe/core/math.hh

@@ -10,13 +10,20 @@ Eigen::Vector3d smat_vec(const Eigen::VectorXd& mat, const Eigen::Vector3d& vec,
 
 Eigen::VectorXd multipole_derivative(int k, int l, const Eigen::Vector3d& Rji,
                                      Eigen::VectorXd Mkj,
-                                     std::vector<Eigen::MatrixXi>& Tk_coeffs);
+                                     std::vector<Eigen::MatrixXi>& Tk_coeffs,
+                                     double damping_factor = 0.0, double alpha_i = 0.0,
+                                     double alpha_j = 0.0);
 
 double T(const Eigen::Vector3d& Rij, int x, int y, int z,
-         std::vector<Eigen::MatrixXi>& Cijn);
+         std::vector<Eigen::MatrixXi>& Cijn, double damping_factor = 0.0,
+         double alpha_i = 0.0, double alpha_j = 0.0);
 
 Eigen::VectorXd Tk_tensor(int k, const Eigen::Vector3d& Rij,
-                          std::vector<Eigen::MatrixXi>& Tk_coeffs);
+                          std::vector<Eigen::MatrixXi>& Tk_coeffs,
+                          double damping_factor = 0.0, double alpha_i = 0.0,
+                          double alpha_j = 0.0);
+
+std::vector<double> thole_screening_factors(double v, int k);
 
 std::vector<Eigen::MatrixXi> Tk_coefficients(int max_order);
 
@@ -37,4 +44,4 @@ std::vector<double> prefactors_nuclei(unsigned k);
 
 int multipole_components(int k);
 
-}  // namespace libcppe
+}  // namespace libcppe

cppe/core/pe_options.hh

@@ -39,6 +39,13 @@ struct PeOptions {
   double diis_start_norm =
         1.0;  ///< maximal residual norm for which DIIS is being enabled
 
+  bool damp_induced             = false;   ///< Enable Thole damping for induced moments
+  double damping_factor_induced = 2.1304;  ///< damping factor for induced moments
+  bool damp_multipole =
+        false;  ///< Enable Thole damping for electric fields created by multipole moments
+  double damping_factor_multipole = 2.1304;  ///< damping factor for electric
+                                             ///< fields created by multipole moments
+
   bool pe_border = false;          ///< Activate border options for sites in proximity
                                    ///< to the QM/MM border
   BorderOptions border_options{};  ///< Options for QM/MM border

cppe/core/potential.hh

@@ -59,6 +59,8 @@ class Polarizability {
     m_values                                      = new_values;
   }
 
+  double get_isotropic_value() { return (m_values[0] + m_values[3] + m_values[5]) / 3.0; }
+
   Eigen::VectorXd get_values_vec() {
     return Eigen::Map<Eigen::VectorXd>(m_values.data(), m_values.size());
   }

cppe/python_iface/export_fields.cc

@@ -16,7 +16,7 @@ void export_fields(py::module& m) {
 
   py::class_<libcppe::MultipoleFields> mul_fields(
         m, "MultipoleFields", "Electric fields created by multipoles");
-  mul_fields.def(py::init<std::vector<libcppe::Potential>>())
+  mul_fields.def(py::init<std::vector<libcppe::Potential>, libcppe::PeOptions>())
         .def("compute", &libcppe::MultipoleFields::compute);
 
   py::class_<libcppe::InducedMoments> ind_moments(m, "InducedMoments");
@@ -28,4 +28,4 @@ void export_fields(py::module& m) {
              py::arg("total_fields"), py::arg("make_guess"));
 
   m.def("get_polarizable_sites", &libcppe::get_polarizable_sites);
-}
+}

cppe/python_iface/export_math.cc

@@ -11,8 +11,12 @@ using namespace libcppe;
 void export_math(py::module& m) {
   m.def("smat_vec", &smat_vec);
   m.def("multipole_derivative", &multipole_derivative);
-  m.def("T", &T);
-  m.def("Tk_tensor", &Tk_tensor);
+  m.def("T", &T, py::arg("Rij"), py::arg("x"), py::arg("y"), py::arg("z"),
+        py::arg("Cijn"), py::arg("damping_factor") = 0.0, py::arg("alpha_i") = 0.0,
+        py::arg("alpha_j") = 0.0);
+  m.def("Tk_tensor", &Tk_tensor, py::arg("k"), py::arg("Rij"), py::arg("Tk_coeffs"),
+        py::arg("damping_factor") = 0.0, py::arg("alpha_i") = 0.0,
+        py::arg("alpha_j") = 0.0);
   m.def("Tk_coefficients", &Tk_coefficients);
   m.def("xyz2idx", &xyz2idx);
   m.def("factorial", &factorial);
@@ -23,4 +27,4 @@ void export_math(py::module& m) {
         "Prefactors for the multipole components (alphabetical order)");
   m.def("prefactors_nuclei", &prefactors_nuclei);
   m.def("multipole_components", &multipole_components);
-}
+}

cppe/python_iface/export_options.cc

@@ -29,7 +29,13 @@ void export_options(py::module& m) {
         .def_readwrite("do_diis", &libcppe::PeOptions::do_diis)
         .def_readwrite("diis_start_norm", &libcppe::PeOptions::diis_start_norm)
         .def_readwrite("maxiter", &libcppe::PeOptions::maxiter)
+        .def_readwrite("damp_induced", &libcppe::PeOptions::damp_induced)
+        .def_readwrite("damping_factor_induced",
+                       &libcppe::PeOptions::damping_factor_induced)
+        .def_readwrite("damp_multipole", &libcppe::PeOptions::damp_multipole)
+        .def_readwrite("damping_factor_multipole",
+                       &libcppe::PeOptions::damping_factor_multipole)
 
         .def_readwrite("pe_border", &libcppe::PeOptions::pe_border)
         .def_readwrite("border_options", &libcppe::PeOptions::border_options);
-}
+}