ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

enable and apply clang-format

Browse Source
This commit is contained in:
Axel Kohlmeyer
2022-03-17 19:44:53 -04:00
parent 58030009da
commit 85e581eb2a
5 changed files with 535 additions and 548 deletions
Download Patch File Download Diff File Expand all files Collapse all files

264
src/MOLECULE/dihedral_charmm.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories https://www.lammps.org/, Sandia National Laboratories
@ -32,13 +31,13 @@
#include <cmath> #include <cmath>
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst; using MathConst::DEG2RAD;
#define TOLERANCE 0.05 static constexpr double TOLERANCE = 0.05;
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
DihedralCharmm::DihedralCharmm(LAMMPS *lmp) : Dihedral(lmp) DihedralCharmm::DihedralCharmm(LAMMPS *_lmp) : Dihedral(_lmp)
{ {
weightflag = 0; weightflag = 0;
writedata = 1; writedata = 1;
@ -63,19 +62,19 @@ DihedralCharmm::~DihedralCharmm()
void DihedralCharmm::compute(int eflag, int vflag) void DihedralCharmm::compute(int eflag, int vflag)
{ {
int i1,i2,i3,i4,i,m,n,type; int i1, i2, i3, i4, i, m, n, type;
double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,vb2xm,vb2ym,vb2zm; double vb1x, vb1y, vb1z, vb2x, vb2y, vb2z, vb3x, vb3y, vb3z, vb2xm, vb2ym, vb2zm;
double edihedral,f1[3],f2[3],f3[3],f4[3]; double edihedral, f1[3], f2[3], f3[3], f4[3];
double ax,ay,az,bx,by,bz,rasq,rbsq,rgsq,rg,rginv,ra2inv,rb2inv,rabinv; double ax, ay, az, bx, by, bz, rasq, rbsq, rgsq, rg, rginv, ra2inv, rb2inv, rabinv;
double df,df1,ddf1,fg,hg,fga,hgb,gaa,gbb; double df, df1, ddf1, fg, hg, fga, hgb, gaa, gbb;
double dtfx,dtfy,dtfz,dtgx,dtgy,dtgz,dthx,dthy,dthz; double dtfx, dtfy, dtfz, dtgx, dtgy, dtgz, dthx, dthy, dthz;
double c,s,p,sx2,sy2,sz2; double c, s, p, sx2, sy2, sz2;
int itype,jtype; int itype, jtype;
double delx,dely,delz,rsq,r2inv,r6inv; double delx, dely, delz, rsq, r2inv, r6inv;
double forcecoul,forcelj,fpair,ecoul,evdwl; double forcecoul, forcelj, fpair, ecoul, evdwl;
edihedral = evdwl = ecoul = 0.0; edihedral = evdwl = ecoul = 0.0;
ev_init(eflag,vflag); ev_init(eflag, vflag);
// insure pair->ev_tally() will use 1-4 virial contribution // insure pair->ev_tally() will use 1-4 virial contribution
@ -121,31 +120,30 @@ void DihedralCharmm::compute(int eflag, int vflag)
vb3y = x[i4][1] - x[i3][1]; vb3y = x[i4][1] - x[i3][1];
vb3z = x[i4][2] - x[i3][2]; vb3z = x[i4][2] - x[i3][2];
ax = vb1y*vb2zm - vb1z*vb2ym; ax = vb1y * vb2zm - vb1z * vb2ym;
ay = vb1z*vb2xm - vb1x*vb2zm; ay = vb1z * vb2xm - vb1x * vb2zm;
az = vb1x*vb2ym - vb1y*vb2xm; az = vb1x * vb2ym - vb1y * vb2xm;
bx = vb3y*vb2zm - vb3z*vb2ym; bx = vb3y * vb2zm - vb3z * vb2ym;
by = vb3z*vb2xm - vb3x*vb2zm; by = vb3z * vb2xm - vb3x * vb2zm;
bz = vb3x*vb2ym - vb3y*vb2xm; bz = vb3x * vb2ym - vb3y * vb2xm;
rasq = ax*ax + ay*ay + az*az; rasq = ax * ax + ay * ay + az * az;
rbsq = bx*bx + by*by + bz*bz; rbsq = bx * bx + by * by + bz * bz;
rgsq = vb2xm*vb2xm + vb2ym*vb2ym + vb2zm*vb2zm; rgsq = vb2xm * vb2xm + vb2ym * vb2ym + vb2zm * vb2zm;
rg = sqrt(rgsq); rg = sqrt(rgsq);
rginv = ra2inv = rb2inv = 0.0; rginv = ra2inv = rb2inv = 0.0;
if (rg > 0) rginv = 1.0/rg; if (rg > 0) rginv = 1.0 / rg;
if (rasq > 0) ra2inv = 1.0/rasq; if (rasq > 0) ra2inv = 1.0 / rasq;
if (rbsq > 0) rb2inv = 1.0/rbsq; if (rbsq > 0) rb2inv = 1.0 / rbsq;
rabinv = sqrt(ra2inv*rb2inv); rabinv = sqrt(ra2inv * rb2inv);
c = (ax*bx + ay*by + az*bz)*rabinv; c = (ax * bx + ay * by + az * bz) * rabinv;
s = rg*rabinv*(ax*vb3x + ay*vb3y + az*vb3z); s = rg * rabinv * (ax * vb3x + ay * vb3y + az * vb3z);
// error check // error check
if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE)) if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE)) problem(FLERR, i1, i2, i3, i4);
problem(FLERR, i1, i2, i3, i4);
if (c > 1.0) c = 1.0; if (c > 1.0) c = 1.0;
if (c < -1.0) c = -1.0; if (c < -1.0) c = -1.0;
@ -155,13 +153,13 @@ void DihedralCharmm::compute(int eflag, int vflag)
ddf1 = df1 = 0.0; ddf1 = df1 = 0.0;
for (i = 0; i < m; i++) { for (i = 0; i < m; i++) {
ddf1 = p*c - df1*s; ddf1 = p * c - df1 * s;
df1 = p*s + df1*c; df1 = p * s + df1 * c;
p = ddf1; p = ddf1;
} }
p = p*cos_shift[type] + df1*sin_shift[type]; p = p * cos_shift[type] + df1 * sin_shift[type];
df1 = df1*cos_shift[type] - ddf1*sin_shift[type]; df1 = df1 * cos_shift[type] - ddf1 * sin_shift[type];
df1 *= -m; df1 *= -m;
p += 1.0; p += 1.0;
@ -172,40 +170,40 @@ void DihedralCharmm::compute(int eflag, int vflag)
if (eflag) edihedral = k[type] * p; if (eflag) edihedral = k[type] * p;
fg = vb1x*vb2xm + vb1y*vb2ym + vb1z*vb2zm; fg = vb1x * vb2xm + vb1y * vb2ym + vb1z * vb2zm;
hg = vb3x*vb2xm + vb3y*vb2ym + vb3z*vb2zm; hg = vb3x * vb2xm + vb3y * vb2ym + vb3z * vb2zm;
fga = fg*ra2inv*rginv; fga = fg * ra2inv * rginv;
hgb = hg*rb2inv*rginv; hgb = hg * rb2inv * rginv;
gaa = -ra2inv*rg; gaa = -ra2inv * rg;
gbb = rb2inv*rg; gbb = rb2inv * rg;
dtfx = gaa*ax; dtfx = gaa * ax;
dtfy = gaa*ay; dtfy = gaa * ay;
dtfz = gaa*az; dtfz = gaa * az;
dtgx = fga*ax - hgb*bx; dtgx = fga * ax - hgb * bx;
dtgy = fga*ay - hgb*by; dtgy = fga * ay - hgb * by;
dtgz = fga*az - hgb*bz; dtgz = fga * az - hgb * bz;
dthx = gbb*bx; dthx = gbb * bx;
dthy = gbb*by; dthy = gbb * by;
dthz = gbb*bz; dthz = gbb * bz;
df = -k[type] * df1; df = -k[type] * df1;
sx2 = df*dtgx; sx2 = df * dtgx;
sy2 = df*dtgy; sy2 = df * dtgy;
sz2 = df*dtgz; sz2 = df * dtgz;
f1[0] = df*dtfx; f1[0] = df * dtfx;
f1[1] = df*dtfy; f1[1] = df * dtfy;
f1[2] = df*dtfz; f1[2] = df * dtfz;
f2[0] = sx2 - f1[0]; f2[0] = sx2 - f1[0];
f2[1] = sy2 - f1[1]; f2[1] = sy2 - f1[1];
f2[2] = sz2 - f1[2]; f2[2] = sz2 - f1[2];
f4[0] = df*dthx; f4[0] = df * dthx;
f4[1] = df*dthy; f4[1] = df * dthy;
f4[2] = df*dthz; f4[2] = df * dthz;
f3[0] = -sx2 - f4[0]; f3[0] = -sx2 - f4[0];
f3[1] = -sy2 - f4[1]; f3[1] = -sy2 - f4[1];
@ -238,8 +236,8 @@ void DihedralCharmm::compute(int eflag, int vflag)
} }
if (evflag) if (evflag)
ev_tally(i1,i2,i3,i4,nlocal,newton_bond,edihedral,f1,f3,f4, ev_tally(i1, i2, i3, i4, nlocal, newton_bond, edihedral, f1, f3, f4, vb1x, vb1y, vb1z, vb2x,
vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z); vb2y, vb2z, vb3x, vb3y, vb3z);
// 1-4 LJ and Coulomb interactions // 1-4 LJ and Coulomb interactions
// tally energy/virial in pair, using newton_bond as newton flag // tally energy/virial in pair, using newton_bond as newton flag
@ -251,34 +249,36 @@ void DihedralCharmm::compute(int eflag, int vflag)
delx = x[i1][0] - x[i4][0]; delx = x[i1][0] - x[i4][0];
dely = x[i1][1] - x[i4][1]; dely = x[i1][1] - x[i4][1];
delz = x[i1][2] - x[i4][2]; delz = x[i1][2] - x[i4][2];
rsq = delx*delx + dely*dely + delz*delz; rsq = delx * delx + dely * dely + delz * delz;
r2inv = 1.0/rsq; r2inv = 1.0 / rsq;
r6inv = r2inv*r2inv*r2inv; r6inv = r2inv * r2inv * r2inv;
if (implicit) forcecoul = qqrd2e * q[i1]*q[i4]*r2inv; if (implicit)
else forcecoul = qqrd2e * q[i1]*q[i4]*sqrt(r2inv); forcecoul = qqrd2e * q[i1] * q[i4] * r2inv;
forcelj = r6inv * (lj14_1[itype][jtype]*r6inv - lj14_2[itype][jtype]); else
fpair = weight[type] * (forcelj+forcecoul)*r2inv; forcecoul = qqrd2e * q[i1] * q[i4] * sqrt(r2inv);
forcelj = r6inv * (lj14_1[itype][jtype] * r6inv - lj14_2[itype][jtype]);
fpair = weight[type] * (forcelj + forcecoul) * r2inv;
if (eflag) { if (eflag) {
ecoul = weight[type] * forcecoul; ecoul = weight[type] * forcecoul;
evdwl = r6inv * (lj14_3[itype][jtype]*r6inv - lj14_4[itype][jtype]); evdwl = r6inv * (lj14_3[itype][jtype] * r6inv - lj14_4[itype][jtype]);
evdwl *= weight[type]; evdwl *= weight[type];
} }
if (newton_bond || i1 < nlocal) { if (newton_bond || i1 < nlocal) {
f[i1][0] += delx*fpair; f[i1][0] += delx * fpair;
f[i1][1] += dely*fpair; f[i1][1] += dely * fpair;
f[i1][2] += delz*fpair; f[i1][2] += delz * fpair;
} }
if (newton_bond || i4 < nlocal) { if (newton_bond || i4 < nlocal) {
f[i4][0] -= delx*fpair; f[i4][0] -= delx * fpair;
f[i4][1] -= dely*fpair; f[i4][1] -= dely * fpair;
f[i4][2] -= delz*fpair; f[i4][2] -= delz * fpair;
} }
if (evflag) force->pair->ev_tally(i1,i4,nlocal,newton_bond, if (evflag)
evdwl,ecoul,fpair,delx,dely,delz); force->pair->ev_tally(i1, i4, nlocal, newton_bond, evdwl, ecoul, fpair, delx, dely, delz);
} }
} }
} }
@ -288,17 +288,17 @@ void DihedralCharmm::compute(int eflag, int vflag)
void DihedralCharmm::allocate() void DihedralCharmm::allocate()
{ {
allocated = 1; allocated = 1;
int n = atom->ndihedraltypes; const int np1 = atom->ndihedraltypes + 1;
memory->create(k,n+1,"dihedral:k"); memory->create(k, np1, "dihedral:k");
memory->create(multiplicity,n+1,"dihedral:multiplicity"); memory->create(multiplicity, np1, "dihedral:multiplicity");
memory->create(shift,n+1,"dihedral:shift"); memory->create(shift, np1, "dihedral:shift");
memory->create(cos_shift,n+1,"dihedral:cos_shift"); memory->create(cos_shift, np1, "dihedral:cos_shift");
memory->create(sin_shift,n+1,"dihedral:sin_shift"); memory->create(sin_shift, np1, "dihedral:sin_shift");
memory->create(weight,n+1,"dihedral:weight"); memory->create(weight, np1, "dihedral:weight");
memory->create(setflag,n+1,"dihedral:setflag"); memory->create(setflag, np1, "dihedral:setflag");
for (int i = 1; i <= n; i++) setflag[i] = 0; for (int i = 1; i < np1; i++) setflag[i] = 0;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -307,40 +307,40 @@ void DihedralCharmm::allocate()
void DihedralCharmm::coeff(int narg, char **arg) void DihedralCharmm::coeff(int narg, char **arg)
{ {
if (narg != 5) error->all(FLERR,"Incorrect args for dihedral coefficients"); if (narg != 5) error->all(FLERR, "Incorrect args for dihedral coefficients");
if (!allocated) allocate(); if (!allocated) allocate();
int ilo,ihi; int ilo, ihi;
utils::bounds(FLERR,arg[0],1,atom->ndihedraltypes,ilo,ihi,error); utils::bounds(FLERR, arg[0], 1, atom->ndihedraltypes, ilo, ihi, error);
// require integer values of shift for backwards compatibility // require integer values of shift for backwards compatibility
// arbitrary phase angle shift could be allowed, but would break // arbitrary phase angle shift could be allowed, but would break
// backwards compatibility and is probably not needed // backwards compatibility and is probably not needed
double k_one = utils::numeric(FLERR,arg[1],false,lmp); double k_one = utils::numeric(FLERR, arg[1], false, lmp);
int multiplicity_one = utils::inumeric(FLERR,arg[2],false,lmp); int multiplicity_one = utils::inumeric(FLERR, arg[2], false, lmp);
int shift_one = utils::inumeric(FLERR,arg[3],false,lmp); int shift_one = utils::inumeric(FLERR, arg[3], false, lmp);
double weight_one = utils::numeric(FLERR,arg[4],false,lmp); double weight_one = utils::numeric(FLERR, arg[4], false, lmp);
if (multiplicity_one < 0) if (multiplicity_one < 0)
error->all(FLERR,"Incorrect multiplicity arg for dihedral coefficients"); error->all(FLERR, "Incorrect multiplicity arg for dihedral coefficients");
if (weight_one < 0.0 || weight_one > 1.0) if (weight_one < 0.0 || weight_one > 1.0)
error->all(FLERR,"Incorrect weight arg for dihedral coefficients"); error->all(FLERR, "Incorrect weight arg for dihedral coefficients");
if (weight_one > 0.0) weightflag=1; if (weight_one > 0.0) weightflag = 1;
int count = 0; int count = 0;
for (int i = ilo; i <= ihi; i++) { for (int i = ilo; i <= ihi; i++) {
k[i] = k_one; k[i] = k_one;
shift[i] = shift_one; shift[i] = shift_one;
cos_shift[i] = cos(MY_PI*shift_one/180.0); cos_shift[i] = cos(DEG2RAD * shift_one);
sin_shift[i] = sin(MY_PI*shift_one/180.0); sin_shift[i] = sin(DEG2RAD * shift_one);
multiplicity[i] = multiplicity_one; multiplicity[i] = multiplicity_one;
weight[i] = weight_one; weight[i] = weight_one;
setflag[i] = 1; setflag[i] = 1;
count++; count++;
} }
if (count == 0) error->all(FLERR,"Incorrect args for dihedral coefficients"); if (count == 0) error->all(FLERR, "Incorrect args for dihedral coefficients");
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -349,14 +349,12 @@ void DihedralCharmm::coeff(int narg, char **arg)
void DihedralCharmm::init_style() void DihedralCharmm::init_style()
{ {
if (utils::strmatch(update->integrate_style,"^respa")) { if (utils::strmatch(update->integrate_style, "^respa")) {
Respa *r = (Respa *) update->integrate; Respa *r = (Respa *) update->integrate;
if (r->level_pair >= 0 && (r->level_pair != r->level_dihedral)) if (r->level_pair >= 0 && (r->level_pair != r->level_dihedral))
error->all(FLERR,"Dihedral style charmm must be set to same" error->all(FLERR, "Dihedral style charmm must be set to same r-RESPA level as 'pair'");
" r-RESPA level as 'pair'");
if (r->level_outer >= 0 && (r->level_outer != r->level_dihedral)) if (r->level_outer >= 0 && (r->level_outer != r->level_dihedral))
error->all(FLERR,"Dihedral style charmm must be set to same" error->all(FLERR, "Dihedral style charmm must be set to same r-RESPA level as 'outer'");
" r-RESPA level as 'outer'");
} }
// insure use of CHARMM pair_style if any weight factors are non-zero // insure use of CHARMM pair_style if any weight factors are non-zero
@ -366,19 +364,20 @@ void DihedralCharmm::init_style()
if (weightflag) { if (weightflag) {
if ((force->special_lj[3] != 0.0) || (force->special_coul[3] != 0.0)) if ((force->special_lj[3] != 0.0) || (force->special_coul[3] != 0.0))
error->all(FLERR,"Must use 'special_bonds charmm' with" error->all(FLERR,
" dihedral style charmm for use with CHARMM pair styles"); "Must use 'special_bonds charmm' with dihedral "
"style charmm for use with CHARMM pair styles");
int itmp; int itmp;
if (force->pair == nullptr) if (force->pair == nullptr)
error->all(FLERR,"Dihedral charmm is incompatible with Pair style"); error->all(FLERR, "Dihedral charmm is incompatible with Pair style");
lj14_1 = (double **) force->pair->extract("lj14_1",itmp); lj14_1 = (double **) force->pair->extract("lj14_1", itmp);
lj14_2 = (double **) force->pair->extract("lj14_2",itmp); lj14_2 = (double **) force->pair->extract("lj14_2", itmp);
lj14_3 = (double **) force->pair->extract("lj14_3",itmp); lj14_3 = (double **) force->pair->extract("lj14_3", itmp);
lj14_4 = (double **) force->pair->extract("lj14_4",itmp); lj14_4 = (double **) force->pair->extract("lj14_4", itmp);
int *ptr = (int *) force->pair->extract("implicit",itmp); int *ptr = (int *) force->pair->extract("implicit", itmp);
if (!lj14_1 || !lj14_2 || !lj14_3 || !lj14_4 || !ptr) if (!lj14_1 || !lj14_2 || !lj14_3 || !lj14_4 || !ptr)
error->all(FLERR,"Dihedral charmm is incompatible with Pair style"); error->all(FLERR, "Dihedral charmm is incompatible with Pair style");
implicit = *ptr; implicit = *ptr;
} }
} }
@ -389,11 +388,11 @@ void DihedralCharmm::init_style()
void DihedralCharmm::write_restart(FILE *fp) void DihedralCharmm::write_restart(FILE *fp)
{ {
fwrite(&k[1],sizeof(double),atom->ndihedraltypes,fp); fwrite(&k[1], sizeof(double), atom->ndihedraltypes, fp);
fwrite(&multiplicity[1],sizeof(int),atom->ndihedraltypes,fp); fwrite(&multiplicity[1], sizeof(int), atom->ndihedraltypes, fp);
fwrite(&shift[1],sizeof(int),atom->ndihedraltypes,fp); fwrite(&shift[1], sizeof(int), atom->ndihedraltypes, fp);
fwrite(&weight[1],sizeof(double),atom->ndihedraltypes,fp); fwrite(&weight[1], sizeof(double), atom->ndihedraltypes, fp);
fwrite(&weightflag,sizeof(int),1,fp); fwrite(&weightflag, sizeof(int), 1, fp);
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -405,22 +404,22 @@ void DihedralCharmm::read_restart(FILE *fp)
allocate(); allocate();
if (comm->me == 0) { if (comm->me == 0) {
utils::sfread(FLERR,&k[1],sizeof(double),atom->ndihedraltypes,fp,nullptr,error); utils::sfread(FLERR, &k[1], sizeof(double), atom->ndihedraltypes, fp, nullptr, error);
utils::sfread(FLERR,&multiplicity[1],sizeof(int),atom->ndihedraltypes,fp,nullptr,error); utils::sfread(FLERR, &multiplicity[1], sizeof(int), atom->ndihedraltypes, fp, nullptr, error);
utils::sfread(FLERR,&shift[1],sizeof(int),atom->ndihedraltypes,fp,nullptr,error); utils::sfread(FLERR, &shift[1], sizeof(int), atom->ndihedraltypes, fp, nullptr, error);
utils::sfread(FLERR,&weight[1],sizeof(double),atom->ndihedraltypes,fp,nullptr,error); utils::sfread(FLERR, &weight[1], sizeof(double), atom->ndihedraltypes, fp, nullptr, error);
utils::sfread(FLERR,&weightflag,sizeof(int),1,fp,nullptr,error); utils::sfread(FLERR, &weightflag, sizeof(int), 1, fp, nullptr, error);
} }
MPI_Bcast(&k[1],atom->ndihedraltypes,MPI_DOUBLE,0,world); MPI_Bcast(&k[1], atom->ndihedraltypes, MPI_DOUBLE, 0, world);
MPI_Bcast(&multiplicity[1],atom->ndihedraltypes,MPI_INT,0,world); MPI_Bcast(&multiplicity[1], atom->ndihedraltypes, MPI_INT, 0, world);
MPI_Bcast(&shift[1],atom->ndihedraltypes,MPI_INT,0,world); MPI_Bcast(&shift[1], atom->ndihedraltypes, MPI_INT, 0, world);
MPI_Bcast(&weight[1],atom->ndihedraltypes,MPI_DOUBLE,0,world); MPI_Bcast(&weight[1], atom->ndihedraltypes, MPI_DOUBLE, 0, world);
MPI_Bcast(&weightflag,1,MPI_INT,0,world); MPI_Bcast(&weightflag, 1, MPI_INT, 0, world);
for (int i = 1; i <= atom->ndihedraltypes; i++) { for (int i = 1; i <= atom->ndihedraltypes; i++) {
setflag[i] = 1; setflag[i] = 1;
cos_shift[i] = cos(MY_PI*shift[i]/180.0); cos_shift[i] = cos(DEG2RAD * shift[i]);
sin_shift[i] = sin(MY_PI*shift[i]/180.0); sin_shift[i] = sin(DEG2RAD * shift[i]);
} }
} }
@ -431,6 +430,5 @@ void DihedralCharmm::read_restart(FILE *fp)
void DihedralCharmm::write_data(FILE *fp) void DihedralCharmm::write_data(FILE *fp)
{ {
for (int i = 1; i <= atom->ndihedraltypes; i++) for (int i = 1; i <= atom->ndihedraltypes; i++)
fprintf(fp,"%d %g %d %d %g\n",i,k[i],multiplicity[i],shift[i],weight[i]); fprintf(fp, "%d %g %d %d %g\n", i, k[i], multiplicity[i], shift[i], weight[i]);
} }

304
src/MOLECULE/dihedral_charmmfsw.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories https://www.lammps.org/, Sandia National Laboratories
@ -35,13 +34,13 @@
#include <cmath> #include <cmath>
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst; using MathConst::DEG2RAD;
#define TOLERANCE 0.05 static constexpr double TOLERANCE = 0.05;
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
DihedralCharmmfsw::DihedralCharmmfsw(LAMMPS *lmp) : Dihedral(lmp) DihedralCharmmfsw::DihedralCharmmfsw(LAMMPS *_lmp) : Dihedral(_lmp)
{ {
weightflag = 0; weightflag = 0;
writedata = 1; writedata = 1;
@ -66,19 +65,19 @@ DihedralCharmmfsw::~DihedralCharmmfsw()
void DihedralCharmmfsw::compute(int eflag, int vflag) void DihedralCharmmfsw::compute(int eflag, int vflag)
{ {
int i1,i2,i3,i4,i,m,n,type; int i1, i2, i3, i4, i, m, n, type;
double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,vb2xm,vb2ym,vb2zm; double vb1x, vb1y, vb1z, vb2x, vb2y, vb2z, vb3x, vb3y, vb3z, vb2xm, vb2ym, vb2zm;
double edihedral,f1[3],f2[3],f3[3],f4[3]; double edihedral, f1[3], f2[3], f3[3], f4[3];
double ax,ay,az,bx,by,bz,rasq,rbsq,rgsq,rg,rginv,ra2inv,rb2inv,rabinv; double ax, ay, az, bx, by, bz, rasq, rbsq, rgsq, rg, rginv, ra2inv, rb2inv, rabinv;
double df,df1,ddf1,fg,hg,fga,hgb,gaa,gbb; double df, df1, ddf1, fg, hg, fga, hgb, gaa, gbb;
double dtfx,dtfy,dtfz,dtgx,dtgy,dtgz,dthx,dthy,dthz; double dtfx, dtfy, dtfz, dtgx, dtgy, dtgz, dthx, dthy, dthz;
double c,s,p,sx2,sy2,sz2; double c, s, p, sx2, sy2, sz2;
int itype,jtype; int itype, jtype;
double delx,dely,delz,rsq,r2inv,r6inv,r; double delx, dely, delz, rsq, r2inv, r6inv, r;
double forcecoul,forcelj,fpair,ecoul,evdwl; double forcecoul, forcelj, fpair, ecoul, evdwl;
edihedral = evdwl = ecoul = 0.0; edihedral = evdwl = ecoul = 0.0;
ev_init(eflag,vflag); ev_init(eflag, vflag);
// insure pair->ev_tally() will use 1-4 virial contribution // insure pair->ev_tally() will use 1-4 virial contribution
@ -124,31 +123,30 @@ void DihedralCharmmfsw::compute(int eflag, int vflag)
vb3y = x[i4][1] - x[i3][1]; vb3y = x[i4][1] - x[i3][1];
vb3z = x[i4][2] - x[i3][2]; vb3z = x[i4][2] - x[i3][2];
ax = vb1y*vb2zm - vb1z*vb2ym; ax = vb1y * vb2zm - vb1z * vb2ym;
ay = vb1z*vb2xm - vb1x*vb2zm; ay = vb1z * vb2xm - vb1x * vb2zm;
az = vb1x*vb2ym - vb1y*vb2xm; az = vb1x * vb2ym - vb1y * vb2xm;
bx = vb3y*vb2zm - vb3z*vb2ym; bx = vb3y * vb2zm - vb3z * vb2ym;
by = vb3z*vb2xm - vb3x*vb2zm; by = vb3z * vb2xm - vb3x * vb2zm;
bz = vb3x*vb2ym - vb3y*vb2xm; bz = vb3x * vb2ym - vb3y * vb2xm;
rasq = ax*ax + ay*ay + az*az; rasq = ax * ax + ay * ay + az * az;
rbsq = bx*bx + by*by + bz*bz; rbsq = bx * bx + by * by + bz * bz;
rgsq = vb2xm*vb2xm + vb2ym*vb2ym + vb2zm*vb2zm; rgsq = vb2xm * vb2xm + vb2ym * vb2ym + vb2zm * vb2zm;
rg = sqrt(rgsq); rg = sqrt(rgsq);
rginv = ra2inv = rb2inv = 0.0; rginv = ra2inv = rb2inv = 0.0;
if (rg > 0) rginv = 1.0/rg; if (rg > 0) rginv = 1.0 / rg;
if (rasq > 0) ra2inv = 1.0/rasq; if (rasq > 0) ra2inv = 1.0 / rasq;
if (rbsq > 0) rb2inv = 1.0/rbsq; if (rbsq > 0) rb2inv = 1.0 / rbsq;
rabinv = sqrt(ra2inv*rb2inv); rabinv = sqrt(ra2inv * rb2inv);
c = (ax*bx + ay*by + az*bz)*rabinv; c = (ax * bx + ay * by + az * bz) * rabinv;
s = rg*rabinv*(ax*vb3x + ay*vb3y + az*vb3z); s = rg * rabinv * (ax * vb3x + ay * vb3y + az * vb3z);
// error check // error check
if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE)) if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE)) problem(FLERR, i1, i2, i3, i4);
problem(FLERR, i1, i2, i3, i4);
if (c > 1.0) c = 1.0; if (c > 1.0) c = 1.0;
if (c < -1.0) c = -1.0; if (c < -1.0) c = -1.0;
@ -158,13 +156,13 @@ void DihedralCharmmfsw::compute(int eflag, int vflag)
ddf1 = df1 = 0.0; ddf1 = df1 = 0.0;
for (i = 0; i < m; i++) { for (i = 0; i < m; i++) {
ddf1 = p*c - df1*s; ddf1 = p * c - df1 * s;
df1 = p*s + df1*c; df1 = p * s + df1 * c;
p = ddf1; p = ddf1;
} }
p = p*cos_shift[type] + df1*sin_shift[type]; p = p * cos_shift[type] + df1 * sin_shift[type];
df1 = df1*cos_shift[type] - ddf1*sin_shift[type]; df1 = df1 * cos_shift[type] - ddf1 * sin_shift[type];
df1 *= -m; df1 *= -m;
p += 1.0; p += 1.0;
@ -175,40 +173,40 @@ void DihedralCharmmfsw::compute(int eflag, int vflag)
if (eflag) edihedral = k[type] * p; if (eflag) edihedral = k[type] * p;
fg = vb1x*vb2xm + vb1y*vb2ym + vb1z*vb2zm; fg = vb1x * vb2xm + vb1y * vb2ym + vb1z * vb2zm;
hg = vb3x*vb2xm + vb3y*vb2ym + vb3z*vb2zm; hg = vb3x * vb2xm + vb3y * vb2ym + vb3z * vb2zm;
fga = fg*ra2inv*rginv; fga = fg * ra2inv * rginv;
hgb = hg*rb2inv*rginv; hgb = hg * rb2inv * rginv;
gaa = -ra2inv*rg; gaa = -ra2inv * rg;
gbb = rb2inv*rg; gbb = rb2inv * rg;
dtfx = gaa*ax; dtfx = gaa * ax;
dtfy = gaa*ay; dtfy = gaa * ay;
dtfz = gaa*az; dtfz = gaa * az;
dtgx = fga*ax - hgb*bx; dtgx = fga * ax - hgb * bx;
dtgy = fga*ay - hgb*by; dtgy = fga * ay - hgb * by;
dtgz = fga*az - hgb*bz; dtgz = fga * az - hgb * bz;
dthx = gbb*bx; dthx = gbb * bx;
dthy = gbb*by; dthy = gbb * by;
dthz = gbb*bz; dthz = gbb * bz;
df = -k[type] * df1; df = -k[type] * df1;
sx2 = df*dtgx; sx2 = df * dtgx;
sy2 = df*dtgy; sy2 = df * dtgy;
sz2 = df*dtgz; sz2 = df * dtgz;
f1[0] = df*dtfx; f1[0] = df * dtfx;
f1[1] = df*dtfy; f1[1] = df * dtfy;
f1[2] = df*dtfz; f1[2] = df * dtfz;
f2[0] = sx2 - f1[0]; f2[0] = sx2 - f1[0];
f2[1] = sy2 - f1[1]; f2[1] = sy2 - f1[1];
f2[2] = sz2 - f1[2]; f2[2] = sz2 - f1[2];
f4[0] = df*dthx; f4[0] = df * dthx;
f4[1] = df*dthy; f4[1] = df * dthy;
f4[2] = df*dthz; f4[2] = df * dthz;
f3[0] = -sx2 - f4[0]; f3[0] = -sx2 - f4[0];
f3[1] = -sy2 - f4[1]; f3[1] = -sy2 - f4[1];
@ -241,8 +239,8 @@ void DihedralCharmmfsw::compute(int eflag, int vflag)
} }
if (evflag) if (evflag)
ev_tally(i1,i2,i3,i4,nlocal,newton_bond,edihedral,f1,f3,f4, ev_tally(i1, i2, i3, i4, nlocal, newton_bond, edihedral, f1, f3, f4, vb1x, vb1y, vb1z, vb2x,
vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z); vb2y, vb2z, vb3x, vb3y, vb3z);
// 1-4 LJ and Coulomb interactions // 1-4 LJ and Coulomb interactions
// tally energy/virial in pair, using newton_bond as newton flag // tally energy/virial in pair, using newton_bond as newton flag
@ -254,9 +252,9 @@ void DihedralCharmmfsw::compute(int eflag, int vflag)
delx = x[i1][0] - x[i4][0]; delx = x[i1][0] - x[i4][0];
dely = x[i1][1] - x[i4][1]; dely = x[i1][1] - x[i4][1];
delz = x[i1][2] - x[i4][2]; delz = x[i1][2] - x[i4][2];
rsq = delx*delx + dely*dely + delz*delz; rsq = delx * delx + dely * dely + delz * delz;
r2inv = 1.0/rsq; r2inv = 1.0 / rsq;
r6inv = r2inv*r2inv*r2inv; r6inv = r2inv * r2inv * r2inv;
// modifying coul and LJ force and energies to apply // modifying coul and LJ force and energies to apply
// force_shift and force_switch as in CHARMM pairwise // force_shift and force_switch as in CHARMM pairwise
@ -264,39 +262,42 @@ void DihedralCharmmfsw::compute(int eflag, int vflag)
// for r < cut_inner, so switching not applied // for r < cut_inner, so switching not applied
r = sqrt(rsq); r = sqrt(rsq);
if (implicit) forcecoul = qqrd2e * q[i1]*q[i4]*r2inv; if (implicit)
else if (dihedflag) forcecoul = qqrd2e * q[i1]*q[i4]*sqrt(r2inv); forcecoul = qqrd2e * q[i1] * q[i4] * r2inv;
else forcecoul = qqrd2e * q[i1]*q[i4]*(sqrt(r2inv) - else if (dihedflag)
r*cut_coulinv14*cut_coulinv14); forcecoul = qqrd2e * q[i1] * q[i4] * sqrt(r2inv);
forcelj = r6inv * (lj14_1[itype][jtype]*r6inv - lj14_2[itype][jtype]); else
fpair = weight[type] * (forcelj+forcecoul)*r2inv; forcecoul = qqrd2e * q[i1] * q[i4] * (sqrt(r2inv) - r * cut_coulinv14 * cut_coulinv14);
forcelj = r6inv * (lj14_1[itype][jtype] * r6inv - lj14_2[itype][jtype]);
fpair = weight[type] * (forcelj + forcecoul) * r2inv;
if (eflag) { if (eflag) {
if (dihedflag) ecoul = weight[type] * forcecoul; if (dihedflag)
else ecoul = weight[type] * qqrd2e * q[i1]*q[i4] * ecoul = weight[type] * forcecoul;
(sqrt(r2inv) + r*cut_coulinv14*cut_coulinv14 - else
2.0*cut_coulinv14); ecoul = weight[type] * qqrd2e * q[i1] * q[i4] *
evdwl14_12 = r6inv*lj14_3[itype][jtype]*r6inv - (sqrt(r2inv) + r * cut_coulinv14 * cut_coulinv14 - 2.0 * cut_coulinv14);
lj14_3[itype][jtype]*cut_lj_inner6inv*cut_lj6inv; evdwl14_12 = r6inv * lj14_3[itype][jtype] * r6inv -
evdwl14_6 = -lj14_4[itype][jtype]*r6inv + lj14_3[itype][jtype] * cut_lj_inner6inv * cut_lj6inv;
lj14_4[itype][jtype]*cut_lj_inner3inv*cut_lj3inv; evdwl14_6 =
-lj14_4[itype][jtype] * r6inv + lj14_4[itype][jtype] * cut_lj_inner3inv * cut_lj3inv;
evdwl = evdwl14_12 + evdwl14_6; evdwl = evdwl14_12 + evdwl14_6;
evdwl *= weight[type]; evdwl *= weight[type];
} }
if (newton_bond || i1 < nlocal) { if (newton_bond || i1 < nlocal) {
f[i1][0] += delx*fpair; f[i1][0] += delx * fpair;
f[i1][1] += dely*fpair; f[i1][1] += dely * fpair;
f[i1][2] += delz*fpair; f[i1][2] += delz * fpair;
} }
if (newton_bond || i4 < nlocal) { if (newton_bond || i4 < nlocal) {
f[i4][0] -= delx*fpair; f[i4][0] -= delx * fpair;
f[i4][1] -= dely*fpair; f[i4][1] -= dely * fpair;
f[i4][2] -= delz*fpair; f[i4][2] -= delz * fpair;
} }
if (evflag) force->pair->ev_tally(i1,i4,nlocal,newton_bond, if (evflag)
evdwl,ecoul,fpair,delx,dely,delz); force->pair->ev_tally(i1, i4, nlocal, newton_bond, evdwl, ecoul, fpair, delx, dely, delz);
} }
} }
} }
@ -306,17 +307,17 @@ void DihedralCharmmfsw::compute(int eflag, int vflag)
void DihedralCharmmfsw::allocate() void DihedralCharmmfsw::allocate()
{ {
allocated = 1; allocated = 1;
int n = atom->ndihedraltypes; const int np1 = atom->ndihedraltypes + 1;
memory->create(k,n+1,"dihedral:k"); memory->create(k, np1, "dihedral:k");
memory->create(multiplicity,n+1,"dihedral:multiplicity"); memory->create(multiplicity, np1, "dihedral:multiplicity");
memory->create(shift,n+1,"dihedral:shift"); memory->create(shift, np1, "dihedral:shift");
memory->create(cos_shift,n+1,"dihedral:cos_shift"); memory->create(cos_shift, np1, "dihedral:cos_shift");
memory->create(sin_shift,n+1,"dihedral:sin_shift"); memory->create(sin_shift, np1, "dihedral:sin_shift");
memory->create(weight,n+1,"dihedral:weight"); memory->create(weight, np1, "dihedral:weight");
memory->create(setflag,n+1,"dihedral:setflag"); memory->create(setflag, np1, "dihedral:setflag");
for (int i = 1; i <= n; i++) setflag[i] = 0; for (int i = 1; i < np1; i++) setflag[i] = 0;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -325,40 +326,40 @@ void DihedralCharmmfsw::allocate()
void DihedralCharmmfsw::coeff(int narg, char **arg) void DihedralCharmmfsw::coeff(int narg, char **arg)
{ {
if (narg != 5) error->all(FLERR,"Incorrect args for dihedral coefficients"); if (narg != 5) error->all(FLERR, "Incorrect args for dihedral coefficients");
if (!allocated) allocate(); if (!allocated) allocate();
int ilo,ihi; int ilo, ihi;
utils::bounds(FLERR,arg[0],1,atom->ndihedraltypes,ilo,ihi,error); utils::bounds(FLERR, arg[0], 1, atom->ndihedraltypes, ilo, ihi, error);
// require integer values of shift for backwards compatibility // require integer values of shift for backwards compatibility
// arbitrary phase angle shift could be allowed, but would break // arbitrary phase angle shift could be allowed, but would break
// backwards compatibility and is probably not needed // backwards compatibility and is probably not needed
double k_one = utils::numeric(FLERR,arg[1],false,lmp); double k_one = utils::numeric(FLERR, arg[1], false, lmp);
int multiplicity_one = utils::inumeric(FLERR,arg[2],false,lmp); int multiplicity_one = utils::inumeric(FLERR, arg[2], false, lmp);
int shift_one = utils::inumeric(FLERR,arg[3],false,lmp); int shift_one = utils::inumeric(FLERR, arg[3], false, lmp);
double weight_one = utils::numeric(FLERR,arg[4],false,lmp); double weight_one = utils::numeric(FLERR, arg[4], false, lmp);
if (multiplicity_one < 0) if (multiplicity_one < 0)
error->all(FLERR,"Incorrect multiplicity arg for dihedral coefficients"); error->all(FLERR, "Incorrect multiplicity arg for dihedral coefficients");
if (weight_one < 0.0 || weight_one > 1.0) if (weight_one < 0.0 || weight_one > 1.0)
error->all(FLERR,"Incorrect weight arg for dihedral coefficients"); error->all(FLERR, "Incorrect weight arg for dihedral coefficients");
if (weight_one > 0.0) weightflag=1; if (weight_one > 0.0) weightflag = 1;
int count = 0; int count = 0;
for (int i = ilo; i <= ihi; i++) { for (int i = ilo; i <= ihi; i++) {
k[i] = k_one; k[i] = k_one;
shift[i] = shift_one; shift[i] = shift_one;
cos_shift[i] = cos(MY_PI*shift_one/180.0); cos_shift[i] = cos(DEG2RAD * shift_one);
sin_shift[i] = sin(MY_PI*shift_one/180.0); sin_shift[i] = sin(DEG2RAD * shift_one);
multiplicity[i] = multiplicity_one; multiplicity[i] = multiplicity_one;
weight[i] = weight_one; weight[i] = weight_one;
setflag[i] = 1; setflag[i] = 1;
count++; count++;
} }
if (count == 0) error->all(FLERR,"Incorrect args for dihedral coefficients"); if (count == 0) error->all(FLERR, "Incorrect args for dihedral coefficients");
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -367,14 +368,12 @@ void DihedralCharmmfsw::coeff(int narg, char **arg)
void DihedralCharmmfsw::init_style() void DihedralCharmmfsw::init_style()
{ {
if (utils::strmatch(update->integrate_style,"^respa")) { if (utils::strmatch(update->integrate_style, "^respa")) {
Respa *r = (Respa *) update->integrate; Respa *r = (Respa *) update->integrate;
if (r->level_pair >= 0 && (r->level_pair != r->level_dihedral)) if (r->level_pair >= 0 && (r->level_pair != r->level_dihedral))
error->all(FLERR,"Dihedral style charmmfsw must be set to same" error->all(FLERR, "Dihedral style charmmfsw must be set to same r-RESPA level as 'pair'");
" r-RESPA level as 'pair'");
if (r->level_outer >= 0 && (r->level_outer != r->level_dihedral)) if (r->level_outer >= 0 && (r->level_outer != r->level_dihedral))
error->all(FLERR,"Dihedral style charmmfsw must be set to same" error->all(FLERR, "Dihedral style charmmfsw must be set to same r-RESPA level as 'outer'");
" r-RESPA level as 'outer'");
} }
// insure use of CHARMM pair_style if any weight factors are non-zero // insure use of CHARMM pair_style if any weight factors are non-zero
@ -384,19 +383,20 @@ void DihedralCharmmfsw::init_style()
if (weightflag) { if (weightflag) {
if ((force->special_lj[3] != 0.0) || (force->special_coul[3] != 0.0)) if ((force->special_lj[3] != 0.0) || (force->special_coul[3] != 0.0))
error->all(FLERR,"Must use 'special_bonds charmm' with" error->all(FLERR,
" dihedral style charmm for use with CHARMM pair styles"); "Must use 'special_bonds charmm' with dihedral "
"style charmm for use with CHARMM pair styles");
int itmp; int itmp;
if (force->pair == nullptr) if (force->pair == nullptr)
error->all(FLERR,"Dihedral charmmfsw is incompatible with Pair style"); error->all(FLERR, "Dihedral charmmfsw is incompatible with Pair style");
lj14_1 = (double **) force->pair->extract("lj14_1",itmp); lj14_1 = (double **) force->pair->extract("lj14_1", itmp);
lj14_2 = (double **) force->pair->extract("lj14_2",itmp); lj14_2 = (double **) force->pair->extract("lj14_2", itmp);
lj14_3 = (double **) force->pair->extract("lj14_3",itmp); lj14_3 = (double **) force->pair->extract("lj14_3", itmp);
lj14_4 = (double **) force->pair->extract("lj14_4",itmp); lj14_4 = (double **) force->pair->extract("lj14_4", itmp);
int *ptr = (int *) force->pair->extract("implicit",itmp); int *ptr = (int *) force->pair->extract("implicit", itmp);
if (!lj14_1 || !lj14_2 || !lj14_3 || !lj14_4 || !ptr) if (!lj14_1 || !lj14_2 || !lj14_3 || !lj14_4 || !ptr)
error->all(FLERR,"Dihedral charmmfsw is incompatible with Pair style"); error->all(FLERR, "Dihedral charmmfsw is incompatible with Pair style");
implicit = *ptr; implicit = *ptr;
} }
@ -404,25 +404,24 @@ void DihedralCharmmfsw::init_style()
// to 1/4 dihedral atoms to match CHARMM pairwise interactions // to 1/4 dihedral atoms to match CHARMM pairwise interactions
int itmp; int itmp;
int *p_dihedflag = (int *) force->pair->extract("dihedflag",itmp); int *p_dihedflag = (int *) force->pair->extract("dihedflag", itmp);
double *p_cutljinner = (double *) force->pair->extract("cut_lj_inner",itmp); double *p_cutljinner = (double *) force->pair->extract("cut_lj_inner", itmp);
double *p_cutlj = (double *) force->pair->extract("cut_lj",itmp); double *p_cutlj = (double *) force->pair->extract("cut_lj", itmp);
double *p_cutcoul = (double *) force->pair->extract("cut_coul",itmp); double *p_cutcoul = (double *) force->pair->extract("cut_coul", itmp);
if (p_cutcoul == nullptr || p_cutljinner == nullptr || if (p_cutcoul == nullptr || p_cutljinner == nullptr || p_cutlj == nullptr ||
p_cutlj == nullptr || p_dihedflag == nullptr) p_dihedflag == nullptr)
error->all(FLERR,"Dihedral charmmfsw is incompatible with Pair style"); error->all(FLERR, "Dihedral charmmfsw is incompatible with Pair style");
dihedflag = *p_dihedflag; dihedflag = *p_dihedflag;
cut_coul14 = *p_cutcoul; cut_coul14 = *p_cutcoul;
cut_lj_inner14 = *p_cutljinner; cut_lj_inner14 = *p_cutljinner;
cut_lj14 = *p_cutlj; cut_lj14 = *p_cutlj;
cut_coulinv14 = 1/cut_coul14; cut_coulinv14 = 1 / cut_coul14;
cut_lj_inner3inv = (1/cut_lj_inner14) * (1/cut_lj_inner14) * cut_lj_inner3inv = (1 / cut_lj_inner14) * (1 / cut_lj_inner14) * (1 / cut_lj_inner14);
(1/cut_lj_inner14);
cut_lj_inner6inv = cut_lj_inner3inv * cut_lj_inner3inv; cut_lj_inner6inv = cut_lj_inner3inv * cut_lj_inner3inv;
cut_lj3inv = (1/cut_lj14) * (1/cut_lj14) * (1/cut_lj14); cut_lj3inv = (1 / cut_lj14) * (1 / cut_lj14) * (1 / cut_lj14);
cut_lj6inv = cut_lj3inv * cut_lj3inv; cut_lj6inv = cut_lj3inv * cut_lj3inv;
} }
@ -432,11 +431,11 @@ void DihedralCharmmfsw::init_style()
void DihedralCharmmfsw::write_restart(FILE *fp) void DihedralCharmmfsw::write_restart(FILE *fp)
{ {
fwrite(&k[1],sizeof(double),atom->ndihedraltypes,fp); fwrite(&k[1], sizeof(double), atom->ndihedraltypes, fp);
fwrite(&multiplicity[1],sizeof(int),atom->ndihedraltypes,fp); fwrite(&multiplicity[1], sizeof(int), atom->ndihedraltypes, fp);
fwrite(&shift[1],sizeof(int),atom->ndihedraltypes,fp); fwrite(&shift[1], sizeof(int), atom->ndihedraltypes, fp);
fwrite(&weight[1],sizeof(double),atom->ndihedraltypes,fp); fwrite(&weight[1], sizeof(double), atom->ndihedraltypes, fp);
fwrite(&weightflag,sizeof(int),1,fp); fwrite(&weightflag, sizeof(int), 1, fp);
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -448,22 +447,22 @@ void DihedralCharmmfsw::read_restart(FILE *fp)
allocate(); allocate();
if (comm->me == 0) { if (comm->me == 0) {
utils::sfread(FLERR,&k[1],sizeof(double),atom->ndihedraltypes,fp,nullptr,error); utils::sfread(FLERR, &k[1], sizeof(double), atom->ndihedraltypes, fp, nullptr, error);
utils::sfread(FLERR,&multiplicity[1],sizeof(int),atom->ndihedraltypes,fp,nullptr,error); utils::sfread(FLERR, &multiplicity[1], sizeof(int), atom->ndihedraltypes, fp, nullptr, error);
utils::sfread(FLERR,&shift[1],sizeof(int),atom->ndihedraltypes,fp,nullptr,error); utils::sfread(FLERR, &shift[1], sizeof(int), atom->ndihedraltypes, fp, nullptr, error);
utils::sfread(FLERR,&weight[1],sizeof(double),atom->ndihedraltypes,fp,nullptr,error); utils::sfread(FLERR, &weight[1], sizeof(double), atom->ndihedraltypes, fp, nullptr, error);
utils::sfread(FLERR,&weightflag,sizeof(int),1,fp,nullptr,error); utils::sfread(FLERR, &weightflag, sizeof(int), 1, fp, nullptr, error);
} }
MPI_Bcast(&k[1],atom->ndihedraltypes,MPI_DOUBLE,0,world); MPI_Bcast(&k[1], atom->ndihedraltypes, MPI_DOUBLE, 0, world);
MPI_Bcast(&multiplicity[1],atom->ndihedraltypes,MPI_INT,0,world); MPI_Bcast(&multiplicity[1], atom->ndihedraltypes, MPI_INT, 0, world);
MPI_Bcast(&shift[1],atom->ndihedraltypes,MPI_INT,0,world); MPI_Bcast(&shift[1], atom->ndihedraltypes, MPI_INT, 0, world);
MPI_Bcast(&weight[1],atom->ndihedraltypes,MPI_DOUBLE,0,world); MPI_Bcast(&weight[1], atom->ndihedraltypes, MPI_DOUBLE, 0, world);
MPI_Bcast(&weightflag,1,MPI_INT,0,world); MPI_Bcast(&weightflag, 1, MPI_INT, 0, world);
for (int i = 1; i <= atom->ndihedraltypes; i++) { for (int i = 1; i <= atom->ndihedraltypes; i++) {
setflag[i] = 1; setflag[i] = 1;
cos_shift[i] = cos(MY_PI*shift[i]/180.0); cos_shift[i] = cos(DEG2RAD * shift[i]);
sin_shift[i] = sin(MY_PI*shift[i]/180.0); sin_shift[i] = sin(DEG2RAD * shift[i]);
} }
} }
@ -474,6 +473,5 @@ void DihedralCharmmfsw::read_restart(FILE *fp)
void DihedralCharmmfsw::write_data(FILE *fp) void DihedralCharmmfsw::write_data(FILE *fp)
{ {
for (int i = 1; i <= atom->ndihedraltypes; i++) for (int i = 1; i <= atom->ndihedraltypes; i++)
fprintf(fp,"%d %g %d %d %g\n",i,k[i],multiplicity[i],shift[i],weight[i]); fprintf(fp, "%d %g %d %d %g\n", i, k[i], multiplicity[i], shift[i], weight[i]);
} }

170
src/MOLECULE/dihedral_harmonic.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories https://www.lammps.org/, Sandia National Laboratories
@ -29,12 +28,11 @@
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
#define TOLERANCE 0.05 static constexpr double TOLERANCE = 0.05;
#define SMALL 0.001
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
DihedralHarmonic::DihedralHarmonic(LAMMPS *lmp) : Dihedral(lmp) DihedralHarmonic::DihedralHarmonic(LAMMPS *_lmp) : Dihedral(_lmp)
{ {
writedata = 1; writedata = 1;
} }
@ -57,16 +55,16 @@ DihedralHarmonic::~DihedralHarmonic()
void DihedralHarmonic::compute(int eflag, int vflag) void DihedralHarmonic::compute(int eflag, int vflag)
{ {
int i1,i2,i3,i4,i,m,n,type; int i1, i2, i3, i4, i, m, n, type;
double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,vb2xm,vb2ym,vb2zm; double vb1x, vb1y, vb1z, vb2x, vb2y, vb2z, vb3x, vb3y, vb3z, vb2xm, vb2ym, vb2zm;
double edihedral,f1[3],f2[3],f3[3],f4[3]; double edihedral, f1[3], f2[3], f3[3], f4[3];
double ax,ay,az,bx,by,bz,rasq,rbsq,rgsq,rg,rginv,ra2inv,rb2inv,rabinv; double ax, ay, az, bx, by, bz, rasq, rbsq, rgsq, rg, rginv, ra2inv, rb2inv, rabinv;
double df,df1,ddf1,fg,hg,fga,hgb,gaa,gbb; double df, df1, ddf1, fg, hg, fga, hgb, gaa, gbb;
double dtfx,dtfy,dtfz,dtgx,dtgy,dtgz,dthx,dthy,dthz; double dtfx, dtfy, dtfz, dtgx, dtgy, dtgz, dthx, dthy, dthz;
double c,s,p,sx2,sy2,sz2; double c, s, p, sx2, sy2, sz2;
edihedral = 0.0; edihedral = 0.0;
ev_init(eflag,vflag); ev_init(eflag, vflag);
double **x = atom->x; double **x = atom->x;
double **f = atom->f; double **f = atom->f;
@ -106,31 +104,30 @@ void DihedralHarmonic::compute(int eflag, int vflag)
// c,s calculation // c,s calculation
ax = vb1y*vb2zm - vb1z*vb2ym; ax = vb1y * vb2zm - vb1z * vb2ym;
ay = vb1z*vb2xm - vb1x*vb2zm; ay = vb1z * vb2xm - vb1x * vb2zm;
az = vb1x*vb2ym - vb1y*vb2xm; az = vb1x * vb2ym - vb1y * vb2xm;
bx = vb3y*vb2zm - vb3z*vb2ym; bx = vb3y * vb2zm - vb3z * vb2ym;
by = vb3z*vb2xm - vb3x*vb2zm; by = vb3z * vb2xm - vb3x * vb2zm;
bz = vb3x*vb2ym - vb3y*vb2xm; bz = vb3x * vb2ym - vb3y * vb2xm;
rasq = ax*ax + ay*ay + az*az; rasq = ax * ax + ay * ay + az * az;
rbsq = bx*bx + by*by + bz*bz; rbsq = bx * bx + by * by + bz * bz;
rgsq = vb2xm*vb2xm + vb2ym*vb2ym + vb2zm*vb2zm; rgsq = vb2xm * vb2xm + vb2ym * vb2ym + vb2zm * vb2zm;
rg = sqrt(rgsq); rg = sqrt(rgsq);
rginv = ra2inv = rb2inv = 0.0; rginv = ra2inv = rb2inv = 0.0;
if (rg > 0) rginv = 1.0/rg; if (rg > 0) rginv = 1.0 / rg;
if (rasq > 0) ra2inv = 1.0/rasq; if (rasq > 0) ra2inv = 1.0 / rasq;
if (rbsq > 0) rb2inv = 1.0/rbsq; if (rbsq > 0) rb2inv = 1.0 / rbsq;
rabinv = sqrt(ra2inv*rb2inv); rabinv = sqrt(ra2inv * rb2inv);
c = (ax*bx + ay*by + az*bz)*rabinv; c = (ax * bx + ay * by + az * bz) * rabinv;
s = rg*rabinv*(ax*vb3x + ay*vb3y + az*vb3z); s = rg * rabinv * (ax * vb3x + ay * vb3y + az * vb3z);
// error check // error check
if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE)) if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE)) problem(FLERR, i1, i2, i3, i4);
problem(FLERR, i1, i2, i3, i4);
if (c > 1.0) c = 1.0; if (c > 1.0) c = 1.0;
if (c < -1.0) c = -1.0; if (c < -1.0) c = -1.0;
@ -140,13 +137,13 @@ void DihedralHarmonic::compute(int eflag, int vflag)
ddf1 = df1 = 0.0; ddf1 = df1 = 0.0;
for (i = 0; i < m; i++) { for (i = 0; i < m; i++) {
ddf1 = p*c - df1*s; ddf1 = p * c - df1 * s;
df1 = p*s + df1*c; df1 = p * s + df1 * c;
p = ddf1; p = ddf1;
} }
p = p*cos_shift[type] + df1*sin_shift[type]; p = p * cos_shift[type] + df1 * sin_shift[type];
df1 = df1*cos_shift[type] - ddf1*sin_shift[type]; df1 = df1 * cos_shift[type] - ddf1 * sin_shift[type];
df1 *= -m; df1 *= -m;
p += 1.0; p += 1.0;
@ -157,40 +154,40 @@ void DihedralHarmonic::compute(int eflag, int vflag)
if (eflag) edihedral = k[type] * p; if (eflag) edihedral = k[type] * p;
fg = vb1x*vb2xm + vb1y*vb2ym + vb1z*vb2zm; fg = vb1x * vb2xm + vb1y * vb2ym + vb1z * vb2zm;
hg = vb3x*vb2xm + vb3y*vb2ym + vb3z*vb2zm; hg = vb3x * vb2xm + vb3y * vb2ym + vb3z * vb2zm;
fga = fg*ra2inv*rginv; fga = fg * ra2inv * rginv;
hgb = hg*rb2inv*rginv; hgb = hg * rb2inv * rginv;
gaa = -ra2inv*rg; gaa = -ra2inv * rg;
gbb = rb2inv*rg; gbb = rb2inv * rg;
dtfx = gaa*ax; dtfx = gaa * ax;
dtfy = gaa*ay; dtfy = gaa * ay;
dtfz = gaa*az; dtfz = gaa * az;
dtgx = fga*ax - hgb*bx; dtgx = fga * ax - hgb * bx;
dtgy = fga*ay - hgb*by; dtgy = fga * ay - hgb * by;
dtgz = fga*az - hgb*bz; dtgz = fga * az - hgb * bz;
dthx = gbb*bx; dthx = gbb * bx;
dthy = gbb*by; dthy = gbb * by;
dthz = gbb*bz; dthz = gbb * bz;
df = -k[type] * df1; df = -k[type] * df1;
sx2 = df*dtgx; sx2 = df * dtgx;
sy2 = df*dtgy; sy2 = df * dtgy;
sz2 = df*dtgz; sz2 = df * dtgz;
f1[0] = df*dtfx; f1[0] = df * dtfx;
f1[1] = df*dtfy; f1[1] = df * dtfy;
f1[2] = df*dtfz; f1[2] = df * dtfz;
f2[0] = sx2 - f1[0]; f2[0] = sx2 - f1[0];
f2[1] = sy2 - f1[1]; f2[1] = sy2 - f1[1];
f2[2] = sz2 - f1[2]; f2[2] = sz2 - f1[2];
f4[0] = df*dthx; f4[0] = df * dthx;
f4[1] = df*dthy; f4[1] = df * dthy;
f4[2] = df*dthz; f4[2] = df * dthz;
f3[0] = -sx2 - f4[0]; f3[0] = -sx2 - f4[0];
f3[1] = -sy2 - f4[1]; f3[1] = -sy2 - f4[1];
@ -223,8 +220,8 @@ void DihedralHarmonic::compute(int eflag, int vflag)
} }
if (evflag) if (evflag)
ev_tally(i1,i2,i3,i4,nlocal,newton_bond,edihedral,f1,f3,f4, ev_tally(i1, i2, i3, i4, nlocal, newton_bond, edihedral, f1, f3, f4, vb1x, vb1y, vb1z, vb2x,
vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z); vb2y, vb2z, vb3x, vb3y, vb3z);
} }
} }
@ -233,16 +230,16 @@ void DihedralHarmonic::compute(int eflag, int vflag)
void DihedralHarmonic::allocate() void DihedralHarmonic::allocate()
{ {
allocated = 1; allocated = 1;
int n = atom->ndihedraltypes; const int np1 = atom->ndihedraltypes + 1;
memory->create(k,n+1,"dihedral:k"); memory->create(k, np1, "dihedral:k");
memory->create(sign,n+1,"dihedral:sign"); memory->create(sign, np1, "dihedral:sign");
memory->create(multiplicity,n+1,"dihedral:multiplicity"); memory->create(multiplicity, np1, "dihedral:multiplicity");
memory->create(cos_shift,n+1,"dihedral:cos_shift"); memory->create(cos_shift, np1, "dihedral:cos_shift");
memory->create(sin_shift,n+1,"dihedral:sin_shift"); memory->create(sin_shift, np1, "dihedral:sin_shift");
memory->create(setflag,n+1,"dihedral:setflag"); memory->create(setflag, np1, "dihedral:setflag");
for (int i = 1; i <= n; i++) setflag[i] = 0; for (int i = 1; i < np1; i++) setflag[i] = 0;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -251,24 +248,24 @@ void DihedralHarmonic::allocate()
void DihedralHarmonic::coeff(int narg, char **arg) void DihedralHarmonic::coeff(int narg, char **arg)
{ {
if (narg != 4) error->all(FLERR,"Incorrect args for dihedral coefficients"); if (narg != 4) error->all(FLERR, "Incorrect args for dihedral coefficients");
if (!allocated) allocate(); if (!allocated) allocate();
int ilo,ihi; int ilo, ihi;
utils::bounds(FLERR,arg[0],1,atom->ndihedraltypes,ilo,ihi,error); utils::bounds(FLERR, arg[0], 1, atom->ndihedraltypes, ilo, ihi, error);
double k_one = utils::numeric(FLERR,arg[1],false,lmp); double k_one = utils::numeric(FLERR, arg[1], false, lmp);
int sign_one = utils::inumeric(FLERR,arg[2],false,lmp); int sign_one = utils::inumeric(FLERR, arg[2], false, lmp);
int multiplicity_one = utils::inumeric(FLERR,arg[3],false,lmp); int multiplicity_one = utils::inumeric(FLERR, arg[3], false, lmp);
// require sign = +/- 1 for backwards compatibility // require sign = +/- 1 for backwards compatibility
// arbitrary phase angle shift could be allowed, but would break // arbitrary phase angle shift could be allowed, but would break
// backwards compatibility and is probably not needed // backwards compatibility and is probably not needed
if (sign_one != -1 && sign_one != 1) if (sign_one != -1 && sign_one != 1)
error->all(FLERR,"Incorrect sign arg for dihedral coefficients"); error->all(FLERR, "Incorrect sign arg for dihedral coefficients");
if (multiplicity_one < 0) if (multiplicity_one < 0)
error->all(FLERR,"Incorrect multiplicity arg for dihedral coefficients"); error->all(FLERR, "Incorrect multiplicity arg for dihedral coefficients");
int count = 0; int count = 0;
for (int i = ilo; i <= ihi; i++) { for (int i = ilo; i <= ihi; i++) {
@ -286,7 +283,7 @@ void DihedralHarmonic::coeff(int narg, char **arg)
count++; count++;
} }
if (count == 0) error->all(FLERR,"Incorrect args for dihedral coefficients"); if (count == 0) error->all(FLERR, "Incorrect args for dihedral coefficients");
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -295,9 +292,9 @@ void DihedralHarmonic::coeff(int narg, char **arg)
void DihedralHarmonic::write_restart(FILE *fp) void DihedralHarmonic::write_restart(FILE *fp)
{ {
fwrite(&k[1],sizeof(double),atom->ndihedraltypes,fp); fwrite(&k[1], sizeof(double), atom->ndihedraltypes, fp);
fwrite(&sign[1],sizeof(int),atom->ndihedraltypes,fp); fwrite(&sign[1], sizeof(int), atom->ndihedraltypes, fp);
fwrite(&multiplicity[1],sizeof(int),atom->ndihedraltypes,fp); fwrite(&multiplicity[1], sizeof(int), atom->ndihedraltypes, fp);
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -309,13 +306,13 @@ void DihedralHarmonic::read_restart(FILE *fp)
allocate(); allocate();
if (comm->me == 0) { if (comm->me == 0) {
utils::sfread(FLERR,&k[1],sizeof(double),atom->ndihedraltypes,fp,nullptr,error); utils::sfread(FLERR, &k[1], sizeof(double), atom->ndihedraltypes, fp, nullptr, error);
utils::sfread(FLERR,&sign[1],sizeof(int),atom->ndihedraltypes,fp,nullptr,error); utils::sfread(FLERR, &sign[1], sizeof(int), atom->ndihedraltypes, fp, nullptr, error);
utils::sfread(FLERR,&multiplicity[1],sizeof(int),atom->ndihedraltypes,fp,nullptr,error); utils::sfread(FLERR, &multiplicity[1], sizeof(int), atom->ndihedraltypes, fp, nullptr, error);
} }
MPI_Bcast(&k[1],atom->ndihedraltypes,MPI_DOUBLE,0,world); MPI_Bcast(&k[1], atom->ndihedraltypes, MPI_DOUBLE, 0, world);
MPI_Bcast(&sign[1],atom->ndihedraltypes,MPI_INT,0,world); MPI_Bcast(&sign[1], atom->ndihedraltypes, MPI_INT, 0, world);
MPI_Bcast(&multiplicity[1],atom->ndihedraltypes,MPI_INT,0,world); MPI_Bcast(&multiplicity[1], atom->ndihedraltypes, MPI_INT, 0, world);
for (int i = 1; i <= atom->ndihedraltypes; i++) { for (int i = 1; i <= atom->ndihedraltypes; i++) {
setflag[i] = 1; setflag[i] = 1;
@ -336,6 +333,5 @@ void DihedralHarmonic::read_restart(FILE *fp)
void DihedralHarmonic::write_data(FILE *fp) void DihedralHarmonic::write_data(FILE *fp)
{ {
for (int i = 1; i <= atom->ndihedraltypes; i++) for (int i = 1; i <= atom->ndihedraltypes; i++)
fprintf(fp,"%d %g %d %d\n",i,k[i],sign[i],multiplicity[i]); fprintf(fp, "%d %g %d %d\n", i, k[i], sign[i], multiplicity[i]);
} }

164
src/MOLECULE/dihedral_multi_harmonic.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories https://www.lammps.org/, Sandia National Laboratories
@ -29,12 +28,12 @@
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
#define TOLERANCE 0.05 static constexpr double TOLERANCE = 0.05;
#define SMALL 0.001 static constexpr double SMALL = 0.001;
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
DihedralMultiHarmonic::DihedralMultiHarmonic(LAMMPS *lmp) : Dihedral(lmp) DihedralMultiHarmonic::DihedralMultiHarmonic(LAMMPS *_lmp) : Dihedral(_lmp)
{ {
writedata = 1; writedata = 1;
} }
@ -57,17 +56,17 @@ DihedralMultiHarmonic::~DihedralMultiHarmonic()
void DihedralMultiHarmonic::compute(int eflag, int vflag) void DihedralMultiHarmonic::compute(int eflag, int vflag)
{ {
int i1,i2,i3,i4,n,type; int i1, i2, i3, i4, n, type;
double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,vb2xm,vb2ym,vb2zm; double vb1x, vb1y, vb1z, vb2x, vb2y, vb2z, vb3x, vb3y, vb3z, vb2xm, vb2ym, vb2zm;
double edihedral,f1[3],f2[3],f3[3],f4[3]; double edihedral, f1[3], f2[3], f3[3], f4[3];
double sb1,sb2,sb3,rb1,rb3,c0,b1mag2,b1mag,b2mag2; double sb1, sb2, sb3, rb1, rb3, c0, b1mag2, b1mag, b2mag2;
double b2mag,b3mag2,b3mag,ctmp,r12c1,c1mag,r12c2; double b2mag, b3mag2, b3mag, ctmp, r12c1, c1mag, r12c2;
double c2mag,sc1,sc2,s1,s12,c,p,pd,a,a11,a22; double c2mag, sc1, sc2, s1, s12, c, p, pd, a, a11, a22;
double a33,a12,a13,a23,sx2,sy2,sz2; double a33, a12, a13, a23, sx2, sy2, sz2;
double s2,sin2; double s2, sin2;
edihedral = 0.0; edihedral = 0.0;
ev_init(eflag,vflag); ev_init(eflag, vflag);
double **x = atom->x; double **x = atom->x;
double **f = atom->f; double **f = atom->f;
@ -107,53 +106,52 @@ void DihedralMultiHarmonic::compute(int eflag, int vflag)
// c0 calculation // c0 calculation
sb1 = 1.0 / (vb1x*vb1x + vb1y*vb1y + vb1z*vb1z); sb1 = 1.0 / (vb1x * vb1x + vb1y * vb1y + vb1z * vb1z);
sb2 = 1.0 / (vb2x*vb2x + vb2y*vb2y + vb2z*vb2z); sb2 = 1.0 / (vb2x * vb2x + vb2y * vb2y + vb2z * vb2z);
sb3 = 1.0 / (vb3x*vb3x + vb3y*vb3y + vb3z*vb3z); sb3 = 1.0 / (vb3x * vb3x + vb3y * vb3y + vb3z * vb3z);
rb1 = sqrt(sb1); rb1 = sqrt(sb1);
rb3 = sqrt(sb3); rb3 = sqrt(sb3);
c0 = (vb1x*vb3x + vb1y*vb3y + vb1z*vb3z) * rb1*rb3; c0 = (vb1x * vb3x + vb1y * vb3y + vb1z * vb3z) * rb1 * rb3;
// 1st and 2nd angle // 1st and 2nd angle
b1mag2 = vb1x*vb1x + vb1y*vb1y + vb1z*vb1z; b1mag2 = vb1x * vb1x + vb1y * vb1y + vb1z * vb1z;
b1mag = sqrt(b1mag2); b1mag = sqrt(b1mag2);
b2mag2 = vb2x*vb2x + vb2y*vb2y + vb2z*vb2z; b2mag2 = vb2x * vb2x + vb2y * vb2y + vb2z * vb2z;
b2mag = sqrt(b2mag2); b2mag = sqrt(b2mag2);
b3mag2 = vb3x*vb3x + vb3y*vb3y + vb3z*vb3z; b3mag2 = vb3x * vb3x + vb3y * vb3y + vb3z * vb3z;
b3mag = sqrt(b3mag2); b3mag = sqrt(b3mag2);
ctmp = vb1x*vb2x + vb1y*vb2y + vb1z*vb2z; ctmp = vb1x * vb2x + vb1y * vb2y + vb1z * vb2z;
r12c1 = 1.0 / (b1mag*b2mag); r12c1 = 1.0 / (b1mag * b2mag);
c1mag = ctmp * r12c1; c1mag = ctmp * r12c1;
ctmp = vb2xm*vb3x + vb2ym*vb3y + vb2zm*vb3z; ctmp = vb2xm * vb3x + vb2ym * vb3y + vb2zm * vb3z;
r12c2 = 1.0 / (b2mag*b3mag); r12c2 = 1.0 / (b2mag * b3mag);
c2mag = ctmp * r12c2; c2mag = ctmp * r12c2;
// cos and sin of 2 angles and final c // cos and sin of 2 angles and final c
sin2 = MAX(1.0 - c1mag*c1mag,0.0); sin2 = MAX(1.0 - c1mag * c1mag, 0.0);
sc1 = sqrt(sin2); sc1 = sqrt(sin2);
if (sc1 < SMALL) sc1 = SMALL; if (sc1 < SMALL) sc1 = SMALL;
sc1 = 1.0/sc1; sc1 = 1.0 / sc1;
sin2 = MAX(1.0 - c2mag*c2mag,0.0); sin2 = MAX(1.0 - c2mag * c2mag, 0.0);
sc2 = sqrt(sin2); sc2 = sqrt(sin2);
if (sc2 < SMALL) sc2 = SMALL; if (sc2 < SMALL) sc2 = SMALL;
sc2 = 1.0/sc2; sc2 = 1.0 / sc2;
s1 = sc1 * sc1; s1 = sc1 * sc1;
s2 = sc2 * sc2; s2 = sc2 * sc2;
s12 = sc1 * sc2; s12 = sc1 * sc2;
c = (c0 + c1mag*c2mag) * s12; c = (c0 + c1mag * c2mag) * s12;
// error check // error check
if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE)) if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE)) problem(FLERR, i1, i2, i3, i4);
problem(FLERR, i1, i2, i3, i4);
if (c > 1.0) c = 1.0; if (c > 1.0) c = 1.0;
if (c < -1.0) c = -1.0; if (c < -1.0) c = -1.0;
@ -162,36 +160,36 @@ void DihedralMultiHarmonic::compute(int eflag, int vflag)
// p = sum (i=1,5) a_i * c**(i-1) // p = sum (i=1,5) a_i * c**(i-1)
// pd = dp/dc // pd = dp/dc
p = a1[type] + c*(a2[type] + c*(a3[type] + c*(a4[type] + c*a5[type]))); p = a1[type] + c * (a2[type] + c * (a3[type] + c * (a4[type] + c * a5[type])));
pd = a2[type] + c*(2.0*a3[type] + c*(3.0*a4[type] + c*4.0*a5[type])); pd = a2[type] + c * (2.0 * a3[type] + c * (3.0 * a4[type] + c * 4.0 * a5[type]));
if (eflag) edihedral = p; if (eflag) edihedral = p;
a = pd; a = pd;
c = c * a; c = c * a;
s12 = s12 * a; s12 = s12 * a;
a11 = c*sb1*s1; a11 = c * sb1 * s1;
a22 = -sb2 * (2.0*c0*s12 - c*(s1+s2)); a22 = -sb2 * (2.0 * c0 * s12 - c * (s1 + s2));
a33 = c*sb3*s2; a33 = c * sb3 * s2;
a12 = -r12c1*(c1mag*c*s1 + c2mag*s12); a12 = -r12c1 * (c1mag * c * s1 + c2mag * s12);
a13 = -rb1*rb3*s12; a13 = -rb1 * rb3 * s12;
a23 = r12c2*(c2mag*c*s2 + c1mag*s12); a23 = r12c2 * (c2mag * c * s2 + c1mag * s12);
sx2 = a12*vb1x + a22*vb2x + a23*vb3x; sx2 = a12 * vb1x + a22 * vb2x + a23 * vb3x;
sy2 = a12*vb1y + a22*vb2y + a23*vb3y; sy2 = a12 * vb1y + a22 * vb2y + a23 * vb3y;
sz2 = a12*vb1z + a22*vb2z + a23*vb3z; sz2 = a12 * vb1z + a22 * vb2z + a23 * vb3z;
f1[0] = a11*vb1x + a12*vb2x + a13*vb3x; f1[0] = a11 * vb1x + a12 * vb2x + a13 * vb3x;
f1[1] = a11*vb1y + a12*vb2y + a13*vb3y; f1[1] = a11 * vb1y + a12 * vb2y + a13 * vb3y;
f1[2] = a11*vb1z + a12*vb2z + a13*vb3z; f1[2] = a11 * vb1z + a12 * vb2z + a13 * vb3z;
f2[0] = -sx2 - f1[0]; f2[0] = -sx2 - f1[0];
f2[1] = -sy2 - f1[1]; f2[1] = -sy2 - f1[1];
f2[2] = -sz2 - f1[2]; f2[2] = -sz2 - f1[2];
f4[0] = a13*vb1x + a23*vb2x + a33*vb3x; f4[0] = a13 * vb1x + a23 * vb2x + a33 * vb3x;
f4[1] = a13*vb1y + a23*vb2y + a33*vb3y; f4[1] = a13 * vb1y + a23 * vb2y + a33 * vb3y;
f4[2] = a13*vb1z + a23*vb2z + a33*vb3z; f4[2] = a13 * vb1z + a23 * vb2z + a33 * vb3z;
f3[0] = sx2 - f4[0]; f3[0] = sx2 - f4[0];
f3[1] = sy2 - f4[1]; f3[1] = sy2 - f4[1];
@ -224,8 +222,8 @@ void DihedralMultiHarmonic::compute(int eflag, int vflag)
} }
if (evflag) if (evflag)
ev_tally(i1,i2,i3,i4,nlocal,newton_bond,edihedral,f1,f3,f4, ev_tally(i1, i2, i3, i4, nlocal, newton_bond, edihedral, f1, f3, f4, vb1x, vb1y, vb1z, vb2x,
vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z); vb2y, vb2z, vb3x, vb3y, vb3z);
} }
} }
@ -234,16 +232,16 @@ void DihedralMultiHarmonic::compute(int eflag, int vflag)
void DihedralMultiHarmonic::allocate() void DihedralMultiHarmonic::allocate()
{ {
allocated = 1; allocated = 1;
int n = atom->ndihedraltypes; const int np1 = atom->ndihedraltypes + 1;
memory->create(a1,n+1,"dihedral:a1"); memory->create(a1, np1, "dihedral:a1");
memory->create(a2,n+1,"dihedral:a2"); memory->create(a2, np1, "dihedral:a2");
memory->create(a3,n+1,"dihedral:a3"); memory->create(a3, np1, "dihedral:a3");
memory->create(a4,n+1,"dihedral:a4"); memory->create(a4, np1, "dihedral:a4");
memory->create(a5,n+1,"dihedral:a5"); memory->create(a5, np1, "dihedral:a5");
memory->create(setflag,n+1,"dihedral:setflag"); memory->create(setflag, np1, "dihedral:setflag");
for (int i = 1; i <= n; i++) setflag[i] = 0; for (int i = 1; i < np1; i++) setflag[i] = 0;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -252,17 +250,17 @@ void DihedralMultiHarmonic::allocate()
void DihedralMultiHarmonic::coeff(int narg, char **arg) void DihedralMultiHarmonic::coeff(int narg, char **arg)
{ {
if (narg != 6) error->all(FLERR,"Incorrect args for dihedral coefficients"); if (narg != 6) error->all(FLERR, "Incorrect args for dihedral coefficients");
if (!allocated) allocate(); if (!allocated) allocate();
int ilo,ihi; int ilo, ihi;
utils::bounds(FLERR,arg[0],1,atom->ndihedraltypes,ilo,ihi,error); utils::bounds(FLERR, arg[0], 1, atom->ndihedraltypes, ilo, ihi, error);
double a1_one = utils::numeric(FLERR,arg[1],false,lmp); double a1_one = utils::numeric(FLERR, arg[1], false, lmp);
double a2_one = utils::numeric(FLERR,arg[2],false,lmp); double a2_one = utils::numeric(FLERR, arg[2], false, lmp);
double a3_one = utils::numeric(FLERR,arg[3],false,lmp); double a3_one = utils::numeric(FLERR, arg[3], false, lmp);
double a4_one = utils::numeric(FLERR,arg[4],false,lmp); double a4_one = utils::numeric(FLERR, arg[4], false, lmp);
double a5_one = utils::numeric(FLERR,arg[5],false,lmp); double a5_one = utils::numeric(FLERR, arg[5], false, lmp);
int count = 0; int count = 0;
for (int i = ilo; i <= ihi; i++) { for (int i = ilo; i <= ihi; i++) {
@ -275,7 +273,7 @@ void DihedralMultiHarmonic::coeff(int narg, char **arg)
count++; count++;
} }
if (count == 0) error->all(FLERR,"Incorrect args for dihedral coefficients"); if (count == 0) error->all(FLERR, "Incorrect args for dihedral coefficients");
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -284,11 +282,11 @@ void DihedralMultiHarmonic::coeff(int narg, char **arg)
void DihedralMultiHarmonic::write_restart(FILE *fp) void DihedralMultiHarmonic::write_restart(FILE *fp)
{ {
fwrite(&a1[1],sizeof(double),atom->ndihedraltypes,fp); fwrite(&a1[1], sizeof(double), atom->ndihedraltypes, fp);
fwrite(&a2[1],sizeof(double),atom->ndihedraltypes,fp); fwrite(&a2[1], sizeof(double), atom->ndihedraltypes, fp);
fwrite(&a3[1],sizeof(double),atom->ndihedraltypes,fp); fwrite(&a3[1], sizeof(double), atom->ndihedraltypes, fp);
fwrite(&a4[1],sizeof(double),atom->ndihedraltypes,fp); fwrite(&a4[1], sizeof(double), atom->ndihedraltypes, fp);
fwrite(&a5[1],sizeof(double),atom->ndihedraltypes,fp); fwrite(&a5[1], sizeof(double), atom->ndihedraltypes, fp);
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -300,17 +298,17 @@ void DihedralMultiHarmonic::read_restart(FILE *fp)
allocate(); allocate();
if (comm->me == 0) { if (comm->me == 0) {
utils::sfread(FLERR,&a1[1],sizeof(double),atom->ndihedraltypes,fp,nullptr,error); utils::sfread(FLERR, &a1[1], sizeof(double), atom->ndihedraltypes, fp, nullptr, error);
utils::sfread(FLERR,&a2[1],sizeof(double),atom->ndihedraltypes,fp,nullptr,error); utils::sfread(FLERR, &a2[1], sizeof(double), atom->ndihedraltypes, fp, nullptr, error);
utils::sfread(FLERR,&a3[1],sizeof(double),atom->ndihedraltypes,fp,nullptr,error); utils::sfread(FLERR, &a3[1], sizeof(double), atom->ndihedraltypes, fp, nullptr, error);
utils::sfread(FLERR,&a4[1],sizeof(double),atom->ndihedraltypes,fp,nullptr,error); utils::sfread(FLERR, &a4[1], sizeof(double), atom->ndihedraltypes, fp, nullptr, error);
utils::sfread(FLERR,&a5[1],sizeof(double),atom->ndihedraltypes,fp,nullptr,error); utils::sfread(FLERR, &a5[1], sizeof(double), atom->ndihedraltypes, fp, nullptr, error);
} }
MPI_Bcast(&a1[1],atom->ndihedraltypes,MPI_DOUBLE,0,world); MPI_Bcast(&a1[1], atom->ndihedraltypes, MPI_DOUBLE, 0, world);
MPI_Bcast(&a2[1],atom->ndihedraltypes,MPI_DOUBLE,0,world); MPI_Bcast(&a2[1], atom->ndihedraltypes, MPI_DOUBLE, 0, world);
MPI_Bcast(&a3[1],atom->ndihedraltypes,MPI_DOUBLE,0,world); MPI_Bcast(&a3[1], atom->ndihedraltypes, MPI_DOUBLE, 0, world);
MPI_Bcast(&a4[1],atom->ndihedraltypes,MPI_DOUBLE,0,world); MPI_Bcast(&a4[1], atom->ndihedraltypes, MPI_DOUBLE, 0, world);
MPI_Bcast(&a5[1],atom->ndihedraltypes,MPI_DOUBLE,0,world); MPI_Bcast(&a5[1], atom->ndihedraltypes, MPI_DOUBLE, 0, world);
for (int i = 1; i <= atom->ndihedraltypes; i++) setflag[i] = 1; for (int i = 1; i <= atom->ndihedraltypes; i++) setflag[i] = 1;
} }
@ -322,5 +320,5 @@ void DihedralMultiHarmonic::read_restart(FILE *fp)
void DihedralMultiHarmonic::write_data(FILE *fp) void DihedralMultiHarmonic::write_data(FILE *fp)
{ {
for (int i = 1; i <= atom->ndihedraltypes; i++) for (int i = 1; i <= atom->ndihedraltypes; i++)
fprintf(fp,"%d %g %g %g %g %g\n",i,a1[i],a2[i],a3[i],a4[i],a5[i]); fprintf(fp, "%d %g %g %g %g %g\n", i, a1[i], a2[i], a3[i], a4[i], a5[i]);
} }

181
src/MOLECULE/dihedral_opls.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories https://www.lammps.org/, Sandia National Laboratories
@ -29,13 +28,13 @@
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
#define TOLERANCE 0.05 static constexpr double TOLERANCE = 0.05;
#define SMALL 0.001 static constexpr double SMALL = 0.001;
#define SMALLER 0.00001 static constexpr double SMALLER = 0.00001;
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
DihedralOPLS::DihedralOPLS(LAMMPS *lmp) : Dihedral(lmp) DihedralOPLS::DihedralOPLS(LAMMPS *_lmp) : Dihedral(_lmp)
{ {
writedata = 1; writedata = 1;
} }
@ -57,17 +56,17 @@ DihedralOPLS::~DihedralOPLS()
void DihedralOPLS::compute(int eflag, int vflag) void DihedralOPLS::compute(int eflag, int vflag)
{ {
int i1,i2,i3,i4,n,type; int i1, i2, i3, i4, n, type;
double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,vb2xm,vb2ym,vb2zm; double vb1x, vb1y, vb1z, vb2x, vb2y, vb2z, vb3x, vb3y, vb3z, vb2xm, vb2ym, vb2zm;
double edihedral,f1[3],f2[3],f3[3],f4[3]; double edihedral, f1[3], f2[3], f3[3], f4[3];
double sb1,sb2,sb3,rb1,rb3,c0,b1mag2,b1mag,b2mag2; double sb1, sb2, sb3, rb1, rb3, c0, b1mag2, b1mag, b2mag2;
double b2mag,b3mag2,b3mag,ctmp,r12c1,c1mag,r12c2; double b2mag, b3mag2, b3mag, ctmp, r12c1, c1mag, r12c2;
double c2mag,sc1,sc2,s1,s12,c,p,pd,a,a11,a22; double c2mag, sc1, sc2, s1, s12, c, p, pd, a, a11, a22;
double a33,a12,a13,a23,sx2,sy2,sz2; double a33, a12, a13, a23, sx2, sy2, sz2;
double s2,cx,cy,cz,cmag,dx,phi,si,siinv,sin2; double s2, cx, cy, cz, cmag, dx, phi, si, siinv, sin2;
edihedral = 0.0; edihedral = 0.0;
ev_init(eflag,vflag); ev_init(eflag, vflag);
double **x = atom->x; double **x = atom->x;
double **f = atom->f; double **f = atom->f;
@ -107,59 +106,58 @@ void DihedralOPLS::compute(int eflag, int vflag)
// c0 calculation // c0 calculation
sb1 = 1.0 / (vb1x*vb1x + vb1y*vb1y + vb1z*vb1z); sb1 = 1.0 / (vb1x * vb1x + vb1y * vb1y + vb1z * vb1z);
sb2 = 1.0 / (vb2x*vb2x + vb2y*vb2y + vb2z*vb2z); sb2 = 1.0 / (vb2x * vb2x + vb2y * vb2y + vb2z * vb2z);
sb3 = 1.0 / (vb3x*vb3x + vb3y*vb3y + vb3z*vb3z); sb3 = 1.0 / (vb3x * vb3x + vb3y * vb3y + vb3z * vb3z);
rb1 = sqrt(sb1); rb1 = sqrt(sb1);
rb3 = sqrt(sb3); rb3 = sqrt(sb3);
c0 = (vb1x*vb3x + vb1y*vb3y + vb1z*vb3z) * rb1*rb3; c0 = (vb1x * vb3x + vb1y * vb3y + vb1z * vb3z) * rb1 * rb3;
// 1st and 2nd angle // 1st and 2nd angle
b1mag2 = vb1x*vb1x + vb1y*vb1y + vb1z*vb1z; b1mag2 = vb1x * vb1x + vb1y * vb1y + vb1z * vb1z;
b1mag = sqrt(b1mag2); b1mag = sqrt(b1mag2);
b2mag2 = vb2x*vb2x + vb2y*vb2y + vb2z*vb2z; b2mag2 = vb2x * vb2x + vb2y * vb2y + vb2z * vb2z;
b2mag = sqrt(b2mag2); b2mag = sqrt(b2mag2);
b3mag2 = vb3x*vb3x + vb3y*vb3y + vb3z*vb3z; b3mag2 = vb3x * vb3x + vb3y * vb3y + vb3z * vb3z;
b3mag = sqrt(b3mag2); b3mag = sqrt(b3mag2);
ctmp = vb1x*vb2x + vb1y*vb2y + vb1z*vb2z; ctmp = vb1x * vb2x + vb1y * vb2y + vb1z * vb2z;
r12c1 = 1.0 / (b1mag*b2mag); r12c1 = 1.0 / (b1mag * b2mag);
c1mag = ctmp * r12c1; c1mag = ctmp * r12c1;
ctmp = vb2xm*vb3x + vb2ym*vb3y + vb2zm*vb3z; ctmp = vb2xm * vb3x + vb2ym * vb3y + vb2zm * vb3z;
r12c2 = 1.0 / (b2mag*b3mag); r12c2 = 1.0 / (b2mag * b3mag);
c2mag = ctmp * r12c2; c2mag = ctmp * r12c2;
// cos and sin of 2 angles and final c // cos and sin of 2 angles and final c
sin2 = MAX(1.0 - c1mag*c1mag,0.0); sin2 = MAX(1.0 - c1mag * c1mag, 0.0);
sc1 = sqrt(sin2); sc1 = sqrt(sin2);
if (sc1 < SMALL) sc1 = SMALL; if (sc1 < SMALL) sc1 = SMALL;
sc1 = 1.0/sc1; sc1 = 1.0 / sc1;
sin2 = MAX(1.0 - c2mag*c2mag,0.0); sin2 = MAX(1.0 - c2mag * c2mag, 0.0);
sc2 = sqrt(sin2); sc2 = sqrt(sin2);
if (sc2 < SMALL) sc2 = SMALL; if (sc2 < SMALL) sc2 = SMALL;
sc2 = 1.0/sc2; sc2 = 1.0 / sc2;
s1 = sc1 * sc1; s1 = sc1 * sc1;
s2 = sc2 * sc2; s2 = sc2 * sc2;
s12 = sc1 * sc2; s12 = sc1 * sc2;
c = (c0 + c1mag*c2mag) * s12; c = (c0 + c1mag * c2mag) * s12;
cx = vb1y*vb2z - vb1z*vb2y; cx = vb1y * vb2z - vb1z * vb2y;
cy = vb1z*vb2x - vb1x*vb2z; cy = vb1z * vb2x - vb1x * vb2z;
cz = vb1x*vb2y - vb1y*vb2x; cz = vb1x * vb2y - vb1y * vb2x;
cmag = sqrt(cx*cx + cy*cy + cz*cz); cmag = sqrt(cx * cx + cy * cy + cz * cz);
dx = (cx*vb3x + cy*vb3y + cz*vb3z)/cmag/b3mag; dx = (cx * vb3x + cy * vb3y + cz * vb3z) / cmag / b3mag;
// error check // error check
if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE)) if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE)) problem(FLERR, i1, i2, i3, i4);
problem(FLERR, i1, i2, i3, i4);
if (c > 1.0) c = 1.0; if (c > 1.0) c = 1.0;
if (c < -1.0) c = -1.0; if (c < -1.0) c = -1.0;
@ -172,40 +170,40 @@ void DihedralOPLS::compute(int eflag, int vflag)
if (dx < 0.0) phi *= -1.0; if (dx < 0.0) phi *= -1.0;
si = sin(phi); si = sin(phi);
if (fabs(si) < SMALLER) si = SMALLER; if (fabs(si) < SMALLER) si = SMALLER;
siinv = 1.0/si; siinv = 1.0 / si;
p = k1[type]*(1.0 + c) + k2[type]*(1.0 - cos(2.0*phi)) + p = k1[type] * (1.0 + c) + k2[type] * (1.0 - cos(2.0 * phi)) +
k3[type]*(1.0 + cos(3.0*phi)) + k4[type]*(1.0 - cos(4.0*phi)) ; k3[type] * (1.0 + cos(3.0 * phi)) + k4[type] * (1.0 - cos(4.0 * phi));
pd = k1[type] - 2.0*k2[type]*sin(2.0*phi)*siinv + pd = k1[type] - 2.0 * k2[type] * sin(2.0 * phi) * siinv +
3.0*k3[type]*sin(3.0*phi)*siinv - 4.0*k4[type]*sin(4.0*phi)*siinv; 3.0 * k3[type] * sin(3.0 * phi) * siinv - 4.0 * k4[type] * sin(4.0 * phi) * siinv;
if (eflag) edihedral = p; if (eflag) edihedral = p;
a = pd; a = pd;
c = c * a; c = c * a;
s12 = s12 * a; s12 = s12 * a;
a11 = c*sb1*s1; a11 = c * sb1 * s1;
a22 = -sb2 * (2.0*c0*s12 - c*(s1+s2)); a22 = -sb2 * (2.0 * c0 * s12 - c * (s1 + s2));
a33 = c*sb3*s2; a33 = c * sb3 * s2;
a12 = -r12c1 * (c1mag*c*s1 + c2mag*s12); a12 = -r12c1 * (c1mag * c * s1 + c2mag * s12);
a13 = -rb1*rb3*s12; a13 = -rb1 * rb3 * s12;
a23 = r12c2 * (c2mag*c*s2 + c1mag*s12); a23 = r12c2 * (c2mag * c * s2 + c1mag * s12);
sx2 = a12*vb1x + a22*vb2x + a23*vb3x; sx2 = a12 * vb1x + a22 * vb2x + a23 * vb3x;
sy2 = a12*vb1y + a22*vb2y + a23*vb3y; sy2 = a12 * vb1y + a22 * vb2y + a23 * vb3y;
sz2 = a12*vb1z + a22*vb2z + a23*vb3z; sz2 = a12 * vb1z + a22 * vb2z + a23 * vb3z;
f1[0] = a11*vb1x + a12*vb2x + a13*vb3x; f1[0] = a11 * vb1x + a12 * vb2x + a13 * vb3x;
f1[1] = a11*vb1y + a12*vb2y + a13*vb3y; f1[1] = a11 * vb1y + a12 * vb2y + a13 * vb3y;
f1[2] = a11*vb1z + a12*vb2z + a13*vb3z; f1[2] = a11 * vb1z + a12 * vb2z + a13 * vb3z;
f2[0] = -sx2 - f1[0]; f2[0] = -sx2 - f1[0];
f2[1] = -sy2 - f1[1]; f2[1] = -sy2 - f1[1];
f2[2] = -sz2 - f1[2]; f2[2] = -sz2 - f1[2];
f4[0] = a13*vb1x + a23*vb2x + a33*vb3x; f4[0] = a13 * vb1x + a23 * vb2x + a33 * vb3x;
f4[1] = a13*vb1y + a23*vb2y + a33*vb3y; f4[1] = a13 * vb1y + a23 * vb2y + a33 * vb3y;
f4[2] = a13*vb1z + a23*vb2z + a33*vb3z; f4[2] = a13 * vb1z + a23 * vb2z + a33 * vb3z;
f3[0] = sx2 - f4[0]; f3[0] = sx2 - f4[0];
f3[1] = sy2 - f4[1]; f3[1] = sy2 - f4[1];
@ -238,8 +236,8 @@ void DihedralOPLS::compute(int eflag, int vflag)
} }
if (evflag) if (evflag)
ev_tally(i1,i2,i3,i4,nlocal,newton_bond,edihedral,f1,f3,f4, ev_tally(i1, i2, i3, i4, nlocal, newton_bond, edihedral, f1, f3, f4, vb1x, vb1y, vb1z, vb2x,
vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z); vb2y, vb2z, vb3x, vb3y, vb3z);
} }
} }
@ -248,15 +246,15 @@ void DihedralOPLS::compute(int eflag, int vflag)
void DihedralOPLS::allocate() void DihedralOPLS::allocate()
{ {
allocated = 1; allocated = 1;
int n = atom->ndihedraltypes; const int np1 = atom->ndihedraltypes + 1;
memory->create(k1,n+1,"dihedral:k1"); memory->create(k1, np1, "dihedral:k1");
memory->create(k2,n+1,"dihedral:k2"); memory->create(k2, np1, "dihedral:k2");
memory->create(k3,n+1,"dihedral:k3"); memory->create(k3, np1, "dihedral:k3");
memory->create(k4,n+1,"dihedral:k4"); memory->create(k4, np1, "dihedral:k4");
memory->create(setflag,n+1,"dihedral:setflag"); memory->create(setflag, np1, "dihedral:setflag");
for (int i = 1; i <= n; i++) setflag[i] = 0; for (int i = 1; i < np1; i++) setflag[i] = 0;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -265,30 +263,30 @@ void DihedralOPLS::allocate()
void DihedralOPLS::coeff(int narg, char **arg) void DihedralOPLS::coeff(int narg, char **arg)
{ {
if (narg != 5) error->all(FLERR,"Incorrect args for dihedral coefficients"); if (narg != 5) error->all(FLERR, "Incorrect args for dihedral coefficients");
if (!allocated) allocate(); if (!allocated) allocate();
int ilo,ihi; int ilo, ihi;
utils::bounds(FLERR,arg[0],1,atom->ndihedraltypes,ilo,ihi,error); utils::bounds(FLERR, arg[0], 1, atom->ndihedraltypes, ilo, ihi, error);
double k1_one = utils::numeric(FLERR,arg[1],false,lmp); double k1_one = utils::numeric(FLERR, arg[1], false, lmp);
double k2_one = utils::numeric(FLERR,arg[2],false,lmp); double k2_one = utils::numeric(FLERR, arg[2], false, lmp);
double k3_one = utils::numeric(FLERR,arg[3],false,lmp); double k3_one = utils::numeric(FLERR, arg[3], false, lmp);
double k4_one = utils::numeric(FLERR,arg[4],false,lmp); double k4_one = utils::numeric(FLERR, arg[4], false, lmp);
// store 1/2 factor with prefactor // store 1/2 factor with prefactor
int count = 0; int count = 0;
for (int i = ilo; i <= ihi; i++) { for (int i = ilo; i <= ihi; i++) {
k1[i] = 0.5*k1_one; k1[i] = 0.5 * k1_one;
k2[i] = 0.5*k2_one; k2[i] = 0.5 * k2_one;
k3[i] = 0.5*k3_one; k3[i] = 0.5 * k3_one;
k4[i] = 0.5*k4_one; k4[i] = 0.5 * k4_one;
setflag[i] = 1; setflag[i] = 1;
count++; count++;
} }
if (count == 0) error->all(FLERR,"Incorrect args for dihedral coefficients"); if (count == 0) error->all(FLERR, "Incorrect args for dihedral coefficients");
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -297,10 +295,10 @@ void DihedralOPLS::coeff(int narg, char **arg)
void DihedralOPLS::write_restart(FILE *fp) void DihedralOPLS::write_restart(FILE *fp)
{ {
fwrite(&k1[1],sizeof(double),atom->ndihedraltypes,fp); fwrite(&k1[1], sizeof(double), atom->ndihedraltypes, fp);
fwrite(&k2[1],sizeof(double),atom->ndihedraltypes,fp); fwrite(&k2[1], sizeof(double), atom->ndihedraltypes, fp);
fwrite(&k3[1],sizeof(double),atom->ndihedraltypes,fp); fwrite(&k3[1], sizeof(double), atom->ndihedraltypes, fp);
fwrite(&k4[1],sizeof(double),atom->ndihedraltypes,fp); fwrite(&k4[1], sizeof(double), atom->ndihedraltypes, fp);
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -312,15 +310,15 @@ void DihedralOPLS::read_restart(FILE *fp)
allocate(); allocate();
if (comm->me == 0) { if (comm->me == 0) {
utils::sfread(FLERR,&k1[1],sizeof(double),atom->ndihedraltypes,fp,nullptr,error); utils::sfread(FLERR, &k1[1], sizeof(double), atom->ndihedraltypes, fp, nullptr, error);
utils::sfread(FLERR,&k2[1],sizeof(double),atom->ndihedraltypes,fp,nullptr,error); utils::sfread(FLERR, &k2[1], sizeof(double), atom->ndihedraltypes, fp, nullptr, error);
utils::sfread(FLERR,&k3[1],sizeof(double),atom->ndihedraltypes,fp,nullptr,error); utils::sfread(FLERR, &k3[1], sizeof(double), atom->ndihedraltypes, fp, nullptr, error);
utils::sfread(FLERR,&k4[1],sizeof(double),atom->ndihedraltypes,fp,nullptr,error); utils::sfread(FLERR, &k4[1], sizeof(double), atom->ndihedraltypes, fp, nullptr, error);
} }
MPI_Bcast(&k1[1],atom->ndihedraltypes,MPI_DOUBLE,0,world); MPI_Bcast(&k1[1], atom->ndihedraltypes, MPI_DOUBLE, 0, world);
MPI_Bcast(&k2[1],atom->ndihedraltypes,MPI_DOUBLE,0,world); MPI_Bcast(&k2[1], atom->ndihedraltypes, MPI_DOUBLE, 0, world);
MPI_Bcast(&k3[1],atom->ndihedraltypes,MPI_DOUBLE,0,world); MPI_Bcast(&k3[1], atom->ndihedraltypes, MPI_DOUBLE, 0, world);
MPI_Bcast(&k4[1],atom->ndihedraltypes,MPI_DOUBLE,0,world); MPI_Bcast(&k4[1], atom->ndihedraltypes, MPI_DOUBLE, 0, world);
for (int i = 1; i <= atom->ndihedraltypes; i++) setflag[i] = 1; for (int i = 1; i <= atom->ndihedraltypes; i++) setflag[i] = 1;
} }
@ -332,6 +330,5 @@ void DihedralOPLS::read_restart(FILE *fp)
void DihedralOPLS::write_data(FILE *fp) void DihedralOPLS::write_data(FILE *fp)
{ {
for (int i = 1; i <= atom->ndihedraltypes; i++) for (int i = 1; i <= atom->ndihedraltypes; i++)
fprintf(fp,"%d %g %g %g %g\n",i,2.0*k1[i],2.0*k2[i],2.0*k3[i],2.0*k4[i]); fprintf(fp, "%d %g %g %g %g\n", i, 2.0 * k1[i], 2.0 * k2[i], 2.0 * k3[i], 2.0 * k4[i]);
} }