zhu_philpott.cpp

Go to the documentation of this file.

/*
** This file is part of eOn.
**
** SPDX-License-Identifier: BSD-3-Clause
**
** Copyright (c) 2010--present, eOn Development Team
** All rights reserved.
**
** Repo:
** https://github.com/TheochemUI/eOn
*/
#include "eon/potentials/Water_Pt/zhu_philpott.hpp"
#include <cassert>
#include <cmath>

 
// #define DEBUG_LEVEL 5
#ifdef DEBUG_LEVEL
#warning zhu_philpott debug: DEBUG_LEVEL defined
#define DEBUG_LEVEL_RETURN(x)                                                  \
  if (x == DEBUG_LEVEL)                                                        \
    continue;
#else
#define DEBUG_LEVEL_RETURN(x)
#endif
 
using namespace forcefields;
 
template <class P>
ZhuPhilpott<P>::ZhuPhilpott()
    : SpceCcl() {
  nPlatinum_ = 0;
  positions_ = 0;
  forces_ = 0;
}
 
template <class P>
ZhuPhilpott<P>::ZhuPhilpott(double cutoff, double switchingWidth)
    : SpceCcl(cutoff, switchingWidth) {
  nPlatinum_ = 0;
  positions_ = 0;
  forces_ = 0;
}
 
template <class P> ZhuPhilpott<P>::ZhuPhilpott(ZhuPhilpott const &zhuPhilpott) {
  operator=(zhuPhilpott);
}
 
template <class P>
void ZhuPhilpott<P>::operator=(ZhuPhilpott const &zhuPhilpott) {
  setPlatinum(zhuPhilpott.nPlatinum_, zhuPhilpott.positions_.data());
}
 
template <class P> ZhuPhilpott<P>::~ZhuPhilpott() = default;

template <class P>
void ZhuPhilpott<P>::computeHH_O_Pt_(const int nWater, const int nPt,
                                     const double r[], double f[],
                                     double &energy, double const b[],
                                     bool const fixed[]) {
  const double (*const rh1)[6] = reinterpret_cast<const double (*)[6]>(r);
  const double (*const rh2)[6] = reinterpret_cast<const double (*)[6]>(&r[3]);
  const double (*const ro)[3] =
      reinterpret_cast<const double (*)[3]>(&r[nWater * 6]);
  const double (*const rPt)[3] =
      reinterpret_cast<const double (*)[3]>(&r[nWater * 9]);
 
  double (*const fh1)[6] = reinterpret_cast<double (*)[6]>(f);
  double (*const fh2)[6] = reinterpret_cast<double (*)[6]>(&f[3]);
  double (*const fo)[3] = reinterpret_cast<double (*)[3]>(&f[nWater * 6]);
  double (*const fPt)[3] = reinterpret_cast<double (*)[3]>(&f[nWater * 9]);
 
  if (fixed) {
    bool const(*const xh1)[2] = reinterpret_cast<bool const(*)[2]>(fixed);
    bool const(*const xh2)[2] = reinterpret_cast<bool const(*)[2]>(&fixed[1]);
    bool const(*const xo)[1] =
        reinterpret_cast<bool const(*)[1]>(&fixed[nWater * 2]);
    bool const *const xPt = &fixed[nWater * 3];
    computeTemplate(nWater, rh1, rh2, ro, fh1, fh2, fo, nPt, rPt, fPt, energy,
                    b, xh1, xh2, xo, xPt);
  } else {
    bool const(*const xh1)[2] = 0;
    bool const(*const xh2)[2] = 0;
    bool const(*const xo)[1] = 0;
    computeTemplate(nWater, rh1, rh2, ro, fh1, fh2, fo, nPt, rPt, fPt, energy,
                    b, xh1, xh2, xo, 0);
  };
}

template <class P>
void ZhuPhilpott<P>::computeHH_O_(const int nWater, const double r[],
                                  double f[], double &energy, double const b[],
                                  bool const fixed[]) {
  const double (*const rh1)[6] = reinterpret_cast<const double (*)[6]>(r);
  const double (*const rh2)[6] = reinterpret_cast<const double (*)[6]>(&r[3]);
  const double (*const ro)[3] =
      reinterpret_cast<const double (*)[3]>(&r[nWater * 6]);
  const double (*const rPt)[3] =
      reinterpret_cast<const double (*)[3]>(positions_.data());
 
  double (*const fh1)[6] = reinterpret_cast<double (*)[6]>(f);
  double (*const fh2)[6] = reinterpret_cast<double (*)[6]>(&f[3]);
  double (*const fo)[3] = reinterpret_cast<double (*)[3]>(&f[nWater * 6]);
  double (*const fPt)[3] = reinterpret_cast<double (*)[3]>(forces_.data());
 
  if (fixed) {
    bool const(*const xh1)[2] = reinterpret_cast<bool const(*)[2]>(fixed);
    bool const(*const xh2)[2] = reinterpret_cast<bool const(*)[2]>(&fixed[1]);
    bool const(*const xo)[1] =
        reinterpret_cast<bool const(*)[1]>(&fixed[nWater * 2]);
    computeTemplate(nWater, rh1, rh2, ro, fh1, fh2, fo, nPlatinum_, rPt, fPt,
                    energy, b, xh1, xh2, xo, 0);
  } else {
    bool const(*const xh1)[2] = 0;
    bool const(*const xh2)[2] = 0;
    bool const(*const xo)[1] = 0;
    computeTemplate(nWater, rh1, rh2, ro, fh1, fh2, fo, nPlatinum_, rPt, fPt,
                    energy, b, xh1, xh2, xo, 0);
  };
}

