ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
9cfb8228470fad0f8afa01fb6bf185710c9d96c4
lammps/src/INTERLAYER/pair_drip.cpp
Axel Kohlmeyer 340207988c fix a couple more bugs like in 5246cedda6
2021-10-05 10:36:25 -04:00

987 lines
30 KiB
C++
Raw Blame History

/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Mingjian Wen (University of Minnesota)
e-mail: wenxx151@umn.edu, wenxx151@gmail.com
This implements the DRIP model as described in
M. Wen, S. Carr, S. Fang, E. Kaxiras, and E. B. Tadmor,
Phys. Rev. B, 98, 235404 (2018).
------------------------------------------------------------------------- */
#include "pair_drip.h"
#include "atom.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "memory.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "neighbor.h"
#include "potential_file_reader.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
#define MAXLINE 1024
#define DELTA 4
#define HALF 0.5
// inline functions
static inline double dot(double const *x, double const *y)
{
return x[0] * y[0] + x[1] * y[1] + x[2] * y[2];
}
static inline void mat_dot_vec(PairDRIP::V3 const *X, double const *y, double *const z)
{
for (int k = 0; k < 3; k++) { z[k] = X[k][0] * y[0] + X[k][1] * y[1] + X[k][2] * y[2]; }
}
/* ---------------------------------------------------------------------- */
PairDRIP::PairDRIP(LAMMPS *lmp) : Pair(lmp)
{
single_enable = 0;
restartinfo = 0;
one_coeff = 1;
manybody_flag = 1;
centroidstressflag = CENTROID_NOTAVAIL;
unit_convert_flag = utils::get_supported_conversions(utils::ENERGY);
params = nullptr;
nearest3neigh = nullptr;
cutmax = 0.0;
}
/* ---------------------------------------------------------------------- */
PairDRIP::~PairDRIP()
{
if (allocated) {
memory->destroy(setflag);
memory->destroy(cutsq);
}
memory->destroy(params);
memory->destroy(elem2param);
memory->destroy(nearest3neigh);
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void PairDRIP::init_style()
{
if (force->newton_pair == 0) error->all(FLERR, "Pair style drip requires newton pair on");
if (!atom->molecule_flag) error->all(FLERR, "Pair style drip requires atom attribute molecule");
// need a full neighbor list, including neighbors of ghosts
int irequest = neighbor->request(this, instance_me);
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->full = 1;
neighbor->requests[irequest]->ghost = 1;
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
void PairDRIP::allocate()
{
allocated = 1;
int n = atom->ntypes + 1;
memory->create(setflag, n, n, "pair:setflag");
memory->create(cutsq, n, n, "pair:cutsq");
map = new int[n];
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairDRIP::settings(int narg, char ** /* arg */)
{
if (narg != 0) error->all(FLERR, "Illegal pair_style command");
if (!utils::strmatch(force->pair_style, "^hybrid/overlay"))
error->all(FLERR, "Pair style drip must be used as sub-style with hybrid/overlay");
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairDRIP::coeff(int narg, char **arg)
{
if (!allocated) allocate();
map_element2type(narg - 3, arg + 3);
read_file(arg[2]);
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairDRIP::init_one(int i, int j)
{
if (setflag[i][j] == 0) error->all(FLERR, "All pair coeffs are not set");
int itype = map[i];
int jtype = map[j];
int iparam_ij = elem2param[itype][jtype];
Param &p = params[iparam_ij];
// max cutoff is the main cutoff plus the normal cutoff such that
return p.rcut + p.ncut;
}
/* ----------------------------------------------------------------------
read DRIP file
------------------------------------------------------------------------- */
void PairDRIP::read_file(char *filename)
{
memory->sfree(params);
params = nullptr;
nparams = maxparam = 0;
// open file on proc 0
if (comm->me == 0) {
PotentialFileReader reader(lmp, filename, "drip", unit_convert_flag);
char *line;
// transparently convert units for supported conversions
int unit_convert = reader.get_unit_convert();
double conversion_factor = utils::get_conversion_factor(utils::ENERGY, unit_convert);
while ((line = reader.next_line(NPARAMS_PER_LINE))) {
try {
ValueTokenizer values(line);
std::string iname = values.next_string();
std::string jname = values.next_string();
// ielement,jelement = 1st args
// if both args are in element list, then parse this line
// else skip to next entry in file
int ielement, jelement;
for (ielement = 0; ielement < nelements; ielement++)
if (iname == elements[ielement]) break;
if (ielement == nelements) continue;
for (jelement = 0; jelement < nelements; jelement++)
if (jname == elements[jelement]) break;
if (jelement == nelements) continue;
// expand storage, if needed
if (nparams == maxparam) {
maxparam += DELTA;
params = (Param *) memory->srealloc(params, maxparam * sizeof(Param), "pair:params");
// make certain all addional allocated storage is initialized
// to avoid false positives when checking with valgrind
memset(params + nparams, 0, DELTA * sizeof(Param));
}
params[nparams].ielement = ielement;
params[nparams].jelement = jelement;
params[nparams].C0 = values.next_double();
params[nparams].C2 = values.next_double();
params[nparams].C4 = values.next_double();
params[nparams].C = values.next_double();
params[nparams].delta = values.next_double();
params[nparams].lambda = values.next_double();
params[nparams].A = values.next_double();
params[nparams].z0 = values.next_double();
params[nparams].B = values.next_double();
params[nparams].eta = values.next_double();
params[nparams].rhocut = values.next_double();
params[nparams].rcut = values.next_double();
params[nparams].ncut = values.next_double();
} catch (TokenizerException &e) {
error->one(FLERR, e.what());
}
if (unit_convert) {
params[nparams].C0 *= conversion_factor;
params[nparams].C2 *= conversion_factor;
params[nparams].C4 *= conversion_factor;
params[nparams].C *= conversion_factor;
params[nparams].A *= conversion_factor;
params[nparams].B *= conversion_factor;
}
// convenient precomputations
params[nparams].rhocutsq = params[nparams].rhocut * params[nparams].rhocut;
params[nparams].rcutsq = params[nparams].rcut * params[nparams].rcut;
params[nparams].ncutsq = params[nparams].ncut * params[nparams].ncut;
nparams++;
}
}
MPI_Bcast(&nparams, 1, MPI_INT, 0, world);
MPI_Bcast(&maxparam, 1, MPI_INT, 0, world);
if (comm->me != 0) {
params = (Param *) memory->srealloc(params, maxparam * sizeof(Param), "pair:params");
}
MPI_Bcast(params, maxparam * sizeof(Param), MPI_BYTE, 0, world);
memory->destroy(elem2param);
memory->create(elem2param, nelements, nelements, "pair:elem2param");
for (int i = 0; i < nelements; i++) {
for (int j = 0; j < nelements; j++) {
int n = -1;
for (int m = 0; m < nparams; m++) {
if (i == params[m].ielement && j == params[m].jelement) {
if (n >= 0) error->all(FLERR, "DRIP potential file has duplicate entry");
n = m;
}
}
if (n < 0) error->all(FLERR, "Potential file is missing an entry");
elem2param[i][j] = n;
}
}
}
/* ---------------------------------------------------------------------- */
void PairDRIP::compute(int eflag, int vflag)
{
int i, j, ii, jj, inum, jnum, itype, jtype;
double xtmp, ytmp, ztmp, delx, dely, delz, evdwl, rsq;
int *ilist, *jlist, *numneigh, **firstneigh;
double ni[3];
double dni_dri[3][3], dni_drnb1[3][3];
double dni_drnb2[3][3], dni_drnb3[3][3];
ev_init(eflag, vflag);
double **x = atom->x;
double **f = atom->f;
int *type = atom->type;
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
find_nearest3neigh();
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
if (nearest3neigh[i][0] == -1) { continue; }
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype = map[type[i]];
jlist = firstneigh[i];
jnum = numneigh[i];
// normal and its derivatives w.r.t. atom i and its 3 nearest neighbors
calc_normal(i, ni, dni_dri, dni_drnb1, dni_drnb2, dni_drnb3);
double fi[3] = {0., 0., 0.};
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
if (nearest3neigh[j][0] == -1) { continue; }
jtype = map[type[j]];
delx = x[j][0] - xtmp;
dely = x[j][1] - ytmp;
delz = x[j][2] - ztmp;
rsq = delx * delx + dely * dely + delz * delz;
int iparam_ij = elem2param[itype][jtype];
Param &p = params[iparam_ij];
double rcutsq = p.rcutsq;
// only include the interaction between different layers
if (rsq < rcutsq && atom->molecule[i] != atom->molecule[j]) {
double fj[3] = {0., 0., 0.};
double rvec[3] = {delx, dely, delz};
double phi_attr = calc_attractive(p, rsq, rvec, fi, fj);
double phi_repul = calc_repulsive(i, j, p, rsq, rvec, ni, dni_dri, dni_drnb1, dni_drnb2,
dni_drnb3, fi, fj);
if (eflag)
evdwl = HALF * (phi_repul + phi_attr);
else
evdwl = 0.0;
if (evflag) ev_tally(i, j, nlocal, newton_pair, evdwl, 0.0, 0, 0, 0, 0);
f[j][0] += fj[0];
f[j][1] += fj[1];
f[j][2] += fj[2];
if (vflag_either) v_tally2_newton(j, fj, x[j]);
}
} //loop over jj
f[i][0] += fi[0];
f[i][1] += fi[1];
f[i][2] += fi[2];
if (vflag_either) v_tally2_newton(i, fi, x[i]);
} // loop over ii
if (vflag_fdotr) virial_fdotr_compute();
}
/* ----------------------------------------------------------------------
Attractive part, i.e. the r^(-6) part
------------------------------------------------------------------------- */
double PairDRIP::calc_attractive(Param &p, double const rsq, double const *rvec, double *const fi,
double *const fj)
{
double const z0 = p.z0;
double const A = p.A;
double const cutoff = p.rcut;
double const r = sqrt(rsq);
double roz0_sq = rsq / (z0 * z0);
double dtp;
double tp = tap(r, cutoff, dtp);
double r6 = A / (roz0_sq * roz0_sq * roz0_sq);
double dr6 = -6 * r6 / r;
double phi = -r6 * tp;
double fpair = -HALF * (r6 * dtp + dr6 * tp);
fi[0] += rvec[0] * fpair / r;
fi[1] += rvec[1] * fpair / r;
fi[2] += rvec[2] * fpair / r;
fj[0] -= rvec[0] * fpair / r;
fj[1] -= rvec[1] * fpair / r;
fj[2] -= rvec[2] * fpair / r;
return phi;
}
/* ----------------------------------------------------------------------
Repulsive part that depends on transverse distance and dihedral angle
------------------------------------------------------------------------- */
double PairDRIP::calc_repulsive(int const i, int const j, Param &p, double const rsq,
double const *rvec, double const *ni, V3 const *dni_dri,
V3 const *dni_drnb1, V3 const *dni_drnb2, V3 const *dni_drnb3,
double *const fi, double *const fj)
{
double **f = atom->f;
double **x = atom->x;
double C0 = p.C0;
double C2 = p.C2;
double C4 = p.C4;
double C = p.C;
double delta = p.delta;
double lambda = p.lambda;
double z0 = p.z0;
double cutoff = p.rcut;
// nearest 3 neighbors of atoms i and j
int nbi1 = nearest3neigh[i][0];
int nbi2 = nearest3neigh[i][1];
int nbi3 = nearest3neigh[i][2];
int nbj1 = nearest3neigh[j][0];
int nbj2 = nearest3neigh[j][1];
int nbj3 = nearest3neigh[j][2];
double fnbi1[3];
double fnbi2[3];
double fnbi3[3];
double fnbj1[3];
double fnbj2[3];
double fnbj3[3];
V3 dgij_dri;
V3 dgij_drj;
V3 dgij_drk1;
V3 dgij_drk2;
V3 dgij_drk3;
V3 dgij_drl1;
V3 dgij_drl2;
V3 dgij_drl3;
V3 drhosqij_dri;
V3 drhosqij_drj;
V3 drhosqij_drnb1;
V3 drhosqij_drnb2;
V3 drhosqij_drnb3;
double r = sqrt(rsq);
// derivative of rhosq w.r.t. atoms i j and the nearests 3 neighs of i
get_drhosqij(rvec, ni, dni_dri, dni_drnb1, dni_drnb2, dni_drnb3, drhosqij_dri, drhosqij_drj,
drhosqij_drnb1, drhosqij_drnb2, drhosqij_drnb3);
// transverse decay function f(rho) and its derivative w.r.t. rhosq
double rhosqij;
double dtdij;
double tdij = td(C0, C2, C4, delta, rvec, r, ni, rhosqij, dtdij);
// dihedral angle function and its derivateives
double dgij_drhosq;
double gij = dihedral(i, j, p, rhosqij, dgij_drhosq, dgij_dri, dgij_drj, dgij_drk1, dgij_drk2,
dgij_drk3, dgij_drl1, dgij_drl2, dgij_drl3);
double V2 = C + tdij + gij;
// tap part
double dtp;
double tp = tap(r, cutoff, dtp);
// exponential part
double V1 = exp(-lambda * (r - z0));
double dV1 = -V1 * lambda;
// total energy
double phi = tp * V1 * V2;
for (int k = 0; k < 3; k++) {
// forces due to derivatives of tap and V1
double tmp = HALF * (dtp * V1 + tp * dV1) * V2 * rvec[k] / r;
fi[k] += tmp;
fj[k] -= tmp;
// contributions from transverse decay part tdij and the dihedral part gij
// derivative of V2 contribute to atoms i, j
fi[k] -= HALF * tp * V1 * ((dtdij + dgij_drhosq) * drhosqij_dri[k] + dgij_dri[k]);
fj[k] -= HALF * tp * V1 * ((dtdij + dgij_drhosq) * drhosqij_drj[k] + dgij_drj[k]);
// derivative of V2 contribute to nearest 3 neighs of atom i
fnbi1[k] = -HALF * tp * V1 * ((dtdij + dgij_drhosq) * drhosqij_drnb1[k] + dgij_drk1[k]);
fnbi2[k] = -HALF * tp * V1 * ((dtdij + dgij_drhosq) * drhosqij_drnb2[k] + dgij_drk2[k]);
fnbi3[k] = -HALF * tp * V1 * ((dtdij + dgij_drhosq) * drhosqij_drnb3[k] + dgij_drk3[k]);
// derivative of V2 contribute to nearest 3 neighs of atom j
fnbj1[k] = -HALF * tp * V1 * dgij_drl1[k];
fnbj2[k] = -HALF * tp * V1 * dgij_drl2[k];
fnbj3[k] = -HALF * tp * V1 * dgij_drl3[k];
}
for (int k = 0; k < 3; k++) {
f[nbi1][k] += fnbi1[k];
f[nbi2][k] += fnbi2[k];
f[nbi3][k] += fnbi3[k];
f[nbj1][k] += fnbj1[k];
f[nbj2][k] += fnbj2[k];
f[nbj3][k] += fnbj3[k];
}
if (vflag_either) {
v_tally2_newton(nbi1, fnbi1, x[nbi1]);
v_tally2_newton(nbi2, fnbi2, x[nbi2]);
v_tally2_newton(nbi3, fnbi3, x[nbi3]);
v_tally2_newton(nbj1, fnbj1, x[nbj1]);
v_tally2_newton(nbj2, fnbj2, x[nbj2]);
v_tally2_newton(nbj3, fnbj3, x[nbj3]);
}
return phi;
}
/* ---------------------------------------------------------------------- */
void PairDRIP::find_nearest3neigh()
{
int i, j, ii, jj, allnum, inum, jnum, itype, jtype, size;
double xtmp, ytmp, ztmp, delx, dely, delz, rsq;
int *ilist, *jlist, *numneigh, **firstneigh;
double **x = atom->x;
int *type = atom->type;
allnum = list->inum + list->gnum;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
size = allnum;
memory->destroy(nearest3neigh);
memory->create(nearest3neigh, size, 3, "pair:nearest3neigh");
for (ii = 0; ii < allnum; ii++) {
i = ilist[ii];
// If "NULL" used in pair_coeff, i could be larger than allnum
if (i >= size) {
size = i + 1;
memory->grow(nearest3neigh, size, 3, "pair:nearest3neigh");
}
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype = map[type[i]];
jlist = firstneigh[i];
jnum = numneigh[i];
// init nb1 to be the 1st nearest neigh, nb3 the 3rd nearest
int nb1 = -1;
int nb2 = -1;
int nb3 = -1;
double nb1_rsq = 1.0e10 + 1;
double nb2_rsq = 2.0e10;
double nb3_rsq = 3.0e10;
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
jtype = map[type[j]];
delx = x[j][0] - xtmp;
dely = x[j][1] - ytmp;
delz = x[j][2] - ztmp;
rsq = delx * delx + dely * dely + delz * delz;
int iparam_ij = elem2param[itype][jtype];
double ncutsq = params[iparam_ij].ncutsq;
if (rsq < ncutsq && atom->molecule[i] == atom->molecule[j]) {
// find the 3 nearest neigh
if (rsq < nb1_rsq) {
nb3 = nb2;
nb2 = nb1;
nb1 = j;
nb3_rsq = nb2_rsq;
nb2_rsq = nb1_rsq;
nb1_rsq = rsq;
} else if (rsq < nb2_rsq) {
nb3 = nb2;
nb2 = j;
nb3_rsq = nb2_rsq;
nb2_rsq = rsq;
} else if (rsq < nb3_rsq) {
nb3 = j;
nb3_rsq = rsq;
}
}
} // loop over jj
// store neighbors to be used later to compute normal
if (nb3_rsq >= 1.0e10) {
if (i < inum) {
error->one(FLERR,
"No enough neighbors to construct normal. Check the "
"configuration to see whether atoms fly away.");
} else {
// This only happens for ghost atoms that are near the boundary of the
// domain (i.e. r > r_cut + n_cut). These ghost atoms will not be
// the i j atoms in the compute function, but only neighbors of j atoms.
// It is allowed not to have three neighbors for these atoms, since
// their normals are not needed.
nearest3neigh[i][0] = -1;
nearest3neigh[i][1] = -1;
nearest3neigh[i][2] = -1;
}
} else {
nearest3neigh[i][0] = nb1;
nearest3neigh[i][1] = nb2;
nearest3neigh[i][2] = nb3;
}
} // loop over ii
}
/* ---------------------------------------------------------------------- */
void PairDRIP::calc_normal(int const i, double *const normal, V3 *const dn_dri, V3 *const dn_drk1,
V3 *const dn_drk2, V3 *const dn_drk3)
{
int k1 = nearest3neigh[i][0];
int k2 = nearest3neigh[i][1];
int k3 = nearest3neigh[i][2];
// normal does not depend on i, setting to zero
for (int j = 0; j < 3; j++) {
for (int k = 0; k < 3; k++) { dn_dri[j][k] = 0.0; }
}
// get normal and derives of normal w.r.t to its 3 nearest neighbors
double **x = atom->x;
deriv_cross(x[k1], x[k2], x[k3], normal, dn_drk1, dn_drk2, dn_drk3);
}
/* ---------------------------------------------------------------------- */
void PairDRIP::get_drhosqij(double const *rij, double const *ni, V3 const *dni_dri,
V3 const *dni_drn1, V3 const *dni_drn2, V3 const *dni_drn3,
double *const drhosq_dri, double *const drhosq_drj,
double *const drhosq_drn1, double *const drhosq_drn2,
double *const drhosq_drn3)
{
int k;
double ni_dot_rij = 0;
double dni_dri_dot_rij[3];
double dni_drn1_dot_rij[3];
double dni_drn2_dot_rij[3];
double dni_drn3_dot_rij[3];
ni_dot_rij = dot(ni, rij);
mat_dot_vec(dni_dri, rij, dni_dri_dot_rij);
mat_dot_vec(dni_drn1, rij, dni_drn1_dot_rij);
mat_dot_vec(dni_drn2, rij, dni_drn2_dot_rij);
mat_dot_vec(dni_drn3, rij, dni_drn3_dot_rij);
for (k = 0; k < 3; k++) {
drhosq_dri[k] = -2 * rij[k] - 2 * ni_dot_rij * (-ni[k] + dni_dri_dot_rij[k]);
drhosq_drj[k] = 2 * rij[k] - 2 * ni_dot_rij * ni[k];
drhosq_drn1[k] = -2 * ni_dot_rij * dni_drn1_dot_rij[k];
drhosq_drn2[k] = -2 * ni_dot_rij * dni_drn2_dot_rij[k];
drhosq_drn3[k] = -2 * ni_dot_rij * dni_drn3_dot_rij[k];
}
}
/* ----------------------------------------------------------------------
derivartive of transverse decay function f(rho) w.r.t. rho
------------------------------------------------------------------------- */
double PairDRIP::td(double C0, double C2, double C4, double delta, double const *const rvec,
double r, const double *const n, double &rho_sq, double &dtd)
{
double n_dot_r = dot(n, rvec);
rho_sq = r * r - n_dot_r * n_dot_r;
// in case n is [0, 0, 1] and rho_sq is negative due to numerical error
if (rho_sq < 0) { rho_sq = 0; }
double del_sq = delta * delta;
double rod_sq = rho_sq / del_sq;
double td = exp(-rod_sq) * (C0 + rod_sq * (C2 + rod_sq * C4));
dtd = -td / del_sq + exp(-rod_sq) * (C2 + 2 * C4 * rod_sq) / del_sq;
return td;
}
/* ----------------------------------------------------------------------
derivartive of dihedral angle func gij w.r.t rho, and atom positions
------------------------------------------------------------------------- */
double PairDRIP::dihedral(const int i, const int j, Param &p, double const rhosq, double &d_drhosq,
double *const d_dri, double *const d_drj, double *const d_drk1,
double *const d_drk2, double *const d_drk3, double *const d_drl1,
double *const d_drl2, double *const d_drl3)
{
double **x = atom->x;
// get parameter
double B = p.B;
double eta = p.eta;
double cut_rhosq = p.rhocutsq;
// local vars
double cos_kl[3][3]; // cos_omega_k1ijl1, cos_omega_k1ijl2 ...
double d_dcos_kl[3][3]; // deriv of dihedral w.r.t to cos_omega_kijl
double dcos_kl[3][3][4][3]; // 4 indicates k, i, j, l. e.g. dcoskl[0][1][0]
// means dcos_omega_k1ijl2 / drk
// if larger than cutoff of rho, return 0
if (rhosq >= cut_rhosq) {
d_drhosq = 0;
for (int dim = 0; dim < 3; dim++) {
d_dri[dim] = 0;
d_drj[dim] = 0;
d_drk1[dim] = 0;
d_drk2[dim] = 0;
d_drk3[dim] = 0;
d_drl1[dim] = 0;
d_drl2[dim] = 0;
d_drl3[dim] = 0;
}
double dihe = 0.0;
return dihe;
}
// 3 neighs of atoms i and j
int k[3];
int l[3];
for (int m = 0; m < 3; m++) {
k[m] = nearest3neigh[i][m];
l[m] = nearest3neigh[j][m];
}
// cos_omega_kijl and the derivatives w.r.t coordinates
for (int m = 0; m < 3; m++) {
for (int n = 0; n < 3; n++) {
cos_kl[m][n] = deriv_cos_omega(x[k[m]], x[i], x[j], x[l[n]], dcos_kl[m][n][0],
dcos_kl[m][n][1], dcos_kl[m][n][2], dcos_kl[m][n][3]);
}
}
double epart1 = exp(-eta * cos_kl[0][0] * cos_kl[0][1] * cos_kl[0][2]);
double epart2 = exp(-eta * cos_kl[1][0] * cos_kl[1][1] * cos_kl[1][2]);
double epart3 = exp(-eta * cos_kl[2][0] * cos_kl[2][1] * cos_kl[2][2]);
double D2 = epart1 + epart2 + epart3;
// cutoff function
double d_drhosq_tap;
double D0 = B * tap_rho(rhosq, cut_rhosq, d_drhosq_tap);
// dihedral energy
double dihe = D0 * D2;
// deriv of dihedral w.r.t rhosq
d_drhosq = B * d_drhosq_tap * D2;
// deriv of dihedral w.r.t cos_omega_kijl
d_dcos_kl[0][0] = -D0 * epart1 * eta * cos_kl[0][1] * cos_kl[0][2];
d_dcos_kl[0][1] = -D0 * epart1 * eta * cos_kl[0][0] * cos_kl[0][2];
d_dcos_kl[0][2] = -D0 * epart1 * eta * cos_kl[0][0] * cos_kl[0][1];
d_dcos_kl[1][0] = -D0 * epart2 * eta * cos_kl[1][1] * cos_kl[1][2];
d_dcos_kl[1][1] = -D0 * epart2 * eta * cos_kl[1][0] * cos_kl[1][2];
d_dcos_kl[1][2] = -D0 * epart2 * eta * cos_kl[1][0] * cos_kl[1][1];
d_dcos_kl[2][0] = -D0 * epart3 * eta * cos_kl[2][1] * cos_kl[2][2];
d_dcos_kl[2][1] = -D0 * epart3 * eta * cos_kl[2][0] * cos_kl[2][2];
d_dcos_kl[2][2] = -D0 * epart3 * eta * cos_kl[2][0] * cos_kl[2][1];
// initialization to be zero and later add values
for (int dim = 0; dim < 3; dim++) {
d_drk1[dim] = 0.;
d_drk2[dim] = 0.;
d_drk3[dim] = 0.;
d_dri[dim] = 0.;
d_drj[dim] = 0.;
d_drl1[dim] = 0.;
d_drl2[dim] = 0.;
d_drl3[dim] = 0.;
}
for (int m = 0; m < 3; m++) {
for (int dim = 0; dim < 3; dim++) {
d_drk1[dim] += d_dcos_kl[0][m] * dcos_kl[0][m][0][dim];
d_drk2[dim] += d_dcos_kl[1][m] * dcos_kl[1][m][0][dim];
d_drk3[dim] += d_dcos_kl[2][m] * dcos_kl[2][m][0][dim];
d_drl1[dim] += d_dcos_kl[m][0] * dcos_kl[m][0][3][dim];
d_drl2[dim] += d_dcos_kl[m][1] * dcos_kl[m][1][3][dim];
d_drl3[dim] += d_dcos_kl[m][2] * dcos_kl[m][2][3][dim];
}
for (int n = 0; n < 3; n++) {
for (int dim = 0; dim < 3; dim++) {
d_dri[dim] += d_dcos_kl[m][n] * dcos_kl[m][n][1][dim];
d_drj[dim] += d_dcos_kl[m][n] * dcos_kl[m][n][2][dim];
}
}
}
return dihe;
}
/* ----------------------------------------------------------------------
compute cos(omega_kijl) and the derivateives
------------------------------------------------------------------------- */
double PairDRIP::deriv_cos_omega(double const *rk, double const *ri, double const *rj,
double const *rl, double *const dcos_drk, double *const dcos_dri,
double *const dcos_drj, double *const dcos_drl)
{
double ejik[3];
double eijl[3];
double tmp1[3];
double tmp2[3];
double dejik_dri[3][3];
double dejik_drj[3][3];
double dejik_drk[3][3];
double deijl_dri[3][3];
double deijl_drj[3][3];
double deijl_drl[3][3];
// ejik and derivatives
// Note the returned dejik_dri ... are actually the transpose
deriv_cross(ri, rj, rk, ejik, dejik_dri, dejik_drj, dejik_drk);
// flip sign
// deriv_cross computes rij cross rik, here we need rji cross rik
for (int m = 0; m < 3; m++) {
ejik[m] = -ejik[m];
for (int n = 0; n < 3; n++) {
dejik_dri[m][n] = -dejik_dri[m][n];
dejik_drj[m][n] = -dejik_drj[m][n];
dejik_drk[m][n] = -dejik_drk[m][n];
}
}
// eijl and derivatives
deriv_cross(rj, ri, rl, eijl, deijl_drj, deijl_dri, deijl_drl);
// flip sign
for (int m = 0; m < 3; m++) {
eijl[m] = -eijl[m];
for (int n = 0; n < 3; n++) {
deijl_drj[m][n] = -deijl_drj[m][n];
deijl_dri[m][n] = -deijl_dri[m][n];
deijl_drl[m][n] = -deijl_drl[m][n];
}
}
// dcos_drk
mat_dot_vec(dejik_drk, eijl, dcos_drk);
// dcos_dri
mat_dot_vec(dejik_dri, eijl, tmp1);
mat_dot_vec(deijl_dri, ejik, tmp2);
for (int m = 0; m < 3; m++) { dcos_dri[m] = tmp1[m] + tmp2[m]; }
// dcos_drj
mat_dot_vec(dejik_drj, eijl, tmp1);
mat_dot_vec(deijl_drj, ejik, tmp2);
for (int m = 0; m < 3; m++) { dcos_drj[m] = tmp1[m] + tmp2[m]; }
// dcos drl
mat_dot_vec(deijl_drl, ejik, dcos_drl);
// cos_oemga_kijl
double cos_omega = dot(ejik, eijl);
return cos_omega;
}
/* ---------------------------------------------------------------------- */
double PairDRIP::tap(double r, double cutoff, double &dtap)
{
double t;
double r_min = 0;
if (r <= r_min) {
t = 1;
dtap = 0;
} else {
double roc = (r - r_min) / (cutoff - r_min);
double roc_sq = roc * roc;
t = roc_sq * roc_sq * (-35.0 + 84.0 * roc + roc_sq * (-70.0 + 20.0 * roc)) + 1;
dtap =
roc_sq * roc / (cutoff - r_min) * (-140.0 + 420.0 * roc + roc_sq * (-420.0 + 140.0 * roc));
}
return t;
}
/* ---------------------------------------------------------------------- */
double PairDRIP::tap_rho(double rhosq, double cut_rhosq, double &drhosq)
{
double roc_sq;
double roc;
double t;
roc_sq = rhosq / cut_rhosq;
roc = sqrt(roc_sq);
t = roc_sq * roc_sq * (-35.0 + 84.0 * roc + roc_sq * (-70.0 + 20.0 * roc)) + 1;
// Note this dtap/drho_sq not dtap/drho
drhosq = roc_sq / cut_rhosq * (-70.0 + 210.0 * roc + roc_sq * (-210.0 + 70.0 * roc));
return t;
}
/* ----------------------------------------------------------------------
Compute the normalized cross product of two vector rkl, rkm, and the
derivates w.r.t rk, rl, rm.
Note, the returned dcross_drk, dcross_drl, and dcross_drm are actually the
transpose.
------------------------------------------------------------------------- */
void PairDRIP::deriv_cross(double const *rk, double const *rl, double const *rm,
double *const cross, V3 *const dcross_drk, V3 *const dcross_drl,
V3 *const dcross_drm)
{
double x[3];
double y[3];
double p[3];
double q;
double q_cubic;
double d_invq_d_x0;
double d_invq_d_x1;
double d_invq_d_x2;
double d_invq_d_y0;
double d_invq_d_y1;
double d_invq_d_y2;
int i, j;
// get x = rkl and y = rkm
for (i = 0; i < 3; i++) {
x[i] = rl[i] - rk[i];
y[i] = rm[i] - rk[i];
}
// cross product
p[0] = x[1] * y[2] - x[2] * y[1];
p[1] = x[2] * y[0] - x[0] * y[2];
p[2] = x[0] * y[1] - x[1] * y[0];
q = sqrt(p[0] * p[0] + p[1] * p[1] + p[2] * p[2]);
// normalized cross
cross[0] = p[0] / q;
cross[1] = p[1] / q;
cross[2] = p[2] / q;
// compute derivatives
// derivative of inverse q (i.e. 1/q) w.r.t x and y
q_cubic = q * q * q;
d_invq_d_x0 = (+p[1] * y[2] - p[2] * y[1]) / q_cubic;
d_invq_d_x1 = (-p[0] * y[2] + p[2] * y[0]) / q_cubic;
d_invq_d_x2 = (p[0] * y[1] - p[1] * y[0]) / q_cubic;
d_invq_d_y0 = (-p[1] * x[2] + p[2] * x[1]) / q_cubic;
d_invq_d_y1 = (p[0] * x[2] - p[2] * x[0]) / q_cubic;
d_invq_d_y2 = (-p[0] * x[1] + p[1] * x[0]) / q_cubic;
// dcross/drl transposed
dcross_drl[0][0] = p[0] * d_invq_d_x0;
dcross_drl[0][1] = -y[2] / q + p[1] * d_invq_d_x0;
dcross_drl[0][2] = y[1] / q + p[2] * d_invq_d_x0;
dcross_drl[1][0] = y[2] / q + p[0] * d_invq_d_x1;
dcross_drl[1][1] = p[1] * d_invq_d_x1;
dcross_drl[1][2] = -y[0] / q + p[2] * d_invq_d_x1;
dcross_drl[2][0] = -y[1] / q + p[0] * d_invq_d_x2;
dcross_drl[2][1] = y[0] / q + p[1] * d_invq_d_x2;
dcross_drl[2][2] = p[2] * d_invq_d_x2;
// dcross/drm transposed
dcross_drm[0][0] = p[0] * d_invq_d_y0;
dcross_drm[0][1] = x[2] / q + p[1] * d_invq_d_y0;
dcross_drm[0][2] = -x[1] / q + p[2] * d_invq_d_y0;
dcross_drm[1][0] = -x[2] / q + p[0] * d_invq_d_y1;
dcross_drm[1][1] = p[1] * d_invq_d_y1;
dcross_drm[1][2] = x[0] / q + p[2] * d_invq_d_y1;
dcross_drm[2][0] = x[1] / q + p[0] * d_invq_d_y2;
dcross_drm[2][1] = -x[0] / q + p[1] * d_invq_d_y2;
dcross_drm[2][2] = p[2] * d_invq_d_y2;
// dcross/drk transposed
for (i = 0; i < 3; i++) {
for (j = 0; j < 3; j++) { dcross_drk[i][j] = -(dcross_drl[i][j] + dcross_drm[i][j]); }
}
}
Reference in New Issue View Git Blame Copy Permalink