potential_base.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/potential_base.hpp"
#include <cassert>
#include <cmath>
#include <cstdlib>
#include <iostream>

using namespace forcefields;

PotentialBase::PotentialBase()
    : cutoff_(6.5),
      switchingWidth_(2.0) {
  periods_[0] = 0.0;
  periods_[1] = 0.0;
  periods_[2] = 0.0;
}
 
PotentialBase::PotentialBase(double cutoff, double switchingWidth)
    : cutoff_(cutoff),
      switchingWidth_(switchingWidth) {
  if (switchingWidth_ > cutoff_) {
    std::cerr << "Error: getSwitchingWidth() > getCutoff()" << std::endl;
    std::exit(EXIT_FAILURE);
  };
  periods_[0] = 0.0;
  periods_[1] = 0.0;
  periods_[2] = 0.0;
}
 
double PotentialBase::getCutoff() const { return cutoff_; }
 
void PotentialBase::setCutoff(double cutoff) { cutoff_ = cutoff; }

double PotentialBase::getSwitchingWidth() const { return switchingWidth_; }
 
void PotentialBase::setSwitchingWidth(double width) { switchingWidth_ = width; }
// @}

double PotentialBase::applyPeriodicity0(double r, double const period) {
  if (not std::isinf(period)) {
    double n = r / period + 0.5;
    // This is slightly more efficient than using function floor().
    int m = static_cast<int>(n);
    if (n < 0)
      --m;
    r -= m * period;
  };
  return r;
}

void PotentialBase::applyPeriodicity0(double r[], double const periods[]) {
  for (int i = 0; i < 3; ++i)
    r[i] = applyPeriodicity0(r[i], periods[i]);
}

double PotentialBase::isoscelesBase(double length, double angle) {
  return 2.0 * length * std::sin(angle / 2.0);
}

void PotentialBase::restrainLength(const double R1[], const double R2[],
                                   double F1[], double F2[], double &u,
                                   const double k, const double req) {
  double f, fx;
  double r[3], r1, r2;
  distance(R2, R1, r, r1, r2);
  double d = r1 - req;
  f = -2.0 * k * d;
  u += k * d * d;
  fx = f * r[0] / r1;
  F1[0] -= fx;
  F2[0] += fx;
  fx = f * r[1] / r1;
  F1[1] -= fx;
  F2[1] += fx;
  fx = f * r[2] / r1;
  F1[2] -= fx;
  F2[2] += fx;
}

void PotentialBase::restrainAngle(double const r1[], double const r2[],
                                  double const r3[], double f1[], double f2[],
                                  double f3[], double &u, const double k,
                                  const double aeq) {
  double r12[3], r12_1, r12_2, r23[3], r23_1, r23_2, cosa, a, m, d12, d23;
  distance(r2, r1, r12, r12_1, r12_2);
  distance(r3, r2, r23, r23_1, r23_2);
  cosa = -dotProduct(r12, r23) / r12_1 / r23_1;
  a = std::acos(cosa);
  u += k * (a - aeq) * (a - aeq);
  m = -2 * k * (a - aeq);               // moment
  m /= -1 * std::sqrt(1 - cosa * cosa); //*d(a)/d(cos a)
  for (int i = 0; i < 3; ++i) {
    d12 = -cosa * r12[i] / r12_2 - r23[i] / r12_1 / r23_1; // d(cos a)/d(r12)
    d23 = -cosa * r23[i] / r23_2 - r12[i] / r12_1 / r23_1; // d(cos a)/d(r23)
    f1[i] -= m * d12;
    f2[i] += m * (d12 - d23);
    f3[i] += m * d23;
  };
}
 
const double PotentialBase::ONE_OVER_4_PI_EPSILON0 =
    14.399644532010862; // eV e^2 / Angstrom

void PotentialBase::calculateCentre(double const r1[], double const r2[],
                                    double rc[]) {
  for (int i = 0; i < 3; ++i) {
    rc[i] = (r1[i] + unBreak1(r2[i], r1[i], i)) * 0.5;
  };
}

void PotentialBase::calculateCentre(double const r1[], double const r2[],
                                    double const r3[], double rc[]) {
  for (int i = 0; i < 3; ++i) {
    rc[i] =
        (r1[i] + unBreak1(r2[i], r1[i], i) + unBreak1(r3[i], r1[i], i)) / 3.0;
  };
}

void PotentialBase::calculateWeightedCentre(double const w1, double const w2,
                                            double const w3, double const r1[],
                                            double const r2[],
                                            double const r3[], double rc[]) {
  assert(std::fabs(w1 + w2 + w3 - 1.0) < 1e-9);
  for (int i = 0; i < 3; ++i) {
    rc[i] = w1 * r1[i] + w2 * unBreak1(r2[i], r1[i], i) +
            w3 * unBreak1(r3[i], r1[i], i);
  };
}