template <class P> int ZhuPhilpott<P>::nPlatinum() const { return nPlatinum_; }

template <class P>
void ZhuPhilpott<P>::setPlatinum(int nPlatinum, double const positions[]) {
  assert(nPlatinum >= 0);
  nPlatinum_ = nPlatinum;
  int const n = nPlatinum * 3;
  positions_.assign(positions, positions + n);
  forces_.assign(n, 0.0);
}

template <class P> char const *ZhuPhilpott<P>::getName() {
  return "ZhuPhilpott";
}
 
template <class P>
template <int H, int O, int H3, int O3>
void ZhuPhilpott<P>::computeTemplate(
    const int nWater, const double (*const rh1)[H3],
    const double (*const rh2)[H3], const double (*const ro)[O3],
    double (*const fh1)[H3], double (*const fh2)[H3], double (*const fo)[O3],
    const int nPt, const double rPt[][3], double fPt[][3], double &energy,
    double const b[], bool const (*const xh1)[H], bool const (*const xh2)[H],
    bool const (*const xo)[O], bool const xPt[]) {
  for (int i = 0; i < nWater; ++i) {
    for (int a = 0; a < 3; ++a) {
      fh1[i][a] = 0.0;
      fh2[i][a] = 0.0;
      fo[i][a] = 0.0;
    };
  };
  for (int i = 0; i < nPt; ++i) {
    for (int a = 0; a < 3; ++a)
      fPt[i][a] = 0.0;
  };
  energy = 0.0;
  setPeriodicity(b);
  for (int i = nWater - 1; i >= 0; --i) {
    double centre1[3] = {0};
    calculateCentre(rh1[i], rh2[i], ro[i], centre1);
    Water w1(rh1[i], rh2[i], ro[i], centre1, fh1[i], fh2[i], fo[i]);
    intramolecular(w1, energy);
    assert(not std::isnan(energy));
    DEBUG_LEVEL_RETURN(0)
    // Two next lines for interaction with platinum
    interactWithCorePt(w1, nPt, rPt, fPt,
                       energy); // called Vw-core in Zhu and Philpott
    assert(not std::isnan(energy));
    DEBUG_LEVEL_RETURN(1)
    interactWithImage(w1, w1, energy); // called Vw-cond in Zhu and Philpott
    DEBUG_LEVEL_RETURN(2)
    assert(not std::isnan(energy));
    for (int j = i - 1; j >= 0; --j) {
      bool areFixed = false;
      if (xh1 and xh2 and xo) {
        areFixed = xh1[i][0] and xh2[i][0] and xo[i][0];
        // check if all the atoms of molecule j are fixed.
        areFixed &= xh1[j][0] and xh2[j][0] and xo[j][0];
      };
      // if both molecules are fixed skip force calculation
      if (not areFixed) {
        double centre2[3] = {0};
        calculateCentre(rh1[j], rh2[j], ro[j], centre2);
        Water w2(rh1[j], rh2[j], ro[j], centre2, fh1[j], fh2[j], fo[j]);
        lennardJonesWithCutoff(w1, w2, energy);
        DEBUG_LEVEL_RETURN(3)
        coulombWithCutoff(w1, w2, energy, 1.0);
        DEBUG_LEVEL_RETURN(4)
        // Two next lines for interaction with platinum
        // called Vw-cond in Zhu and Philpott
        interactWithImage(w1, w2, energy); // w1 interacts with image of w2
        interactWithImage(w2, w1, energy); // w2 interacts with image of w1
        DEBUG_LEVEL_RETURN(5)
      };
    };
#ifndef NDEBUG
    for (int a = 0; a < 3; a++) {
      assert(not std::isnan(fh1[i][a]));
      assert(not std::isnan(fh2[i][a]));
      assert(not std::isnan(fo[i][a]));
    };
#endif
  };
  assert(not std::isnan(energy) and not std::isinf(energy));
}

