Program Listing for File atomlist_potential.hpp

↰ Return to documentation for file (source/pele/atomlist_potential.hpp)

#ifndef _ATOMLIST_POTENTIAL_H_
#define _ATOMLIST_POTENTIAL_H_

#include <iostream>

#include "array.hpp"
#include "pairwise_potential_interface.hpp"

namespace pele {

template <typename pairwise_interaction, typename distance_policy>
class AtomListPotential : public PairwisePotentialInterface {
 protected:
  std::shared_ptr<pairwise_interaction> _interaction;
  std::shared_ptr<distance_policy> _dist;
  Array<size_t> _atoms1;
  Array<size_t> _atoms2;
  bool _one_list;
  static const size_t _ndim = distance_policy::_ndim;

  AtomListPotential(std::shared_ptr<pairwise_interaction> interaction,
                    std::shared_ptr<distance_policy> dist,
                    Array<size_t> &atoms1, Array<size_t> &atoms2,
                    const Array<double> radii)
      : PairwisePotentialInterface(radii),
        _interaction(interaction),
        _dist(dist),
        _atoms1(atoms1.copy()),
        _atoms2(atoms2.copy()),
        _one_list(false) {}

  AtomListPotential(std::shared_ptr<pairwise_interaction> interaction,
                    std::shared_ptr<distance_policy> dist,
                    Array<size_t> &atoms1, Array<size_t> &atoms2)
      : _interaction(interaction),
        _dist(dist),
        _atoms1(atoms1.copy()),
        _atoms2(atoms2.copy()),
        _one_list(false) {}

  AtomListPotential(std::shared_ptr<pairwise_interaction> interaction,
                    std::shared_ptr<distance_policy> dist,
                    Array<size_t> &atoms1)
      : _interaction(interaction),
        _dist(dist),
        _atoms1(atoms1.copy()),
        _atoms2(_atoms1),
        _one_list(true) {}

 public:
  virtual inline double get_energy(Array<double> const &x) {
    double e = 0.;
    size_t jstart = 0;
    double dr[_ndim];

    for (size_t i = 0; i < _atoms1.size(); ++i) {
      size_t atom1 = _atoms1[i];
      size_t i1 = _ndim * atom1;
      if (_one_list) jstart = i + 1;
      for (size_t j = jstart; j < _atoms2.size(); ++j) {
        size_t atom2 = _atoms2[j];
        size_t i2 = _ndim * atom2;

        _dist->get_rij(dr, &x[i1], &x[i2]);
        double r2 = 0;
        for (size_t k = 0; k < _ndim; ++k) {
          r2 += dr[k] * dr[k];
        }
        e += _interaction->energy(r2, get_non_additive_cutoff(atom1, atom2));
      }
    }

    return e;
  }

  virtual inline double add_energy_gradient(Array<double> const &x,
                                            Array<double> &grad) {
    if (x.size() != grad.size()) {
      throw std::invalid_argument("the gradient has the wrong size");
    }
    double e = 0.;
    double gij;
    size_t jstart = 0;
    double dr[_ndim];
    ;

    for (size_t i = 0; i < _atoms1.size(); ++i) {
      size_t atom1 = _atoms1[i];
      size_t i1 = _ndim * atom1;
      if (_one_list) {
        jstart = i + 1;
      }
      for (size_t j = jstart; j < _atoms2.size(); ++j) {
        size_t atom2 = _atoms2[j];
        size_t i2 = _ndim * atom2;

        _dist->get_rij(dr, &x[i1], &x[i2]);

        double r2 = 0;
        for (size_t k = 0; k < _ndim; ++k) {
          r2 += dr[k] * dr[k];
        }

        e += _interaction->energy_gradient(
            r2, &gij, get_non_additive_cutoff(atom1, atom2));
        for (size_t k = 0; k < _ndim; ++k) grad[i1 + k] -= gij * dr[k];
        for (size_t k = 0; k < _ndim; ++k) grad[i2 + k] += gij * dr[k];
      }
    }

    return e;
  }

  virtual inline double add_energy_gradient_hessian(Array<double> const &x,
                                                    Array<double> &grad,
                                                    Array<double> &hess) {
    if (x.size() != grad.size()) {
      throw std::invalid_argument("the gradient has the wrong size");
    }
    if (hess.size() != x.size() * x.size()) {
      throw std::invalid_argument("the Hessian has the wrong size");
    }

    double e = 0.;
    double hij, gij;
    size_t jstart = 0;
    double dr[_ndim];
    const size_t N = x.size();

    for (size_t i = 0; i < _atoms1.size(); ++i) {
      size_t atom1 = _atoms1[i];
      size_t i1 = _ndim * atom1;

      if (_one_list) {
        jstart = i + 1;
      }
      for (size_t j = jstart; j < _atoms2.size(); ++j) {
        size_t atom2 = _atoms2[j];
        size_t i2 = _ndim * atom2;

        _dist->get_rij(dr, &x[i1], &x[i2]);
        double r2 = 0;
        for (size_t k = 0; k < _ndim; ++k) {
          r2 += dr[k] * dr[k];
        }

        e += _interaction->energy_gradient_hessian(
            r2, &gij, &hij, get_non_additive_cutoff(atom1, atom2));
        for (size_t k = 0; k < _ndim; ++k) grad[i1 + k] -= gij * dr[k];
        for (size_t k = 0; k < _ndim; ++k) grad[i2 + k] += gij * dr[k];

        for (size_t k = 0; k < _ndim; ++k) {
          // diagonal block - diagonal terms
          double Hii_diag = (hij + gij) * dr[k] * dr[k] / r2 - gij;
          hess[N * (i1 + k) + i1 + k] += Hii_diag;
          hess[N * (i2 + k) + i2 + k] += Hii_diag;
          // off diagonal block - diagonal terms
          double Hij_diag = -Hii_diag;
          hess[N * (i1 + k) + i2 + k] = Hij_diag;
          hess[N * (i2 + k) + i1 + k] = Hij_diag;
          for (size_t l = k + 1; l < _ndim; ++l) {
            // diagonal block - off diagonal terms
            double Hii_off = (hij + gij) * dr[k] * dr[l] / r2;
            hess[N * (i1 + k) + i1 + l] += Hii_off;
            hess[N * (i1 + l) + i1 + k] += Hii_off;
            hess[N * (i2 + k) + i2 + l] += Hii_off;
            hess[N * (i2 + l) + i2 + k] += Hii_off;
            // off diagonal block - off diagonal terms
            double Hij_off = -Hii_off;
            hess[N * (i1 + k) + i2 + l] = Hij_off;
            hess[N * (i1 + l) + i2 + k] = Hij_off;
            hess[N * (i2 + k) + i1 + l] = Hij_off;
            hess[N * (i2 + l) + i1 + k] = Hij_off;
          }
        }
      }
    }

    return e;
  }
};
}  // namespace pele

#endif