void PotentialBase::coulombWithCutoff(const double r1[], const double r2[],
                                      double f1[], double f2[], double &energy,
                                      double qq) {
  double d1, d[3];
  // d: vector distance (r1-r2), d1 distance and d2 squared distance
  distance(r1, r2, d, d1);
  if (d1 < cutoff_) {
    double f, en;
    coulomb(d1, f, en, qq);
    switching(d1, f, en);
    for (int l = 0; l < 3; ++l) {
      f1[l] += f * d[l] / d1;
      f2[l] -= f * d[l] / d1;
    };
    energy += en;
  }
}

void PotentialBase::coulomb(const double r1[], const double r2[], double f1[],
                            double f2[], double &energy, double qq) {
  double d1, d[3], f, en;
  // d: vector distance (r1-r2), d1 distance and d2 squared distance
  distance(r1, r2, d, d1);
  coulomb(d1, f, en, qq);
  for (int l = 0; l < 3; ++l) {
    f1[l] += f * d[l] / d1;
    f2[l] -= f * d[l] / d1;
  };
  energy += en;
}

void PotentialBase::coulomb(double distance, double &force, double &energy,
                            double const qq) {
  energy = ONE_OVER_4_PI_EPSILON0 * qq / distance;
  force = energy / distance;
}

void PotentialBase::spreadForce(double f1[], double f2[], double const fc[]) {
  for (int i = 0; i < 3; ++i) {
    double const f = fc[i] / 2;
    f1[i] += f;
    f2[i] += f;
  };
}

void PotentialBase::spreadForce(double f1[], double f2[], double f3[],
                                double const fc[]) {
  for (int i = 0; i < 3; ++i) {
    double const f = fc[i] / 3;
    f1[i] += f;
    f2[i] += f;
    f3[i] += f;
  };
}

void PotentialBase::spreadWeightedForce(double const w1, double const w2,
                                        double const w3, double f1[],
                                        double f2[], double f3[],
                                        double const fc[]) {
  for (int i = 0; i < 3; ++i) {
    double const f = fc[i];
    f1[i] += f * w1;
    f2[i] += f * w2;
    f3[i] += f * w3;
  };
}

double PotentialBase::epsilon(double const A, double const B) {
  return B * B / A / 4.0;
}

void PotentialBase::lennardJones(double distance, double &force, double &energy,
                                 double const epsilon, double const sigma) {
  double x = sigma * sigma / (distance * distance);
  x *= x * x;
  energy = 4 * epsilon * (x - 1) * x;
  force = 24 * epsilon * (2 * x - 1) * x / distance;
}

void PotentialBase::lennardJones(const double R1[], const double R2[],
                                 double F1[], double F2[], double &E,
                                 double const epsilon, double const sigma) {
  // WARNING: F1 and F2 are incremented.
  double F, Fx, en;
  double R12[3];
  double L, L2;
  distance(R2, R1, R12, L, L2);
  lennardJones(L, F, en, epsilon, sigma);
  E += en;
  Fx = F * R12[0] / L;
  F1[0] -= Fx;
  F2[0] += Fx;
  Fx = F * R12[1] / L;
  F1[1] -= Fx;
  F2[1] += Fx;
  Fx = F * R12[2] / L;
  F1[2] -= Fx;
  F2[2] += Fx;
}

void PotentialBase::lennardJonesWithCutoff(double const r1[], double const r2[],
                                           double f1[], double f2[],
                                           double &energy, double const epsilon,
                                           double const sigma) {
  double r12[3] = {0}, d = 0.0, f, en; // distance
  distance(r2, r1, r12, d);
  if (d < cutoff_) {
    lennardJones(d, f, en, epsilon, sigma);
    switching(d, f, en);
    for (int i = 0; i < 3; ++i) {
      f1[i] -= f * r12[i] / d;
      f2[i] += f * r12[i] / d;
    };
    energy += en;
  };
}

double PotentialBase::sigma(double const A, double const B) {
  return std::pow(A / B, 1.0 / 6.0);
}

double PotentialBase::smithKongEpsilon(double sigma1, double epsilon1,
                                       double sigma2, double epsilon2) {
  double n, d;
  double s1_2 = sigma1 * sigma1, sa6 = s1_2 * s1_2 * s1_2;
  double s2_2 = sigma2 * sigma2, sb6 = s2_2 * s2_2 * s2_2;
  n = 8192.0 * epsilon1 * sa6 * epsilon2 * sb6; // 2^13 = 8192
  d = std::pow(epsilon1 * sa6 * sa6, 1.0 / 13.0) +
      std::pow(epsilon2 * sb6 * sb6, 1.0 / 13.0);
  d = std::pow(d, 13);
  return n / d;
}

double PotentialBase::smithKongSigma(double sigma1, double epsilon1,
                                     double sigma2, double epsilon2) {
  double s1_2 = sigma1 * sigma1;
  double const sa6 = s1_2 * s1_2 * s1_2;
  double s2_2 = sigma2 * sigma2;
  double const sb6 = s2_2 * s2_2 * s2_2;
  double n, d;
  n = std::pow(epsilon1 * sa6 * sa6, 1.0 / 13.0) +
      std::pow(epsilon2 * sb6 * sb6, 1.0 / 13.0);
  n = std::pow(n, 13.0);
  d = 8192.0 * std::sqrt(epsilon1 * sa6 * epsilon2 * sb6); // 2^13 = 8192
  return std::pow(n / d, 1.0 / 6.0);
}