template <class P>
void ZhuPhilpott<P>::interactWithCorePt(Water &water, int const nPt,
                                        double const rPt[][3], double fPt[][3],
                                        double &energy) {
  for (int i = 0; i < nPt; ++i) {
    double r[3], r1; // for distance between two atoms
                     // Pt - O
    distance(water.ro_, rPt[i], r, r1);
    if (r1 <= cutoff_ - switchingWidth_) {
      interactionPtO(water.ro_, rPt[i], water.fo_, fPt[i], energy);
    } else if (r1 < cutoff_) {
      double en = 0;
      double fo[3] = {0}, fPtTmp[3] = {0};
      interactionPtO(water.ro_, rPt[i], fo, fPtTmp, en);
      switching(water.ro_, rPt[i], fo, fPtTmp, en);
      energy += en;
      for (int j = 0; j < 3; j++) {
        water.fo_[j] += fo[j];
        fPt[i][j] += fPtTmp[j];
      };
    };
    // Pt - H1
    distance(water.rh1_, rPt[i], r, r1);
    if (r1 <= cutoff_ - switchingWidth_) {
      interactionPtH(water.rh1_, rPt[i], water.fh1_, fPt[i], energy);
    } else if (r1 < cutoff_) {
      double en = 0;
      double fh[3] = {0}, fPtTmp[3] = {0};
      interactionPtH(water.rh1_, rPt[i], fh, fPtTmp, en);
      switching(water.rh1_, rPt[i], fh, fPtTmp, en);
      energy += en;
      for (int j = 0; j < 3; j++) {
        water.fh1_[j] += fh[j];
        fPt[i][j] += fPtTmp[j];
      };
    };
    // Pt - H2
    distance(water.rh2_, rPt[i], r, r1);
    if (r1 <= cutoff_ - switchingWidth_) {
      interactionPtH(water.rh2_, rPt[i], water.fh2_, fPt[i], energy);
    } else if (r1 < cutoff_) {
      double en = 0;
      double fh[3] = {0}, fPtTmp[3] = {0};
      interactionPtH(water.rh2_, rPt[i], fh, fPtTmp, en);
      switching(water.rh2_, rPt[i], fh, fPtTmp, en);
      energy += en;
      for (int j = 0; j < 3; j++) {
        water.fh2_[j] += fh[j];
        fPt[i][j] += fPtTmp[j];
      };
    };
  };
}
 
template <class P>
void ZhuPhilpott<P>::interactionPtO(double const R1[], double const R2[],
                                    double F1[], double F2[], double &energy) {
  double R12[3], F12[3] = {0}, r1, f1 = 0;
  distance(R1, R2, R12, r1);
  anisotropic(R12, F12, energy, P::epsilonOPt_, P::sigmaOPt_, P::alpha_);
  isotropic10(r1, f1, energy, P::epsilonOPt_, P::sigmaOPt_, P::C10_O_);
  for (int j = 0; j < 3; j++) {
    F1[j] += F12[j] + f1 * R12[j] / r1;
    F2[j] -= F12[j] + f1 * R12[j] / r1;
  };
}
 
template <class P>
void ZhuPhilpott<P>::interactionPtH(double const R1[], double const R2[],
                                    double F1[], double F2[], double &energy) {
  double R12[3], F12[3] = {0}, r1, f1 = 0;
  distance(R1, R2, R12, r1);
  anisotropic(R12, F12, energy, P::epsilonHPt_, P::sigmaHPt_, P::alpha_);
  isotropic10(r1, f1, energy, P::epsilonHPt_, P::sigmaHPt_, P::C10_H_);
  for (int j = 0; j < 3; j++) {
    F1[j] += F12[j] + f1 * R12[j] / r1;
    F2[j] -= F12[j] + f1 * R12[j] / r1;
  };
}

template <class P>
void ZhuPhilpott<P>::anisotropic(const double distance[], double force[],
                                 double &energy, double const epsilon,
                                 double const sigma, double const alpha) {
  double z2, z, a, b, A, B, dE_da, dE_db, rho, rho2, alpha2;
  // WARNING: F1 and F2 are incremented.
  alpha2 = alpha * alpha;
  double f;
  rho2 = distance[0] * distance[0] + distance[1] * distance[1];
  rho = std::sqrt(rho2);
  z = distance[2];
  z2 = z * z;
  double sigma2 = sigma * sigma;
  a = sigma2 / (rho2 * alpha2 + z2);
  b = sigma2 / (rho2 / alpha2 + z2);
  double a2 = a * a, a3 = a2 * a;
  double b2 = b * b;
  A = a3 * a3; // a^6
  B = b * b2;  // b^3
  energy += 4.0 * epsilon * (A - B);
  dE_da = 24.0 * epsilon * a2 * a3; // a^5
  dE_db = -12.0 * epsilon * b2;     // b^2
  double invSigma2 = 1.0 / sigma2;
  f = 2.0 * (dE_da * a2 * alpha2 + dE_db * b2 / alpha2) * invSigma2;
  force[0] += f * distance[0];
  force[1] += f * distance[1];
  f = 2.0 * (dE_da * a2 + dE_db * b2) * invSigma2;
  force[2] += f * distance[2];
}

