MM part for nuclear gradients

cppe/core/cppe_state.cc

@@ -57,7 +57,7 @@ void CppeState::set_potentials(std::vector<Potential> potentials) {
 void CppeState::calculate_static_energies_and_fields() {
  // Electrostatic energy (nuclei-multipoles)
  MultipoleExpansion mexp(m_mol, m_potentials);
- double nuc_mul_energy = mexp.calculate_interaction_energy();
+ double nuc_mul_energy = mexp.interaction_energy();
  m_pe_energy["Electrostatic"]["Nuclear"] = nuc_mul_energy;
 
  // Calculate static fields
@@ -70,6 +70,16 @@ void CppeState::calculate_static_energies_and_fields() {
  m_multipole_fields = mul_fields.compute();
 }
 
+Eigen::MatrixXd CppeState::nuclear_interaction_energy_gradient() {
+ MultipoleExpansion mexp(m_mol, m_potentials);
+ return mexp.nuclear_gradient();
+}
+
+Eigen::MatrixXd CppeState::nuclear_field_gradient() {
+ NuclearFields nfields(m_mol, m_potentials);
+ return nfields.nuclear_gradient();
+}
+
 void CppeState::update_induced_moments(Eigen::VectorXd elec_fields, bool elec_only) {
  Eigen::VectorXd tmp_total_fields = Eigen::VectorXd::Zero(m_polarizable_sites * 3);
  if (elec_only) {

cppe/core/cppe_state.hh

@@ -85,6 +85,9 @@ class CppeState {
  std::string get_energy_summary_string();
 
  Eigen::MatrixXd induced_moments_eef();
+
+ Eigen::MatrixXd nuclear_interaction_energy_gradient();
+ Eigen::MatrixXd nuclear_field_gradient();
 };
 
 } // namespace libcppe

cppe/core/electric_fields.cc

@@ -74,6 +74,35 @@ Eigen::VectorXd NuclearFields::compute() {
  return nuc_fields;
 }
 
+Eigen::MatrixXd NuclearFields::nuclear_gradient() {
+ int natoms = m_mol.size();
+ Eigen::MatrixXd grad = Eigen::MatrixXd::Zero(3 * natoms, 3 * m_n_polsites);
+#pragma omp parallel for
+ for (size_t i = 0; i < m_n_polsites; i++) {
+ size_t site_counter = 3 * i;
+ Potential& potential = m_polsites[i];
+ Eigen::Vector3d site_position = potential.get_site_position();
+ for (int ai = 0; ai < natoms; ++ai) {
+ auto& atom = m_mol[ai];
+ Eigen::Vector3d core_position = atom.get_position();
+ Eigen::Vector3d diff = site_position - core_position;
+ Eigen::VectorXd Tms = tensors::T2(diff);
+ grad(3 * ai + 0, site_counter + 0) = atom.charge * Tms(0);
+ grad(3 * ai + 0, site_counter + 1) = atom.charge * Tms(1);
+ grad(3 * ai + 0, site_counter + 2) = atom.charge * Tms(2);
+
+ grad(3 * ai + 1, site_counter + 0) = atom.charge * Tms(1);
+ grad(3 * ai + 1, site_counter + 1) = atom.charge * Tms(3);
+ grad(3 * ai + 1, site_counter + 2) = atom.charge * Tms(4);
+
+ grad(3 * ai + 2, site_counter + 0) = atom.charge * Tms(2);
+ grad(3 * ai + 2, site_counter + 1) = atom.charge * Tms(4);
+ grad(3 * ai + 2, site_counter + 2) = atom.charge * Tms(5);
+ }
+ }
+ return grad;
+}
+
 Eigen::VectorXd MultipoleFields::compute_tree() {
  int n_sites = m_potentials.size();
  int n_crit = m_options.tree_ncrit;

cppe/core/electric_fields.hh

@@ -32,6 +32,7 @@ class NuclearFields : public ElectricFields {
  NuclearFields(Molecule mol, std::vector<Potential> potentials)
  : ElectricFields(potentials), m_mol(mol){};
  Eigen::VectorXd compute();
+ Eigen::MatrixXd nuclear_gradient();
 };
 
 class MultipoleFields : public ElectricFields {

cppe/core/multipole_expansion.cc

@@ -1,12 +1,13 @@
 #include <Eigen/Dense>
 
 #include "multipole_expansion.hh"
+#include "electric_fields.hh"
 #include "tensors.hh"
 
 namespace libcppe {
 
 // calculates the multipole-nuclei interaction energy through the given order
-double MultipoleExpansion::calculate_interaction_energy() {
+double MultipoleExpansion::interaction_energy() {
  double total_energy = 0.0;
  int npots = m_potentials.size();
 
@@ -31,4 +32,34 @@ double MultipoleExpansion::calculate_interaction_energy() {
  return total_energy;
 }
 
+Eigen::MatrixXd MultipoleExpansion::nuclear_gradient() {
+ int natoms = m_mol.size();
+ int npots = m_potentials.size();
+ Eigen::MatrixXd grad = Eigen::MatrixXd::Zero(natoms, 3);
+
+ #pragma omp parallel for
+ for (size_t i = 0; i < npots; i++) {
+ Potential& potential = m_potentials[i];
+ Eigen::Vector3d site_position = potential.get_site_position();
+ for (auto& multipole : potential.get_multipoles()) {
+ // TODO: refactor in math.cc
+ std::vector<double> pref = prefactors_nuclei(multipole.m_k);
+ Eigen::VectorXd pref_v =
+ Eigen::Map<Eigen::VectorXd>(std::move(pref.data()), pref.size());
+ Eigen::VectorXd mul_v = multipole.get_values_vec();
+ for (int ai = 0; ai < natoms; ++ai) {
+ auto& atom = m_mol[ai];
+ Eigen::Vector3d core_position = atom.get_position();
+ Eigen::Vector3d diff = core_position - site_position;
+ Eigen::Vector3d tmp_grad = multipole_derivative(multipole.m_k, 1, diff, multipole.get_values_vec());
+ #pragma omp critical
+ {
+ grad.row(ai) += -tmp_grad * atom.charge;
+ }
+ }
+ }
+ }
+ return grad;
+}
+
 } // namespace libcppe

cppe/core/multipole_expansion.hh

@@ -16,9 +16,9 @@ class MultipoleExpansion {
  m_potentials(potentials){
 
  };
- ~MultipoleExpansion(){};
 
- double calculate_interaction_energy();
+ double interaction_energy();
+ Eigen::MatrixXd nuclear_gradient();
 };
 
 } // namespace libcppe

cppe/python_iface/export_fields.cc

@@ -3,6 +3,7 @@
 #include <pybind11/stl.h>
 
 #include "../core/electric_fields.hh"
+#include "../core/multipole_expansion.hh"
 #include "../core/molecule.hh"
 #include "../core/pe_options.hh"
 
@@ -34,7 +35,8 @@ void export_fields(py::module& m) {
  py::class_<libcppe::NuclearFields> nuc_fields(m, "NuclearFields",
  "Electric fields created by nuclei");
  nuc_fields.def(py::init<libcppe::Molecule, std::vector<libcppe::Potential>>())
- .def("compute", &libcppe::NuclearFields::compute);
+ .def("compute", &libcppe::NuclearFields::compute)
+ .def("nuclear_gradient", &libcppe::NuclearFields::nuclear_gradient);
 
  py::class_<libcppe::MultipoleFields, std::shared_ptr<libcppe::MultipoleFields>>
  mul_fields(m, "MultipoleFields", "Electric fields created by multipoles");
@@ -71,5 +73,10 @@ void export_fields(py::module& m) {
  .def("apply_diagonal_inverse", &libcppe::BMatrix::apply_diagonal_inverse);
  bmatrix.def_property_readonly("exclusions", &libcppe::BMatrix::get_exclusions);
 
+ py::class_<libcppe::MultipoleExpansion, std::shared_ptr<libcppe::MultipoleExpansion>> multipole_expansion(m, "MultipoleExpansion");
+ multipole_expansion.def(py::init<libcppe::Molecule, std::vector<libcppe::Potential>>())
+ .def("interaction_energy", &libcppe::MultipoleExpansion::interaction_energy)
+ .def("nuclear_gradient", &libcppe::MultipoleExpansion::nuclear_gradient);
+
  m.def("multipole_derivative", &libcppe::multipole_derivative);
 }

cppe/python_iface/export_state.cc

@@ -140,6 +140,10 @@ void export_state(py::module& m) {
  .def("set_potentials", &libcppe::CppeState::set_potentials)
  .def("calculate_static_energies_and_fields",
  &libcppe::CppeState::calculate_static_energies_and_fields)
+ .def("nuclear_interaction_energy_gradient",
+ &libcppe::CppeState::nuclear_interaction_energy_gradient)
+ .def("nuclear_field_gradient",
+ &libcppe::CppeState::nuclear_field_gradient)
  .def("get_induced_moments", &libcppe::CppeState::get_induced_moments)
  .def("induced_moments_eef", &libcppe::CppeState::induced_moments_eef)
  .def_readwrite("energies", &libcppe::CppeState::m_pe_energy)

tests/test_gradients.py

@@ -0,0 +1,99 @@
+import unittest
+import os
+
+import numpy as np
+
+from .cache import cache
+
+import cppe
+from cppe import PotfileReader, MultipoleExpansion, NuclearFields
+from polarizationsolver import tensors
+
+
+prefactors_5p = np.array([1.0, -8.0, 8.0, -1.0]) / 12.0
+multipliers_5p = [-2, -1, 1, 2]
+prefactors_2p = [-0.5, 0.5]
+multipliers_2p = [-1, 1]
+stencils = {
+ "2p": (multipliers_2p, prefactors_2p),
+ "5p": (multipliers_5p, prefactors_5p)
+}
+coords_label = ["x", "y", "z"]
+
+
+def print_callback(output):
+ print("cb: ", output)
+
+
+class TestGradients(unittest.TestCase):
+ dirname = os.path.dirname(__file__)
+ potfile_path = "{}/potfiles/pna_6w.pot".format(dirname)
+ potfile_iso_path = "{}/potfiles/pna_6w_isopol.pot".format(dirname)
+
+ def test_compute_nuclear_interaction_energy_gradient(self):
+ mol = cache.molecule["pna"]
+ options = {"potfile": self.potfile_path}
+
+ p = PotfileReader(self.potfile_path)
+ potentials = p.read()
+
+ atoms = [a.atomic_number for a in mol]
+ coords = np.array([a.position for a in mol])
+
+ grad_nuc = MultipoleExpansion(mol, potentials).nuclear_gradient()
+
+ grad_nuc_fd = grad_nuclear_interaction_energy_fdiff(atoms, coords, potentials)
+ np.testing.assert_allclose(grad_nuc_fd, grad_nuc, atol=1e-10)
+
+ grad_field_fd = grad_nuclear_field_fdiff(atoms, coords, potentials)
+
+ natoms = len(atoms)
+ polsites = cppe.get_polarizable_sites(potentials)
+ grad_field = NuclearFields(mol, potentials).nuclear_gradient().reshape((natoms, 3, len(polsites), 3))
+ np.testing.assert_allclose(grad_field_fd, grad_field, atol=1e-10)
+
+def grad_nuclear_field_fdiff(atoms, coords, potentials, step_au=1e-3):
+ """
+ Computes the finite difference gradient
+ with a 5-point stencil
+ """
+ natoms = len(atoms)
+ polsites = cppe.get_polarizable_sites(potentials)
+ grad = np.zeros((natoms, 3, len(polsites), 3))
+ for i in range(natoms):
+ for c in range(3):
+ print("Computing dE/d{} for atom {}.".format(coords_label[c], i))
+ for f, p in zip(*stencils["5p"]):
+ coords_p = coords.copy()
+ coords_p[i, c] += f * step_au
+ mol_p = cppe.Molecule()
+ for z, coord in zip(atoms, coords_p):
+ mol_p.append(cppe.Atom(z, *coord))
+
+ nf = NuclearFields(mol_p, potentials)
+ en_pert = nf.compute().reshape(len(polsites), 3)
+ grad[i, c] += p * en_pert / step_au
+ return grad
+
+
+def grad_nuclear_interaction_energy_fdiff(atoms, coords, potentials, step_au=1e-3):
+ """
+ Computes the finite difference gradient
+ with a 5-point stencil
+ """
+ natoms = len(atoms)
+ grad = np.zeros((natoms, 3))
+ for i in range(natoms):
+ for c in range(3):
+ print("Computing dE/d{} for atom {}.".format(coords_label[c], i))
+ for f, p in zip(*stencils["5p"]):
+ coords_p = coords.copy()
+ coords_p[i, c] += f * step_au
+ mol_p = cppe.Molecule()
+ for z, coord in zip(atoms, coords_p):
+ mol_p.append(cppe.Atom(z, *coord))
+
+ mexp = MultipoleExpansion(mol_p, potentials)
+ en_pert = mexp.interaction_energy()
+ grad[i, c] += p * en_pert / step_au
+ return grad