void PotentialBase::switching(double const distance, double &force,
                              double &energy) {
  double x = (distance - cutoff_ + switchingWidth_) / switchingWidth_;
  if (x >= 1.0) { // if more: cutoff
    force = 0.0;
    energy = 0.0;
  }
  if (x > 0.0) { // if more: switching zone
    double S = (2 * x - 3) * x * x + 1;
    double dS = 6 * x * (x - 1);
    double const u = energy;
    force = force * S - u * dS / switchingWidth_;
    energy *= S;
  };
  // else: full interaction, leave force and energy as the are.
}

void PotentialBase::switching(double const r1[], double const r2[], double f1[],
                              double f2[], double &energy) {
  double z[3], z1;
  distance(r1, r2, z, z1);
  double x = (z1 - cutoff_ + switchingWidth_) / switchingWidth_;
  assert(x > 0.0);
  assert(x < 1.0);
  double S = (2 * x - 3) * x * x + 1;
  double dS = 6 * x * (x - 1);
  double const u = energy;
  for (int k = 0; k < 3; ++k) {
    f1[k] = f1[k] * S - u * dS / switchingWidth_ * z[k] / z1;
    f2[k] = f2[k] * S + u * dS / switchingWidth_ * z[k] / z1;
  };
  energy *= S;
}

void PotentialBase::applyPeriodicity1(double r[]) {
  applyPeriodicity0(r, periods_);
}

double PotentialBase::applyPeriodicity1(double r, int const axis) {
  return applyPeriodicity0(r, periods_[axis]);
}

void PotentialBase::distance(const double x[], const double y[], double z[],
                             double &z1, double &z2) {
  assert(x);
  assert(y);
  assert(z);
  for (int i = 0; i < 3; ++i)
    z[i] = x[i] - y[i];
  applyPeriodicity1(z);
  z2 = z[0] * z[0] + z[1] * z[1] + z[2] * z[2];
  z1 = std::sqrt(z2);
}

void PotentialBase::distance(const double x[], const double y[], double z[],
                             double &z1) {
  double z2;
  distance(x, y, z, z1, z2);
}

void PotentialBase::distance(const double x[], const double y[], double z[]) {
  double z1;
  distance(x, y, z, z1);
}

void PotentialBase::distance(const double x[], const double y[], double &z1) {
  double z[3];
  distance(x, y, z, z1);
}

void PotentialBase::distance(const double x[], const double y[], Vector3 &z) {
  distance(x, y, z.v, z._1, z._2);
}

void PotentialBase::setPeriodicity(const double periods[]) {
  periods_[0] = periods[0];
  periods_[1] = periods[1];
  periods_[2] = periods[2];
}

double PotentialBase::unBreak0(double const r, double const ref,
                               double const period) {
  return applyPeriodicity0(r - ref, period) + ref;
}

void PotentialBase::unBreak0(double r[], double const ref[],
                             double const periods[]) {
  for (int a = 0; a < 3; ++a)
    r[a] = unBreak0(r[a], ref[a], periods[a]);
}

void PotentialBase::unBreak1(double r[], double const ref[]) {
  unBreak0(r, ref, periods_);
}

double PotentialBase::unBreak1(double const r, double const ref,
                               int const axis) {
  return unBreak0(r, ref, periods_[axis]);
}

void PotentialBase::computePt(int const nAtoms, double positions[],
                              double forces[], double &energy,
                              double const periods[], bool const fixed[]) {
  int const nCoord = 3 * nAtoms;
  for (int i = 0; i < nCoord; ++i)
    forces[i] = 0.0;
  energy = 0.0;
  double const(*r)[3] = reinterpret_cast<double const(*)[3]>(positions);
  double (*f)[3] = reinterpret_cast<double (*)[3]>(forces);
  setPeriodicity(periods); // Set Periodic boundaries. Essential in order for
                           // some functions to work.
  for (int i = nAtoms - 1; i > 0; --i) {
    for (int j = i - 1; j >= 0; --j) {
      if (not fixed[i] and not fixed[j]) {
        lennardJonesWithCutoff(r[i], r[j], f[i], f[j], energy, EPSILON_PT,
                               SIGMA_PT);
      };
    };
  };
}

const double PotentialBase::EPSILON_PT = 0.68165797577788501; // eV
const double PotentialBase::SIGMA_PT = 2.54;                  // Angstrom
/*
 const double PotentialBase::ONE_OVER_4_PI_EPSILON0 =
 unit_system::ONE_OVER_4_PI_EPSILON0; const double
 PotentialBase::EPSILON_PT=65.77*KJ_PER_MOL; const double
 PotentialBase::SIGMA_PT=0.254*NM;
 */