template <class P>
void ZhuPhilpott<P>::isotropic10(double const distance, double &force,
                                 double &energy, double const epsilon,
                                 double const sigma, double const C10) {
  double const s = sigma / distance;
  double s10 = s * s * s * s * s;
  s10 *= s10;
  energy += -4.0 * epsilon * C10 * s10;
  force -= 40.0 * epsilon * C10 * s10 / distance;
}

template <class P>
void ZhuPhilpott<P>::interactWithImage(Water &w1, Water &w2, double &U) {
  // f is used to store the force on the atoms' images. These forces have no use
  // and are at the end discarded. we use the same vector to store of the forces
  // even when they are on different atoms' images.
  double r[4][3], f[3][3] = {{0}};
  for (int k = 0; k < 2; k++) {
    r[0][k] = w2.rh1_[k];
    r[1][k] = w2.rh2_[k];
    r[2][k] = w2.ro_[k];
    r[3][k] = w2.rc_[k];
  };
  // mirror about plane at z=0.
  r[0][2] = -w2.rh1_[2];
  r[1][2] = -w2.rh2_[2];
  r[2][2] = -w2.ro_[2];
  r[3][2] = -w2.rc_[2];
  Water w2m(r[0], r[1], r[2], r[3], f[0], f[1], f[2]);
  // -2 is given as the relative permittivity. This value is not the
  // permittivity of the metal which is infinite in the case of a perfect metal.
  // It is the value (1+epsilon_r)/(1-epsilon_r) which is equal to -1 for a
  // perfect metal. This value is the 'effective' permittivity for the
  // interaction of a charge with a mirror image.
  coulombWithCutoff(w1, w2m, U, -2.0);
  for (int k = 0; k < 2; k++) {
    w2.fh1_[k] += f[0][k];
    w2.fh2_[k] += f[1][k];
    w2.fo_[k] += f[2][k];
  };
  // mirror about plane
  w2.fh1_[2] -= f[0][2];
  w2.fh2_[2] -= f[1][2];
  w2.fo_[2] -= f[2][2];
}

template <class P>
void ZhuPhilpott<P>::coulombWithCutoff(Water &w1, Water &w2, double &U,
                                       double const relativePermittivity) {
  double z[3], z1, z2;
  distance(w1.rc_, w2.rc_, z, z1, z2);
  if (z1 <= cutoff_ - switchingWidth_) {
    coulombFull(w1, w2, U, relativePermittivity);
  } else if (z1 < cutoff_) {
    double f1[3][3] = {{0}}, f2[3][3] = {{0}};
    Water v1(w1, f1[0], f1[1], f1[2]);
    Water v2(w2, f2[0], f2[1], f2[2]);
    double energy = 0.0;
    coulombFull(v1, v2, energy, relativePermittivity);
    ChargeGroup<3> g1 = {v1.rc_, 0, f1};
    ChargeGroup<3> g2 = {v2.rc_, 0, f2};
    switching(g1, g2, energy, cutoff_, switchingWidth_);
    U += energy;
    w1.addForces(v1);
    w2.addForces(v2);
  };
}

template <class P>
void ZhuPhilpott<P>::coulombFull(Water &w1, Water &w2, double &U,
                                 double const relativePermittivity) {
  // q^2/relativePermittivity;
  const double qq2overEr = charge2_ / relativePermittivity;
  // Coulomb interactions between hydrogens
  coulomb(w1.rh1_, w2.rh1_, w1.fh1_, w2.fh1_, U, qq2overEr);
  coulomb(w1.rh1_, w2.rh2_, w1.fh1_, w2.fh2_, U, qq2overEr);
  coulomb(w1.rh2_, w2.rh1_, w1.fh2_, w2.fh1_, U, qq2overEr);
  coulomb(w1.rh2_, w2.rh2_, w1.fh2_, w2.fh2_, U, qq2overEr);
  // interactions between H and O.
  coulomb(w1.ro_, w2.rh1_, w1.fo_, w2.fh1_, U, -2.0 * qq2overEr);
  coulomb(w1.ro_, w2.rh2_, w1.fo_, w2.fh2_, U, -2.0 * qq2overEr);
  coulomb(w1.rh1_, w2.ro_, w1.fh1_, w2.fo_, U, -2.0 * qq2overEr);
  coulomb(w1.rh2_, w2.ro_, w1.fh2_, w2.fo_, U, -2.0 * qq2overEr);
  // interactions between O1, O2
  coulomb(w1.ro_, w2.ro_, w1.fo_, w2.fo_, U, 4.0 * qq2overEr);
}