From 85e581eb2ac70548f1da3d234b3f76edb25c8d17 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 17 Mar 2022 19:44:53 -0400 Subject: [PATCH] enable and apply clang-format --- src/MOLECULE/dihedral_charmm.cpp | 264 ++++++++++---------- src/MOLECULE/dihedral_charmmfsw.cpp | 304 +++++++++++------------ src/MOLECULE/dihedral_harmonic.cpp | 170 +++++++------ src/MOLECULE/dihedral_multi_harmonic.cpp | 164 ++++++------ src/MOLECULE/dihedral_opls.cpp | 181 +++++++------- 5 files changed, 535 insertions(+), 548 deletions(-) diff --git a/src/MOLECULE/dihedral_charmm.cpp b/src/MOLECULE/dihedral_charmm.cpp index 62a3c30434..ee053f4f12 100644 --- a/src/MOLECULE/dihedral_charmm.cpp +++ b/src/MOLECULE/dihedral_charmm.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -32,13 +31,13 @@ #include 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; writedata = 1; @@ -63,19 +62,19 @@ DihedralCharmm::~DihedralCharmm() void DihedralCharmm::compute(int eflag, int vflag) { - int i1,i2,i3,i4,i,m,n,type; - double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,vb2xm,vb2ym,vb2zm; - 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 df,df1,ddf1,fg,hg,fga,hgb,gaa,gbb; - double dtfx,dtfy,dtfz,dtgx,dtgy,dtgz,dthx,dthy,dthz; - double c,s,p,sx2,sy2,sz2; - int itype,jtype; - double delx,dely,delz,rsq,r2inv,r6inv; - double forcecoul,forcelj,fpair,ecoul,evdwl; + int i1, i2, i3, i4, i, m, n, type; + double vb1x, vb1y, vb1z, vb2x, vb2y, vb2z, vb3x, vb3y, vb3z, vb2xm, vb2ym, vb2zm; + 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 df, df1, ddf1, fg, hg, fga, hgb, gaa, gbb; + double dtfx, dtfy, dtfz, dtgx, dtgy, dtgz, dthx, dthy, dthz; + double c, s, p, sx2, sy2, sz2; + int itype, jtype; + double delx, dely, delz, rsq, r2inv, r6inv; + double forcecoul, forcelj, fpair, ecoul, evdwl; edihedral = evdwl = ecoul = 0.0; - ev_init(eflag,vflag); + ev_init(eflag, vflag); // 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]; vb3z = x[i4][2] - x[i3][2]; - ax = vb1y*vb2zm - vb1z*vb2ym; - ay = vb1z*vb2xm - vb1x*vb2zm; - az = vb1x*vb2ym - vb1y*vb2xm; - bx = vb3y*vb2zm - vb3z*vb2ym; - by = vb3z*vb2xm - vb3x*vb2zm; - bz = vb3x*vb2ym - vb3y*vb2xm; + ax = vb1y * vb2zm - vb1z * vb2ym; + ay = vb1z * vb2xm - vb1x * vb2zm; + az = vb1x * vb2ym - vb1y * vb2xm; + bx = vb3y * vb2zm - vb3z * vb2ym; + by = vb3z * vb2xm - vb3x * vb2zm; + bz = vb3x * vb2ym - vb3y * vb2xm; - rasq = ax*ax + ay*ay + az*az; - rbsq = bx*bx + by*by + bz*bz; - rgsq = vb2xm*vb2xm + vb2ym*vb2ym + vb2zm*vb2zm; + rasq = ax * ax + ay * ay + az * az; + rbsq = bx * bx + by * by + bz * bz; + rgsq = vb2xm * vb2xm + vb2ym * vb2ym + vb2zm * vb2zm; rg = sqrt(rgsq); rginv = ra2inv = rb2inv = 0.0; - if (rg > 0) rginv = 1.0/rg; - if (rasq > 0) ra2inv = 1.0/rasq; - if (rbsq > 0) rb2inv = 1.0/rbsq; - rabinv = sqrt(ra2inv*rb2inv); + if (rg > 0) rginv = 1.0 / rg; + if (rasq > 0) ra2inv = 1.0 / rasq; + if (rbsq > 0) rb2inv = 1.0 / rbsq; + rabinv = sqrt(ra2inv * rb2inv); - c = (ax*bx + ay*by + az*bz)*rabinv; - s = rg*rabinv*(ax*vb3x + ay*vb3y + az*vb3z); + c = (ax * bx + ay * by + az * bz) * rabinv; + s = rg * rabinv * (ax * vb3x + ay * vb3y + az * vb3z); // error check - if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE)) - problem(FLERR, i1, i2, i3, i4); + if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE)) problem(FLERR, i1, i2, i3, i4); 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; for (i = 0; i < m; i++) { - ddf1 = p*c - df1*s; - df1 = p*s + df1*c; + ddf1 = p * c - df1 * s; + df1 = p * s + df1 * c; p = ddf1; } - p = p*cos_shift[type] + df1*sin_shift[type]; - df1 = df1*cos_shift[type] - ddf1*sin_shift[type]; + p = p * cos_shift[type] + df1 * sin_shift[type]; + df1 = df1 * cos_shift[type] - ddf1 * sin_shift[type]; df1 *= -m; p += 1.0; @@ -172,40 +170,40 @@ void DihedralCharmm::compute(int eflag, int vflag) if (eflag) edihedral = k[type] * p; - fg = vb1x*vb2xm + vb1y*vb2ym + vb1z*vb2zm; - hg = vb3x*vb2xm + vb3y*vb2ym + vb3z*vb2zm; - fga = fg*ra2inv*rginv; - hgb = hg*rb2inv*rginv; - gaa = -ra2inv*rg; - gbb = rb2inv*rg; + fg = vb1x * vb2xm + vb1y * vb2ym + vb1z * vb2zm; + hg = vb3x * vb2xm + vb3y * vb2ym + vb3z * vb2zm; + fga = fg * ra2inv * rginv; + hgb = hg * rb2inv * rginv; + gaa = -ra2inv * rg; + gbb = rb2inv * rg; - dtfx = gaa*ax; - dtfy = gaa*ay; - dtfz = gaa*az; - dtgx = fga*ax - hgb*bx; - dtgy = fga*ay - hgb*by; - dtgz = fga*az - hgb*bz; - dthx = gbb*bx; - dthy = gbb*by; - dthz = gbb*bz; + dtfx = gaa * ax; + dtfy = gaa * ay; + dtfz = gaa * az; + dtgx = fga * ax - hgb * bx; + dtgy = fga * ay - hgb * by; + dtgz = fga * az - hgb * bz; + dthx = gbb * bx; + dthy = gbb * by; + dthz = gbb * bz; df = -k[type] * df1; - sx2 = df*dtgx; - sy2 = df*dtgy; - sz2 = df*dtgz; + sx2 = df * dtgx; + sy2 = df * dtgy; + sz2 = df * dtgz; - f1[0] = df*dtfx; - f1[1] = df*dtfy; - f1[2] = df*dtfz; + f1[0] = df * dtfx; + f1[1] = df * dtfy; + f1[2] = df * dtfz; f2[0] = sx2 - f1[0]; f2[1] = sy2 - f1[1]; f2[2] = sz2 - f1[2]; - f4[0] = df*dthx; - f4[1] = df*dthy; - f4[2] = df*dthz; + f4[0] = df * dthx; + f4[1] = df * dthy; + f4[2] = df * dthz; f3[0] = -sx2 - f4[0]; f3[1] = -sy2 - f4[1]; @@ -238,8 +236,8 @@ void DihedralCharmm::compute(int eflag, int vflag) } if (evflag) - ev_tally(i1,i2,i3,i4,nlocal,newton_bond,edihedral,f1,f3,f4, - vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z); + ev_tally(i1, i2, i3, i4, nlocal, newton_bond, edihedral, f1, f3, f4, vb1x, vb1y, vb1z, vb2x, + vb2y, vb2z, vb3x, vb3y, vb3z); // 1-4 LJ and Coulomb interactions // 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]; dely = x[i1][1] - x[i4][1]; delz = x[i1][2] - x[i4][2]; - rsq = delx*delx + dely*dely + delz*delz; - r2inv = 1.0/rsq; - r6inv = r2inv*r2inv*r2inv; + rsq = delx * delx + dely * dely + delz * delz; + r2inv = 1.0 / rsq; + r6inv = r2inv * r2inv * r2inv; - if (implicit) forcecoul = qqrd2e * q[i1]*q[i4]*r2inv; - else 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 (implicit) + forcecoul = qqrd2e * q[i1] * q[i4] * r2inv; + else + 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) { 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]; } if (newton_bond || i1 < nlocal) { - f[i1][0] += delx*fpair; - f[i1][1] += dely*fpair; - f[i1][2] += delz*fpair; + f[i1][0] += delx * fpair; + f[i1][1] += dely * fpair; + f[i1][2] += delz * fpair; } if (newton_bond || i4 < nlocal) { - f[i4][0] -= delx*fpair; - f[i4][1] -= dely*fpair; - f[i4][2] -= delz*fpair; + f[i4][0] -= delx * fpair; + f[i4][1] -= dely * fpair; + f[i4][2] -= delz * fpair; } - if (evflag) force->pair->ev_tally(i1,i4,nlocal,newton_bond, - evdwl,ecoul,fpair,delx,dely,delz); + if (evflag) + 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() { allocated = 1; - int n = atom->ndihedraltypes; + const int np1 = atom->ndihedraltypes + 1; - memory->create(k,n+1,"dihedral:k"); - memory->create(multiplicity,n+1,"dihedral:multiplicity"); - memory->create(shift,n+1,"dihedral:shift"); - memory->create(cos_shift,n+1,"dihedral:cos_shift"); - memory->create(sin_shift,n+1,"dihedral:sin_shift"); - memory->create(weight,n+1,"dihedral:weight"); + memory->create(k, np1, "dihedral:k"); + memory->create(multiplicity, np1, "dihedral:multiplicity"); + memory->create(shift, np1, "dihedral:shift"); + memory->create(cos_shift, np1, "dihedral:cos_shift"); + memory->create(sin_shift, np1, "dihedral:sin_shift"); + memory->create(weight, np1, "dihedral:weight"); - memory->create(setflag,n+1,"dihedral:setflag"); - for (int i = 1; i <= n; i++) setflag[i] = 0; + memory->create(setflag, np1, "dihedral:setflag"); + for (int i = 1; i < np1; i++) setflag[i] = 0; } /* ---------------------------------------------------------------------- @@ -307,40 +307,40 @@ void DihedralCharmm::allocate() 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(); - int ilo,ihi; - utils::bounds(FLERR,arg[0],1,atom->ndihedraltypes,ilo,ihi,error); + int ilo, ihi; + utils::bounds(FLERR, arg[0], 1, atom->ndihedraltypes, ilo, ihi, error); // require integer values of shift for backwards compatibility // arbitrary phase angle shift could be allowed, but would break // backwards compatibility and is probably not needed - double k_one = utils::numeric(FLERR,arg[1],false,lmp); - int multiplicity_one = utils::inumeric(FLERR,arg[2],false,lmp); - int shift_one = utils::inumeric(FLERR,arg[3],false,lmp); - double weight_one = utils::numeric(FLERR,arg[4],false,lmp); + double k_one = utils::numeric(FLERR, arg[1], false, lmp); + int multiplicity_one = utils::inumeric(FLERR, arg[2], false, lmp); + int shift_one = utils::inumeric(FLERR, arg[3], false, lmp); + double weight_one = utils::numeric(FLERR, arg[4], false, lmp); 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) - error->all(FLERR,"Incorrect weight arg for dihedral coefficients"); - if (weight_one > 0.0) weightflag=1; + error->all(FLERR, "Incorrect weight arg for dihedral coefficients"); + if (weight_one > 0.0) weightflag = 1; int count = 0; for (int i = ilo; i <= ihi; i++) { k[i] = k_one; shift[i] = shift_one; - cos_shift[i] = cos(MY_PI*shift_one/180.0); - sin_shift[i] = sin(MY_PI*shift_one/180.0); + cos_shift[i] = cos(DEG2RAD * shift_one); + sin_shift[i] = sin(DEG2RAD * shift_one); multiplicity[i] = multiplicity_one; weight[i] = weight_one; setflag[i] = 1; 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() { - if (utils::strmatch(update->integrate_style,"^respa")) { + if (utils::strmatch(update->integrate_style, "^respa")) { Respa *r = (Respa *) update->integrate; if (r->level_pair >= 0 && (r->level_pair != r->level_dihedral)) - error->all(FLERR,"Dihedral style charmm must be set to same" - " r-RESPA level as 'pair'"); + error->all(FLERR, "Dihedral style charmm must be set to same r-RESPA level as 'pair'"); if (r->level_outer >= 0 && (r->level_outer != r->level_dihedral)) - error->all(FLERR,"Dihedral style charmm must be set to same" - " r-RESPA level as 'outer'"); + error->all(FLERR, "Dihedral style charmm must be set to same r-RESPA level as 'outer'"); } // insure use of CHARMM pair_style if any weight factors are non-zero @@ -366,19 +364,20 @@ void DihedralCharmm::init_style() if (weightflag) { if ((force->special_lj[3] != 0.0) || (force->special_coul[3] != 0.0)) - error->all(FLERR,"Must use 'special_bonds charmm' with" - " dihedral style charmm for use with CHARMM pair styles"); + error->all(FLERR, + "Must use 'special_bonds charmm' with dihedral " + "style charmm for use with CHARMM pair styles"); int itmp; if (force->pair == nullptr) - error->all(FLERR,"Dihedral charmm is incompatible with Pair style"); - lj14_1 = (double **) force->pair->extract("lj14_1",itmp); - lj14_2 = (double **) force->pair->extract("lj14_2",itmp); - lj14_3 = (double **) force->pair->extract("lj14_3",itmp); - lj14_4 = (double **) force->pair->extract("lj14_4",itmp); - int *ptr = (int *) force->pair->extract("implicit",itmp); + error->all(FLERR, "Dihedral charmm is incompatible with Pair style"); + lj14_1 = (double **) force->pair->extract("lj14_1", itmp); + lj14_2 = (double **) force->pair->extract("lj14_2", itmp); + lj14_3 = (double **) force->pair->extract("lj14_3", itmp); + lj14_4 = (double **) force->pair->extract("lj14_4", itmp); + int *ptr = (int *) force->pair->extract("implicit", itmp); 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; } } @@ -389,11 +388,11 @@ void DihedralCharmm::init_style() void DihedralCharmm::write_restart(FILE *fp) { - fwrite(&k[1],sizeof(double),atom->ndihedraltypes,fp); - fwrite(&multiplicity[1],sizeof(int),atom->ndihedraltypes,fp); - fwrite(&shift[1],sizeof(int),atom->ndihedraltypes,fp); - fwrite(&weight[1],sizeof(double),atom->ndihedraltypes,fp); - fwrite(&weightflag,sizeof(int),1,fp); + fwrite(&k[1], sizeof(double), atom->ndihedraltypes, fp); + fwrite(&multiplicity[1], sizeof(int), atom->ndihedraltypes, fp); + fwrite(&shift[1], sizeof(int), atom->ndihedraltypes, fp); + fwrite(&weight[1], sizeof(double), atom->ndihedraltypes, fp); + fwrite(&weightflag, sizeof(int), 1, fp); } /* ---------------------------------------------------------------------- @@ -405,22 +404,22 @@ void DihedralCharmm::read_restart(FILE *fp) allocate(); if (comm->me == 0) { - 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,&shift[1],sizeof(int),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, &k[1], sizeof(double), 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, &weight[1], sizeof(double), atom->ndihedraltypes, 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(&multiplicity[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(&weightflag,1,MPI_INT,0,world); + MPI_Bcast(&k[1], atom->ndihedraltypes, MPI_DOUBLE, 0, world); + MPI_Bcast(&multiplicity[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(&weightflag, 1, MPI_INT, 0, world); for (int i = 1; i <= atom->ndihedraltypes; i++) { setflag[i] = 1; - cos_shift[i] = cos(MY_PI*shift[i]/180.0); - sin_shift[i] = sin(MY_PI*shift[i]/180.0); + cos_shift[i] = cos(DEG2RAD * shift[i]); + sin_shift[i] = sin(DEG2RAD * shift[i]); } } @@ -431,6 +430,5 @@ void DihedralCharmm::read_restart(FILE *fp) void DihedralCharmm::write_data(FILE *fp) { 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]); } - diff --git a/src/MOLECULE/dihedral_charmmfsw.cpp b/src/MOLECULE/dihedral_charmmfsw.cpp index b74357ab39..4b7f0b04f3 100644 --- a/src/MOLECULE/dihedral_charmmfsw.cpp +++ b/src/MOLECULE/dihedral_charmmfsw.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -35,13 +34,13 @@ #include 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; writedata = 1; @@ -66,19 +65,19 @@ DihedralCharmmfsw::~DihedralCharmmfsw() void DihedralCharmmfsw::compute(int eflag, int vflag) { - int i1,i2,i3,i4,i,m,n,type; - double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,vb2xm,vb2ym,vb2zm; - 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 df,df1,ddf1,fg,hg,fga,hgb,gaa,gbb; - double dtfx,dtfy,dtfz,dtgx,dtgy,dtgz,dthx,dthy,dthz; - double c,s,p,sx2,sy2,sz2; - int itype,jtype; - double delx,dely,delz,rsq,r2inv,r6inv,r; - double forcecoul,forcelj,fpair,ecoul,evdwl; + int i1, i2, i3, i4, i, m, n, type; + double vb1x, vb1y, vb1z, vb2x, vb2y, vb2z, vb3x, vb3y, vb3z, vb2xm, vb2ym, vb2zm; + 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 df, df1, ddf1, fg, hg, fga, hgb, gaa, gbb; + double dtfx, dtfy, dtfz, dtgx, dtgy, dtgz, dthx, dthy, dthz; + double c, s, p, sx2, sy2, sz2; + int itype, jtype; + double delx, dely, delz, rsq, r2inv, r6inv, r; + double forcecoul, forcelj, fpair, ecoul, evdwl; edihedral = evdwl = ecoul = 0.0; - ev_init(eflag,vflag); + ev_init(eflag, vflag); // 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]; vb3z = x[i4][2] - x[i3][2]; - ax = vb1y*vb2zm - vb1z*vb2ym; - ay = vb1z*vb2xm - vb1x*vb2zm; - az = vb1x*vb2ym - vb1y*vb2xm; - bx = vb3y*vb2zm - vb3z*vb2ym; - by = vb3z*vb2xm - vb3x*vb2zm; - bz = vb3x*vb2ym - vb3y*vb2xm; + ax = vb1y * vb2zm - vb1z * vb2ym; + ay = vb1z * vb2xm - vb1x * vb2zm; + az = vb1x * vb2ym - vb1y * vb2xm; + bx = vb3y * vb2zm - vb3z * vb2ym; + by = vb3z * vb2xm - vb3x * vb2zm; + bz = vb3x * vb2ym - vb3y * vb2xm; - rasq = ax*ax + ay*ay + az*az; - rbsq = bx*bx + by*by + bz*bz; - rgsq = vb2xm*vb2xm + vb2ym*vb2ym + vb2zm*vb2zm; + rasq = ax * ax + ay * ay + az * az; + rbsq = bx * bx + by * by + bz * bz; + rgsq = vb2xm * vb2xm + vb2ym * vb2ym + vb2zm * vb2zm; rg = sqrt(rgsq); rginv = ra2inv = rb2inv = 0.0; - if (rg > 0) rginv = 1.0/rg; - if (rasq > 0) ra2inv = 1.0/rasq; - if (rbsq > 0) rb2inv = 1.0/rbsq; - rabinv = sqrt(ra2inv*rb2inv); + if (rg > 0) rginv = 1.0 / rg; + if (rasq > 0) ra2inv = 1.0 / rasq; + if (rbsq > 0) rb2inv = 1.0 / rbsq; + rabinv = sqrt(ra2inv * rb2inv); - c = (ax*bx + ay*by + az*bz)*rabinv; - s = rg*rabinv*(ax*vb3x + ay*vb3y + az*vb3z); + c = (ax * bx + ay * by + az * bz) * rabinv; + s = rg * rabinv * (ax * vb3x + ay * vb3y + az * vb3z); // error check - if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE)) - problem(FLERR, i1, i2, i3, i4); + if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE)) problem(FLERR, i1, i2, i3, i4); 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; for (i = 0; i < m; i++) { - ddf1 = p*c - df1*s; - df1 = p*s + df1*c; + ddf1 = p * c - df1 * s; + df1 = p * s + df1 * c; p = ddf1; } - p = p*cos_shift[type] + df1*sin_shift[type]; - df1 = df1*cos_shift[type] - ddf1*sin_shift[type]; + p = p * cos_shift[type] + df1 * sin_shift[type]; + df1 = df1 * cos_shift[type] - ddf1 * sin_shift[type]; df1 *= -m; p += 1.0; @@ -175,40 +173,40 @@ void DihedralCharmmfsw::compute(int eflag, int vflag) if (eflag) edihedral = k[type] * p; - fg = vb1x*vb2xm + vb1y*vb2ym + vb1z*vb2zm; - hg = vb3x*vb2xm + vb3y*vb2ym + vb3z*vb2zm; - fga = fg*ra2inv*rginv; - hgb = hg*rb2inv*rginv; - gaa = -ra2inv*rg; - gbb = rb2inv*rg; + fg = vb1x * vb2xm + vb1y * vb2ym + vb1z * vb2zm; + hg = vb3x * vb2xm + vb3y * vb2ym + vb3z * vb2zm; + fga = fg * ra2inv * rginv; + hgb = hg * rb2inv * rginv; + gaa = -ra2inv * rg; + gbb = rb2inv * rg; - dtfx = gaa*ax; - dtfy = gaa*ay; - dtfz = gaa*az; - dtgx = fga*ax - hgb*bx; - dtgy = fga*ay - hgb*by; - dtgz = fga*az - hgb*bz; - dthx = gbb*bx; - dthy = gbb*by; - dthz = gbb*bz; + dtfx = gaa * ax; + dtfy = gaa * ay; + dtfz = gaa * az; + dtgx = fga * ax - hgb * bx; + dtgy = fga * ay - hgb * by; + dtgz = fga * az - hgb * bz; + dthx = gbb * bx; + dthy = gbb * by; + dthz = gbb * bz; df = -k[type] * df1; - sx2 = df*dtgx; - sy2 = df*dtgy; - sz2 = df*dtgz; + sx2 = df * dtgx; + sy2 = df * dtgy; + sz2 = df * dtgz; - f1[0] = df*dtfx; - f1[1] = df*dtfy; - f1[2] = df*dtfz; + f1[0] = df * dtfx; + f1[1] = df * dtfy; + f1[2] = df * dtfz; f2[0] = sx2 - f1[0]; f2[1] = sy2 - f1[1]; f2[2] = sz2 - f1[2]; - f4[0] = df*dthx; - f4[1] = df*dthy; - f4[2] = df*dthz; + f4[0] = df * dthx; + f4[1] = df * dthy; + f4[2] = df * dthz; f3[0] = -sx2 - f4[0]; f3[1] = -sy2 - f4[1]; @@ -241,8 +239,8 @@ void DihedralCharmmfsw::compute(int eflag, int vflag) } if (evflag) - ev_tally(i1,i2,i3,i4,nlocal,newton_bond,edihedral,f1,f3,f4, - vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z); + ev_tally(i1, i2, i3, i4, nlocal, newton_bond, edihedral, f1, f3, f4, vb1x, vb1y, vb1z, vb2x, + vb2y, vb2z, vb3x, vb3y, vb3z); // 1-4 LJ and Coulomb interactions // 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]; dely = x[i1][1] - x[i4][1]; delz = x[i1][2] - x[i4][2]; - rsq = delx*delx + dely*dely + delz*delz; - r2inv = 1.0/rsq; - r6inv = r2inv*r2inv*r2inv; + rsq = delx * delx + dely * dely + delz * delz; + r2inv = 1.0 / rsq; + r6inv = r2inv * r2inv * r2inv; // modifying coul and LJ force and energies to apply // 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 r = sqrt(rsq); - if (implicit) forcecoul = qqrd2e * q[i1]*q[i4]*r2inv; - else if (dihedflag) forcecoul = qqrd2e * q[i1]*q[i4]*sqrt(r2inv); - else 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 (implicit) + forcecoul = qqrd2e * q[i1] * q[i4] * r2inv; + else if (dihedflag) + forcecoul = qqrd2e * q[i1] * q[i4] * sqrt(r2inv); + else + 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 (dihedflag) ecoul = weight[type] * forcecoul; - else ecoul = weight[type] * qqrd2e * q[i1]*q[i4] * - (sqrt(r2inv) + r*cut_coulinv14*cut_coulinv14 - - 2.0*cut_coulinv14); - evdwl14_12 = r6inv*lj14_3[itype][jtype]*r6inv - - lj14_3[itype][jtype]*cut_lj_inner6inv*cut_lj6inv; - evdwl14_6 = -lj14_4[itype][jtype]*r6inv + - lj14_4[itype][jtype]*cut_lj_inner3inv*cut_lj3inv; + if (dihedflag) + ecoul = weight[type] * forcecoul; + else + ecoul = weight[type] * qqrd2e * q[i1] * q[i4] * + (sqrt(r2inv) + r * cut_coulinv14 * cut_coulinv14 - 2.0 * cut_coulinv14); + evdwl14_12 = r6inv * lj14_3[itype][jtype] * r6inv - + lj14_3[itype][jtype] * cut_lj_inner6inv * cut_lj6inv; + evdwl14_6 = + -lj14_4[itype][jtype] * r6inv + lj14_4[itype][jtype] * cut_lj_inner3inv * cut_lj3inv; evdwl = evdwl14_12 + evdwl14_6; evdwl *= weight[type]; } if (newton_bond || i1 < nlocal) { - f[i1][0] += delx*fpair; - f[i1][1] += dely*fpair; - f[i1][2] += delz*fpair; + f[i1][0] += delx * fpair; + f[i1][1] += dely * fpair; + f[i1][2] += delz * fpair; } if (newton_bond || i4 < nlocal) { - f[i4][0] -= delx*fpair; - f[i4][1] -= dely*fpair; - f[i4][2] -= delz*fpair; + f[i4][0] -= delx * fpair; + f[i4][1] -= dely * fpair; + f[i4][2] -= delz * fpair; } - if (evflag) force->pair->ev_tally(i1,i4,nlocal,newton_bond, - evdwl,ecoul,fpair,delx,dely,delz); + if (evflag) + 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() { allocated = 1; - int n = atom->ndihedraltypes; + const int np1 = atom->ndihedraltypes + 1; - memory->create(k,n+1,"dihedral:k"); - memory->create(multiplicity,n+1,"dihedral:multiplicity"); - memory->create(shift,n+1,"dihedral:shift"); - memory->create(cos_shift,n+1,"dihedral:cos_shift"); - memory->create(sin_shift,n+1,"dihedral:sin_shift"); - memory->create(weight,n+1,"dihedral:weight"); + memory->create(k, np1, "dihedral:k"); + memory->create(multiplicity, np1, "dihedral:multiplicity"); + memory->create(shift, np1, "dihedral:shift"); + memory->create(cos_shift, np1, "dihedral:cos_shift"); + memory->create(sin_shift, np1, "dihedral:sin_shift"); + memory->create(weight, np1, "dihedral:weight"); - memory->create(setflag,n+1,"dihedral:setflag"); - for (int i = 1; i <= n; i++) setflag[i] = 0; + memory->create(setflag, np1, "dihedral:setflag"); + for (int i = 1; i < np1; i++) setflag[i] = 0; } /* ---------------------------------------------------------------------- @@ -325,40 +326,40 @@ void DihedralCharmmfsw::allocate() 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(); - int ilo,ihi; - utils::bounds(FLERR,arg[0],1,atom->ndihedraltypes,ilo,ihi,error); + int ilo, ihi; + utils::bounds(FLERR, arg[0], 1, atom->ndihedraltypes, ilo, ihi, error); // require integer values of shift for backwards compatibility // arbitrary phase angle shift could be allowed, but would break // backwards compatibility and is probably not needed - double k_one = utils::numeric(FLERR,arg[1],false,lmp); - int multiplicity_one = utils::inumeric(FLERR,arg[2],false,lmp); - int shift_one = utils::inumeric(FLERR,arg[3],false,lmp); - double weight_one = utils::numeric(FLERR,arg[4],false,lmp); + double k_one = utils::numeric(FLERR, arg[1], false, lmp); + int multiplicity_one = utils::inumeric(FLERR, arg[2], false, lmp); + int shift_one = utils::inumeric(FLERR, arg[3], false, lmp); + double weight_one = utils::numeric(FLERR, arg[4], false, lmp); 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) - error->all(FLERR,"Incorrect weight arg for dihedral coefficients"); - if (weight_one > 0.0) weightflag=1; + error->all(FLERR, "Incorrect weight arg for dihedral coefficients"); + if (weight_one > 0.0) weightflag = 1; int count = 0; for (int i = ilo; i <= ihi; i++) { k[i] = k_one; shift[i] = shift_one; - cos_shift[i] = cos(MY_PI*shift_one/180.0); - sin_shift[i] = sin(MY_PI*shift_one/180.0); + cos_shift[i] = cos(DEG2RAD * shift_one); + sin_shift[i] = sin(DEG2RAD * shift_one); multiplicity[i] = multiplicity_one; weight[i] = weight_one; setflag[i] = 1; 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() { - if (utils::strmatch(update->integrate_style,"^respa")) { + if (utils::strmatch(update->integrate_style, "^respa")) { Respa *r = (Respa *) update->integrate; if (r->level_pair >= 0 && (r->level_pair != r->level_dihedral)) - error->all(FLERR,"Dihedral style charmmfsw must be set to same" - " r-RESPA level as 'pair'"); + error->all(FLERR, "Dihedral style charmmfsw must be set to same r-RESPA level as 'pair'"); if (r->level_outer >= 0 && (r->level_outer != r->level_dihedral)) - error->all(FLERR,"Dihedral style charmmfsw must be set to same" - " r-RESPA level as 'outer'"); + error->all(FLERR, "Dihedral style charmmfsw must be set to same r-RESPA level as 'outer'"); } // insure use of CHARMM pair_style if any weight factors are non-zero @@ -384,19 +383,20 @@ void DihedralCharmmfsw::init_style() if (weightflag) { if ((force->special_lj[3] != 0.0) || (force->special_coul[3] != 0.0)) - error->all(FLERR,"Must use 'special_bonds charmm' with" - " dihedral style charmm for use with CHARMM pair styles"); + error->all(FLERR, + "Must use 'special_bonds charmm' with dihedral " + "style charmm for use with CHARMM pair styles"); int itmp; if (force->pair == nullptr) - error->all(FLERR,"Dihedral charmmfsw is incompatible with Pair style"); - lj14_1 = (double **) force->pair->extract("lj14_1",itmp); - lj14_2 = (double **) force->pair->extract("lj14_2",itmp); - lj14_3 = (double **) force->pair->extract("lj14_3",itmp); - lj14_4 = (double **) force->pair->extract("lj14_4",itmp); - int *ptr = (int *) force->pair->extract("implicit",itmp); + error->all(FLERR, "Dihedral charmmfsw is incompatible with Pair style"); + lj14_1 = (double **) force->pair->extract("lj14_1", itmp); + lj14_2 = (double **) force->pair->extract("lj14_2", itmp); + lj14_3 = (double **) force->pair->extract("lj14_3", itmp); + lj14_4 = (double **) force->pair->extract("lj14_4", itmp); + int *ptr = (int *) force->pair->extract("implicit", itmp); 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; } @@ -404,25 +404,24 @@ void DihedralCharmmfsw::init_style() // to 1/4 dihedral atoms to match CHARMM pairwise interactions int itmp; - int *p_dihedflag = (int *) force->pair->extract("dihedflag",itmp); - double *p_cutljinner = (double *) force->pair->extract("cut_lj_inner",itmp); - double *p_cutlj = (double *) force->pair->extract("cut_lj",itmp); - double *p_cutcoul = (double *) force->pair->extract("cut_coul",itmp); + int *p_dihedflag = (int *) force->pair->extract("dihedflag", itmp); + double *p_cutljinner = (double *) force->pair->extract("cut_lj_inner", itmp); + double *p_cutlj = (double *) force->pair->extract("cut_lj", itmp); + double *p_cutcoul = (double *) force->pair->extract("cut_coul", itmp); - if (p_cutcoul == nullptr || p_cutljinner == nullptr || - p_cutlj == nullptr || p_dihedflag == nullptr) - error->all(FLERR,"Dihedral charmmfsw is incompatible with Pair style"); + if (p_cutcoul == nullptr || p_cutljinner == nullptr || p_cutlj == nullptr || + p_dihedflag == nullptr) + error->all(FLERR, "Dihedral charmmfsw is incompatible with Pair style"); dihedflag = *p_dihedflag; cut_coul14 = *p_cutcoul; cut_lj_inner14 = *p_cutljinner; cut_lj14 = *p_cutlj; - cut_coulinv14 = 1/cut_coul14; - cut_lj_inner3inv = (1/cut_lj_inner14) * (1/cut_lj_inner14) * - (1/cut_lj_inner14); + cut_coulinv14 = 1 / cut_coul14; + cut_lj_inner3inv = (1 / cut_lj_inner14) * (1 / cut_lj_inner14) * (1 / cut_lj_inner14); 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; } @@ -432,11 +431,11 @@ void DihedralCharmmfsw::init_style() void DihedralCharmmfsw::write_restart(FILE *fp) { - fwrite(&k[1],sizeof(double),atom->ndihedraltypes,fp); - fwrite(&multiplicity[1],sizeof(int),atom->ndihedraltypes,fp); - fwrite(&shift[1],sizeof(int),atom->ndihedraltypes,fp); - fwrite(&weight[1],sizeof(double),atom->ndihedraltypes,fp); - fwrite(&weightflag,sizeof(int),1,fp); + fwrite(&k[1], sizeof(double), atom->ndihedraltypes, fp); + fwrite(&multiplicity[1], sizeof(int), atom->ndihedraltypes, fp); + fwrite(&shift[1], sizeof(int), atom->ndihedraltypes, fp); + fwrite(&weight[1], sizeof(double), atom->ndihedraltypes, fp); + fwrite(&weightflag, sizeof(int), 1, fp); } /* ---------------------------------------------------------------------- @@ -448,22 +447,22 @@ void DihedralCharmmfsw::read_restart(FILE *fp) allocate(); if (comm->me == 0) { - 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,&shift[1],sizeof(int),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, &k[1], sizeof(double), 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, &weight[1], sizeof(double), atom->ndihedraltypes, 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(&multiplicity[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(&weightflag,1,MPI_INT,0,world); + MPI_Bcast(&k[1], atom->ndihedraltypes, MPI_DOUBLE, 0, world); + MPI_Bcast(&multiplicity[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(&weightflag, 1, MPI_INT, 0, world); for (int i = 1; i <= atom->ndihedraltypes; i++) { setflag[i] = 1; - cos_shift[i] = cos(MY_PI*shift[i]/180.0); - sin_shift[i] = sin(MY_PI*shift[i]/180.0); + cos_shift[i] = cos(DEG2RAD * shift[i]); + sin_shift[i] = sin(DEG2RAD * shift[i]); } } @@ -474,6 +473,5 @@ void DihedralCharmmfsw::read_restart(FILE *fp) void DihedralCharmmfsw::write_data(FILE *fp) { 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]); } - diff --git a/src/MOLECULE/dihedral_harmonic.cpp b/src/MOLECULE/dihedral_harmonic.cpp index a582987de1..3a4d1c9df0 100644 --- a/src/MOLECULE/dihedral_harmonic.cpp +++ b/src/MOLECULE/dihedral_harmonic.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -29,12 +28,11 @@ using namespace LAMMPS_NS; -#define TOLERANCE 0.05 -#define SMALL 0.001 +static constexpr double TOLERANCE = 0.05; /* ---------------------------------------------------------------------- */ -DihedralHarmonic::DihedralHarmonic(LAMMPS *lmp) : Dihedral(lmp) +DihedralHarmonic::DihedralHarmonic(LAMMPS *_lmp) : Dihedral(_lmp) { writedata = 1; } @@ -57,16 +55,16 @@ DihedralHarmonic::~DihedralHarmonic() void DihedralHarmonic::compute(int eflag, int vflag) { - int i1,i2,i3,i4,i,m,n,type; - double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,vb2xm,vb2ym,vb2zm; - 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 df,df1,ddf1,fg,hg,fga,hgb,gaa,gbb; - double dtfx,dtfy,dtfz,dtgx,dtgy,dtgz,dthx,dthy,dthz; - double c,s,p,sx2,sy2,sz2; + int i1, i2, i3, i4, i, m, n, type; + double vb1x, vb1y, vb1z, vb2x, vb2y, vb2z, vb3x, vb3y, vb3z, vb2xm, vb2ym, vb2zm; + 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 df, df1, ddf1, fg, hg, fga, hgb, gaa, gbb; + double dtfx, dtfy, dtfz, dtgx, dtgy, dtgz, dthx, dthy, dthz; + double c, s, p, sx2, sy2, sz2; edihedral = 0.0; - ev_init(eflag,vflag); + ev_init(eflag, vflag); double **x = atom->x; double **f = atom->f; @@ -106,31 +104,30 @@ void DihedralHarmonic::compute(int eflag, int vflag) // c,s calculation - ax = vb1y*vb2zm - vb1z*vb2ym; - ay = vb1z*vb2xm - vb1x*vb2zm; - az = vb1x*vb2ym - vb1y*vb2xm; - bx = vb3y*vb2zm - vb3z*vb2ym; - by = vb3z*vb2xm - vb3x*vb2zm; - bz = vb3x*vb2ym - vb3y*vb2xm; + ax = vb1y * vb2zm - vb1z * vb2ym; + ay = vb1z * vb2xm - vb1x * vb2zm; + az = vb1x * vb2ym - vb1y * vb2xm; + bx = vb3y * vb2zm - vb3z * vb2ym; + by = vb3z * vb2xm - vb3x * vb2zm; + bz = vb3x * vb2ym - vb3y * vb2xm; - rasq = ax*ax + ay*ay + az*az; - rbsq = bx*bx + by*by + bz*bz; - rgsq = vb2xm*vb2xm + vb2ym*vb2ym + vb2zm*vb2zm; + rasq = ax * ax + ay * ay + az * az; + rbsq = bx * bx + by * by + bz * bz; + rgsq = vb2xm * vb2xm + vb2ym * vb2ym + vb2zm * vb2zm; rg = sqrt(rgsq); rginv = ra2inv = rb2inv = 0.0; - if (rg > 0) rginv = 1.0/rg; - if (rasq > 0) ra2inv = 1.0/rasq; - if (rbsq > 0) rb2inv = 1.0/rbsq; - rabinv = sqrt(ra2inv*rb2inv); + if (rg > 0) rginv = 1.0 / rg; + if (rasq > 0) ra2inv = 1.0 / rasq; + if (rbsq > 0) rb2inv = 1.0 / rbsq; + rabinv = sqrt(ra2inv * rb2inv); - c = (ax*bx + ay*by + az*bz)*rabinv; - s = rg*rabinv*(ax*vb3x + ay*vb3y + az*vb3z); + c = (ax * bx + ay * by + az * bz) * rabinv; + s = rg * rabinv * (ax * vb3x + ay * vb3y + az * vb3z); // error check - if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE)) - problem(FLERR, i1, i2, i3, i4); + if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE)) problem(FLERR, i1, i2, i3, i4); 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; for (i = 0; i < m; i++) { - ddf1 = p*c - df1*s; - df1 = p*s + df1*c; + ddf1 = p * c - df1 * s; + df1 = p * s + df1 * c; p = ddf1; } - p = p*cos_shift[type] + df1*sin_shift[type]; - df1 = df1*cos_shift[type] - ddf1*sin_shift[type]; + p = p * cos_shift[type] + df1 * sin_shift[type]; + df1 = df1 * cos_shift[type] - ddf1 * sin_shift[type]; df1 *= -m; p += 1.0; @@ -157,40 +154,40 @@ void DihedralHarmonic::compute(int eflag, int vflag) if (eflag) edihedral = k[type] * p; - fg = vb1x*vb2xm + vb1y*vb2ym + vb1z*vb2zm; - hg = vb3x*vb2xm + vb3y*vb2ym + vb3z*vb2zm; - fga = fg*ra2inv*rginv; - hgb = hg*rb2inv*rginv; - gaa = -ra2inv*rg; - gbb = rb2inv*rg; + fg = vb1x * vb2xm + vb1y * vb2ym + vb1z * vb2zm; + hg = vb3x * vb2xm + vb3y * vb2ym + vb3z * vb2zm; + fga = fg * ra2inv * rginv; + hgb = hg * rb2inv * rginv; + gaa = -ra2inv * rg; + gbb = rb2inv * rg; - dtfx = gaa*ax; - dtfy = gaa*ay; - dtfz = gaa*az; - dtgx = fga*ax - hgb*bx; - dtgy = fga*ay - hgb*by; - dtgz = fga*az - hgb*bz; - dthx = gbb*bx; - dthy = gbb*by; - dthz = gbb*bz; + dtfx = gaa * ax; + dtfy = gaa * ay; + dtfz = gaa * az; + dtgx = fga * ax - hgb * bx; + dtgy = fga * ay - hgb * by; + dtgz = fga * az - hgb * bz; + dthx = gbb * bx; + dthy = gbb * by; + dthz = gbb * bz; df = -k[type] * df1; - sx2 = df*dtgx; - sy2 = df*dtgy; - sz2 = df*dtgz; + sx2 = df * dtgx; + sy2 = df * dtgy; + sz2 = df * dtgz; - f1[0] = df*dtfx; - f1[1] = df*dtfy; - f1[2] = df*dtfz; + f1[0] = df * dtfx; + f1[1] = df * dtfy; + f1[2] = df * dtfz; f2[0] = sx2 - f1[0]; f2[1] = sy2 - f1[1]; f2[2] = sz2 - f1[2]; - f4[0] = df*dthx; - f4[1] = df*dthy; - f4[2] = df*dthz; + f4[0] = df * dthx; + f4[1] = df * dthy; + f4[2] = df * dthz; f3[0] = -sx2 - f4[0]; f3[1] = -sy2 - f4[1]; @@ -223,8 +220,8 @@ void DihedralHarmonic::compute(int eflag, int vflag) } if (evflag) - ev_tally(i1,i2,i3,i4,nlocal,newton_bond,edihedral,f1,f3,f4, - vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z); + ev_tally(i1, i2, i3, i4, nlocal, newton_bond, edihedral, f1, f3, f4, vb1x, vb1y, vb1z, vb2x, + vb2y, vb2z, vb3x, vb3y, vb3z); } } @@ -233,16 +230,16 @@ void DihedralHarmonic::compute(int eflag, int vflag) void DihedralHarmonic::allocate() { allocated = 1; - int n = atom->ndihedraltypes; + const int np1 = atom->ndihedraltypes + 1; - memory->create(k,n+1,"dihedral:k"); - memory->create(sign,n+1,"dihedral:sign"); - memory->create(multiplicity,n+1,"dihedral:multiplicity"); - memory->create(cos_shift,n+1,"dihedral:cos_shift"); - memory->create(sin_shift,n+1,"dihedral:sin_shift"); + memory->create(k, np1, "dihedral:k"); + memory->create(sign, np1, "dihedral:sign"); + memory->create(multiplicity, np1, "dihedral:multiplicity"); + memory->create(cos_shift, np1, "dihedral:cos_shift"); + memory->create(sin_shift, np1, "dihedral:sin_shift"); - memory->create(setflag,n+1,"dihedral:setflag"); - for (int i = 1; i <= n; i++) setflag[i] = 0; + memory->create(setflag, np1, "dihedral:setflag"); + for (int i = 1; i < np1; i++) setflag[i] = 0; } /* ---------------------------------------------------------------------- @@ -251,24 +248,24 @@ void DihedralHarmonic::allocate() 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(); - int ilo,ihi; - utils::bounds(FLERR,arg[0],1,atom->ndihedraltypes,ilo,ihi,error); + int ilo, ihi; + utils::bounds(FLERR, arg[0], 1, atom->ndihedraltypes, ilo, ihi, error); - double k_one = utils::numeric(FLERR,arg[1],false,lmp); - int sign_one = utils::inumeric(FLERR,arg[2],false,lmp); - int multiplicity_one = utils::inumeric(FLERR,arg[3],false,lmp); + double k_one = utils::numeric(FLERR, arg[1], false, lmp); + int sign_one = utils::inumeric(FLERR, arg[2], false, lmp); + int multiplicity_one = utils::inumeric(FLERR, arg[3], false, lmp); // require sign = +/- 1 for backwards compatibility // arbitrary phase angle shift could be allowed, but would break // backwards compatibility and is probably not needed 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) - error->all(FLERR,"Incorrect multiplicity arg for dihedral coefficients"); + error->all(FLERR, "Incorrect multiplicity arg for dihedral coefficients"); int count = 0; for (int i = ilo; i <= ihi; i++) { @@ -286,7 +283,7 @@ void DihedralHarmonic::coeff(int narg, char **arg) 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) { - fwrite(&k[1],sizeof(double),atom->ndihedraltypes,fp); - fwrite(&sign[1],sizeof(int),atom->ndihedraltypes,fp); - fwrite(&multiplicity[1],sizeof(int),atom->ndihedraltypes,fp); + fwrite(&k[1], sizeof(double), atom->ndihedraltypes, fp); + fwrite(&sign[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(); if (comm->me == 0) { - 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,&multiplicity[1],sizeof(int),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, &multiplicity[1], sizeof(int), atom->ndihedraltypes, fp, nullptr, error); } - MPI_Bcast(&k[1],atom->ndihedraltypes,MPI_DOUBLE,0,world); - MPI_Bcast(&sign[1],atom->ndihedraltypes,MPI_INT,0,world); - MPI_Bcast(&multiplicity[1],atom->ndihedraltypes,MPI_INT,0,world); + MPI_Bcast(&k[1], atom->ndihedraltypes, MPI_DOUBLE, 0, world); + MPI_Bcast(&sign[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++) { setflag[i] = 1; @@ -336,6 +333,5 @@ void DihedralHarmonic::read_restart(FILE *fp) void DihedralHarmonic::write_data(FILE *fp) { 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]); } - diff --git a/src/MOLECULE/dihedral_multi_harmonic.cpp b/src/MOLECULE/dihedral_multi_harmonic.cpp index 5527dd600c..b33a3c2320 100644 --- a/src/MOLECULE/dihedral_multi_harmonic.cpp +++ b/src/MOLECULE/dihedral_multi_harmonic.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -29,12 +28,12 @@ using namespace LAMMPS_NS; -#define TOLERANCE 0.05 -#define SMALL 0.001 +static constexpr double TOLERANCE = 0.05; +static constexpr double SMALL = 0.001; /* ---------------------------------------------------------------------- */ -DihedralMultiHarmonic::DihedralMultiHarmonic(LAMMPS *lmp) : Dihedral(lmp) +DihedralMultiHarmonic::DihedralMultiHarmonic(LAMMPS *_lmp) : Dihedral(_lmp) { writedata = 1; } @@ -57,17 +56,17 @@ DihedralMultiHarmonic::~DihedralMultiHarmonic() void DihedralMultiHarmonic::compute(int eflag, int vflag) { - int i1,i2,i3,i4,n,type; - double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,vb2xm,vb2ym,vb2zm; - double edihedral,f1[3],f2[3],f3[3],f4[3]; - double sb1,sb2,sb3,rb1,rb3,c0,b1mag2,b1mag,b2mag2; - double b2mag,b3mag2,b3mag,ctmp,r12c1,c1mag,r12c2; - double c2mag,sc1,sc2,s1,s12,c,p,pd,a,a11,a22; - double a33,a12,a13,a23,sx2,sy2,sz2; - double s2,sin2; + int i1, i2, i3, i4, n, type; + double vb1x, vb1y, vb1z, vb2x, vb2y, vb2z, vb3x, vb3y, vb3z, vb2xm, vb2ym, vb2zm; + double edihedral, f1[3], f2[3], f3[3], f4[3]; + double sb1, sb2, sb3, rb1, rb3, c0, b1mag2, b1mag, b2mag2; + double b2mag, b3mag2, b3mag, ctmp, r12c1, c1mag, r12c2; + double c2mag, sc1, sc2, s1, s12, c, p, pd, a, a11, a22; + double a33, a12, a13, a23, sx2, sy2, sz2; + double s2, sin2; edihedral = 0.0; - ev_init(eflag,vflag); + ev_init(eflag, vflag); double **x = atom->x; double **f = atom->f; @@ -107,53 +106,52 @@ void DihedralMultiHarmonic::compute(int eflag, int vflag) // c0 calculation - sb1 = 1.0 / (vb1x*vb1x + vb1y*vb1y + vb1z*vb1z); - sb2 = 1.0 / (vb2x*vb2x + vb2y*vb2y + vb2z*vb2z); - sb3 = 1.0 / (vb3x*vb3x + vb3y*vb3y + vb3z*vb3z); + sb1 = 1.0 / (vb1x * vb1x + vb1y * vb1y + vb1z * vb1z); + sb2 = 1.0 / (vb2x * vb2x + vb2y * vb2y + vb2z * vb2z); + sb3 = 1.0 / (vb3x * vb3x + vb3y * vb3y + vb3z * vb3z); rb1 = sqrt(sb1); 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 - b1mag2 = vb1x*vb1x + vb1y*vb1y + vb1z*vb1z; + b1mag2 = vb1x * vb1x + vb1y * vb1y + vb1z * vb1z; b1mag = sqrt(b1mag2); - b2mag2 = vb2x*vb2x + vb2y*vb2y + vb2z*vb2z; + b2mag2 = vb2x * vb2x + vb2y * vb2y + vb2z * vb2z; b2mag = sqrt(b2mag2); - b3mag2 = vb3x*vb3x + vb3y*vb3y + vb3z*vb3z; + b3mag2 = vb3x * vb3x + vb3y * vb3y + vb3z * vb3z; b3mag = sqrt(b3mag2); - ctmp = vb1x*vb2x + vb1y*vb2y + vb1z*vb2z; - r12c1 = 1.0 / (b1mag*b2mag); + ctmp = vb1x * vb2x + vb1y * vb2y + vb1z * vb2z; + r12c1 = 1.0 / (b1mag * b2mag); c1mag = ctmp * r12c1; - ctmp = vb2xm*vb3x + vb2ym*vb3y + vb2zm*vb3z; - r12c2 = 1.0 / (b2mag*b3mag); + ctmp = vb2xm * vb3x + vb2ym * vb3y + vb2zm * vb3z; + r12c2 = 1.0 / (b2mag * b3mag); c2mag = ctmp * r12c2; // 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); 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); if (sc2 < SMALL) sc2 = SMALL; - sc2 = 1.0/sc2; + sc2 = 1.0 / sc2; s1 = sc1 * sc1; s2 = sc2 * sc2; s12 = sc1 * sc2; - c = (c0 + c1mag*c2mag) * s12; + c = (c0 + c1mag * c2mag) * s12; // error check - if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE)) - problem(FLERR, i1, i2, i3, i4); + if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE)) problem(FLERR, i1, i2, i3, i4); 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) // pd = dp/dc - 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])); + 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])); if (eflag) edihedral = p; a = pd; c = c * a; s12 = s12 * a; - a11 = c*sb1*s1; - a22 = -sb2 * (2.0*c0*s12 - c*(s1+s2)); - a33 = c*sb3*s2; - a12 = -r12c1*(c1mag*c*s1 + c2mag*s12); - a13 = -rb1*rb3*s12; - a23 = r12c2*(c2mag*c*s2 + c1mag*s12); + a11 = c * sb1 * s1; + a22 = -sb2 * (2.0 * c0 * s12 - c * (s1 + s2)); + a33 = c * sb3 * s2; + a12 = -r12c1 * (c1mag * c * s1 + c2mag * s12); + a13 = -rb1 * rb3 * s12; + a23 = r12c2 * (c2mag * c * s2 + c1mag * s12); - sx2 = a12*vb1x + a22*vb2x + a23*vb3x; - sy2 = a12*vb1y + a22*vb2y + a23*vb3y; - sz2 = a12*vb1z + a22*vb2z + a23*vb3z; + sx2 = a12 * vb1x + a22 * vb2x + a23 * vb3x; + sy2 = a12 * vb1y + a22 * vb2y + a23 * vb3y; + sz2 = a12 * vb1z + a22 * vb2z + a23 * vb3z; - f1[0] = a11*vb1x + a12*vb2x + a13*vb3x; - f1[1] = a11*vb1y + a12*vb2y + a13*vb3y; - f1[2] = a11*vb1z + a12*vb2z + a13*vb3z; + f1[0] = a11 * vb1x + a12 * vb2x + a13 * vb3x; + f1[1] = a11 * vb1y + a12 * vb2y + a13 * vb3y; + f1[2] = a11 * vb1z + a12 * vb2z + a13 * vb3z; f2[0] = -sx2 - f1[0]; f2[1] = -sy2 - f1[1]; f2[2] = -sz2 - f1[2]; - f4[0] = a13*vb1x + a23*vb2x + a33*vb3x; - f4[1] = a13*vb1y + a23*vb2y + a33*vb3y; - f4[2] = a13*vb1z + a23*vb2z + a33*vb3z; + f4[0] = a13 * vb1x + a23 * vb2x + a33 * vb3x; + f4[1] = a13 * vb1y + a23 * vb2y + a33 * vb3y; + f4[2] = a13 * vb1z + a23 * vb2z + a33 * vb3z; f3[0] = sx2 - f4[0]; f3[1] = sy2 - f4[1]; @@ -224,8 +222,8 @@ void DihedralMultiHarmonic::compute(int eflag, int vflag) } if (evflag) - ev_tally(i1,i2,i3,i4,nlocal,newton_bond,edihedral,f1,f3,f4, - vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z); + ev_tally(i1, i2, i3, i4, nlocal, newton_bond, edihedral, f1, f3, f4, vb1x, vb1y, vb1z, vb2x, + vb2y, vb2z, vb3x, vb3y, vb3z); } } @@ -234,16 +232,16 @@ void DihedralMultiHarmonic::compute(int eflag, int vflag) void DihedralMultiHarmonic::allocate() { allocated = 1; - int n = atom->ndihedraltypes; + const int np1 = atom->ndihedraltypes + 1; - memory->create(a1,n+1,"dihedral:a1"); - memory->create(a2,n+1,"dihedral:a2"); - memory->create(a3,n+1,"dihedral:a3"); - memory->create(a4,n+1,"dihedral:a4"); - memory->create(a5,n+1,"dihedral:a5"); + memory->create(a1, np1, "dihedral:a1"); + memory->create(a2, np1, "dihedral:a2"); + memory->create(a3, np1, "dihedral:a3"); + memory->create(a4, np1, "dihedral:a4"); + memory->create(a5, np1, "dihedral:a5"); - memory->create(setflag,n+1,"dihedral:setflag"); - for (int i = 1; i <= n; i++) setflag[i] = 0; + memory->create(setflag, np1, "dihedral:setflag"); + for (int i = 1; i < np1; i++) setflag[i] = 0; } /* ---------------------------------------------------------------------- @@ -252,17 +250,17 @@ void DihedralMultiHarmonic::allocate() 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(); - int ilo,ihi; - utils::bounds(FLERR,arg[0],1,atom->ndihedraltypes,ilo,ihi,error); + int ilo, ihi; + utils::bounds(FLERR, arg[0], 1, atom->ndihedraltypes, ilo, ihi, error); - double a1_one = utils::numeric(FLERR,arg[1],false,lmp); - double a2_one = utils::numeric(FLERR,arg[2],false,lmp); - double a3_one = utils::numeric(FLERR,arg[3],false,lmp); - double a4_one = utils::numeric(FLERR,arg[4],false,lmp); - double a5_one = utils::numeric(FLERR,arg[5],false,lmp); + double a1_one = utils::numeric(FLERR, arg[1], false, lmp); + double a2_one = utils::numeric(FLERR, arg[2], false, lmp); + double a3_one = utils::numeric(FLERR, arg[3], false, lmp); + double a4_one = utils::numeric(FLERR, arg[4], false, lmp); + double a5_one = utils::numeric(FLERR, arg[5], false, lmp); int count = 0; for (int i = ilo; i <= ihi; i++) { @@ -275,7 +273,7 @@ void DihedralMultiHarmonic::coeff(int narg, char **arg) 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) { - fwrite(&a1[1],sizeof(double),atom->ndihedraltypes,fp); - fwrite(&a2[1],sizeof(double),atom->ndihedraltypes,fp); - fwrite(&a3[1],sizeof(double),atom->ndihedraltypes,fp); - fwrite(&a4[1],sizeof(double),atom->ndihedraltypes,fp); - fwrite(&a5[1],sizeof(double),atom->ndihedraltypes,fp); + fwrite(&a1[1], sizeof(double), atom->ndihedraltypes, fp); + fwrite(&a2[1], sizeof(double), atom->ndihedraltypes, fp); + fwrite(&a3[1], sizeof(double), atom->ndihedraltypes, fp); + fwrite(&a4[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(); if (comm->me == 0) { - 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,&a3[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, &a1[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, &a4[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(&a2[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(&a5[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(&a3[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); 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) { 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]); } diff --git a/src/MOLECULE/dihedral_opls.cpp b/src/MOLECULE/dihedral_opls.cpp index 2f99019ef7..7c60cba923 100644 --- a/src/MOLECULE/dihedral_opls.cpp +++ b/src/MOLECULE/dihedral_opls.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -29,13 +28,13 @@ using namespace LAMMPS_NS; -#define TOLERANCE 0.05 -#define SMALL 0.001 -#define SMALLER 0.00001 +static constexpr double TOLERANCE = 0.05; +static constexpr double SMALL = 0.001; +static constexpr double SMALLER = 0.00001; /* ---------------------------------------------------------------------- */ -DihedralOPLS::DihedralOPLS(LAMMPS *lmp) : Dihedral(lmp) +DihedralOPLS::DihedralOPLS(LAMMPS *_lmp) : Dihedral(_lmp) { writedata = 1; } @@ -57,17 +56,17 @@ DihedralOPLS::~DihedralOPLS() void DihedralOPLS::compute(int eflag, int vflag) { - int i1,i2,i3,i4,n,type; - double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,vb2xm,vb2ym,vb2zm; - double edihedral,f1[3],f2[3],f3[3],f4[3]; - double sb1,sb2,sb3,rb1,rb3,c0,b1mag2,b1mag,b2mag2; - double b2mag,b3mag2,b3mag,ctmp,r12c1,c1mag,r12c2; - double c2mag,sc1,sc2,s1,s12,c,p,pd,a,a11,a22; - double a33,a12,a13,a23,sx2,sy2,sz2; - double s2,cx,cy,cz,cmag,dx,phi,si,siinv,sin2; + int i1, i2, i3, i4, n, type; + double vb1x, vb1y, vb1z, vb2x, vb2y, vb2z, vb3x, vb3y, vb3z, vb2xm, vb2ym, vb2zm; + double edihedral, f1[3], f2[3], f3[3], f4[3]; + double sb1, sb2, sb3, rb1, rb3, c0, b1mag2, b1mag, b2mag2; + double b2mag, b3mag2, b3mag, ctmp, r12c1, c1mag, r12c2; + double c2mag, sc1, sc2, s1, s12, c, p, pd, a, a11, a22; + double a33, a12, a13, a23, sx2, sy2, sz2; + double s2, cx, cy, cz, cmag, dx, phi, si, siinv, sin2; edihedral = 0.0; - ev_init(eflag,vflag); + ev_init(eflag, vflag); double **x = atom->x; double **f = atom->f; @@ -107,59 +106,58 @@ void DihedralOPLS::compute(int eflag, int vflag) // c0 calculation - sb1 = 1.0 / (vb1x*vb1x + vb1y*vb1y + vb1z*vb1z); - sb2 = 1.0 / (vb2x*vb2x + vb2y*vb2y + vb2z*vb2z); - sb3 = 1.0 / (vb3x*vb3x + vb3y*vb3y + vb3z*vb3z); + sb1 = 1.0 / (vb1x * vb1x + vb1y * vb1y + vb1z * vb1z); + sb2 = 1.0 / (vb2x * vb2x + vb2y * vb2y + vb2z * vb2z); + sb3 = 1.0 / (vb3x * vb3x + vb3y * vb3y + vb3z * vb3z); rb1 = sqrt(sb1); 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 - b1mag2 = vb1x*vb1x + vb1y*vb1y + vb1z*vb1z; + b1mag2 = vb1x * vb1x + vb1y * vb1y + vb1z * vb1z; b1mag = sqrt(b1mag2); - b2mag2 = vb2x*vb2x + vb2y*vb2y + vb2z*vb2z; + b2mag2 = vb2x * vb2x + vb2y * vb2y + vb2z * vb2z; b2mag = sqrt(b2mag2); - b3mag2 = vb3x*vb3x + vb3y*vb3y + vb3z*vb3z; + b3mag2 = vb3x * vb3x + vb3y * vb3y + vb3z * vb3z; b3mag = sqrt(b3mag2); - ctmp = vb1x*vb2x + vb1y*vb2y + vb1z*vb2z; - r12c1 = 1.0 / (b1mag*b2mag); + ctmp = vb1x * vb2x + vb1y * vb2y + vb1z * vb2z; + r12c1 = 1.0 / (b1mag * b2mag); c1mag = ctmp * r12c1; - ctmp = vb2xm*vb3x + vb2ym*vb3y + vb2zm*vb3z; - r12c2 = 1.0 / (b2mag*b3mag); + ctmp = vb2xm * vb3x + vb2ym * vb3y + vb2zm * vb3z; + r12c2 = 1.0 / (b2mag * b3mag); c2mag = ctmp * r12c2; // 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); 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); if (sc2 < SMALL) sc2 = SMALL; - sc2 = 1.0/sc2; + sc2 = 1.0 / sc2; s1 = sc1 * sc1; s2 = sc2 * sc2; s12 = sc1 * sc2; - c = (c0 + c1mag*c2mag) * s12; + c = (c0 + c1mag * c2mag) * s12; - cx = vb1y*vb2z - vb1z*vb2y; - cy = vb1z*vb2x - vb1x*vb2z; - cz = vb1x*vb2y - vb1y*vb2x; - cmag = sqrt(cx*cx + cy*cy + cz*cz); - dx = (cx*vb3x + cy*vb3y + cz*vb3z)/cmag/b3mag; + cx = vb1y * vb2z - vb1z * vb2y; + cy = vb1z * vb2x - vb1x * vb2z; + cz = vb1x * vb2y - vb1y * vb2x; + cmag = sqrt(cx * cx + cy * cy + cz * cz); + dx = (cx * vb3x + cy * vb3y + cz * vb3z) / cmag / b3mag; // error check - if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE)) - problem(FLERR, i1, i2, i3, i4); + if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE)) problem(FLERR, i1, i2, i3, i4); 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; si = sin(phi); 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)) + - 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 + - 3.0*k3[type]*sin(3.0*phi)*siinv - 4.0*k4[type]*sin(4.0*phi)*siinv; + 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)); + 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; if (eflag) edihedral = p; a = pd; c = c * a; s12 = s12 * a; - a11 = c*sb1*s1; - a22 = -sb2 * (2.0*c0*s12 - c*(s1+s2)); - a33 = c*sb3*s2; - a12 = -r12c1 * (c1mag*c*s1 + c2mag*s12); - a13 = -rb1*rb3*s12; - a23 = r12c2 * (c2mag*c*s2 + c1mag*s12); + a11 = c * sb1 * s1; + a22 = -sb2 * (2.0 * c0 * s12 - c * (s1 + s2)); + a33 = c * sb3 * s2; + a12 = -r12c1 * (c1mag * c * s1 + c2mag * s12); + a13 = -rb1 * rb3 * s12; + a23 = r12c2 * (c2mag * c * s2 + c1mag * s12); - sx2 = a12*vb1x + a22*vb2x + a23*vb3x; - sy2 = a12*vb1y + a22*vb2y + a23*vb3y; - sz2 = a12*vb1z + a22*vb2z + a23*vb3z; + sx2 = a12 * vb1x + a22 * vb2x + a23 * vb3x; + sy2 = a12 * vb1y + a22 * vb2y + a23 * vb3y; + sz2 = a12 * vb1z + a22 * vb2z + a23 * vb3z; - f1[0] = a11*vb1x + a12*vb2x + a13*vb3x; - f1[1] = a11*vb1y + a12*vb2y + a13*vb3y; - f1[2] = a11*vb1z + a12*vb2z + a13*vb3z; + f1[0] = a11 * vb1x + a12 * vb2x + a13 * vb3x; + f1[1] = a11 * vb1y + a12 * vb2y + a13 * vb3y; + f1[2] = a11 * vb1z + a12 * vb2z + a13 * vb3z; f2[0] = -sx2 - f1[0]; f2[1] = -sy2 - f1[1]; f2[2] = -sz2 - f1[2]; - f4[0] = a13*vb1x + a23*vb2x + a33*vb3x; - f4[1] = a13*vb1y + a23*vb2y + a33*vb3y; - f4[2] = a13*vb1z + a23*vb2z + a33*vb3z; + f4[0] = a13 * vb1x + a23 * vb2x + a33 * vb3x; + f4[1] = a13 * vb1y + a23 * vb2y + a33 * vb3y; + f4[2] = a13 * vb1z + a23 * vb2z + a33 * vb3z; f3[0] = sx2 - f4[0]; f3[1] = sy2 - f4[1]; @@ -238,8 +236,8 @@ void DihedralOPLS::compute(int eflag, int vflag) } if (evflag) - ev_tally(i1,i2,i3,i4,nlocal,newton_bond,edihedral,f1,f3,f4, - vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z); + ev_tally(i1, i2, i3, i4, nlocal, newton_bond, edihedral, f1, f3, f4, vb1x, vb1y, vb1z, vb2x, + vb2y, vb2z, vb3x, vb3y, vb3z); } } @@ -248,15 +246,15 @@ void DihedralOPLS::compute(int eflag, int vflag) void DihedralOPLS::allocate() { allocated = 1; - int n = atom->ndihedraltypes; + const int np1 = atom->ndihedraltypes + 1; - memory->create(k1,n+1,"dihedral:k1"); - memory->create(k2,n+1,"dihedral:k2"); - memory->create(k3,n+1,"dihedral:k3"); - memory->create(k4,n+1,"dihedral:k4"); + memory->create(k1, np1, "dihedral:k1"); + memory->create(k2, np1, "dihedral:k2"); + memory->create(k3, np1, "dihedral:k3"); + memory->create(k4, np1, "dihedral:k4"); - memory->create(setflag,n+1,"dihedral:setflag"); - for (int i = 1; i <= n; i++) setflag[i] = 0; + memory->create(setflag, np1, "dihedral:setflag"); + for (int i = 1; i < np1; i++) setflag[i] = 0; } /* ---------------------------------------------------------------------- @@ -265,30 +263,30 @@ void DihedralOPLS::allocate() 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(); - int ilo,ihi; - utils::bounds(FLERR,arg[0],1,atom->ndihedraltypes,ilo,ihi,error); + int ilo, ihi; + utils::bounds(FLERR, arg[0], 1, atom->ndihedraltypes, ilo, ihi, error); - double k1_one = utils::numeric(FLERR,arg[1],false,lmp); - double k2_one = utils::numeric(FLERR,arg[2],false,lmp); - double k3_one = utils::numeric(FLERR,arg[3],false,lmp); - double k4_one = utils::numeric(FLERR,arg[4],false,lmp); + double k1_one = utils::numeric(FLERR, arg[1], false, lmp); + double k2_one = utils::numeric(FLERR, arg[2], false, lmp); + double k3_one = utils::numeric(FLERR, arg[3], false, lmp); + double k4_one = utils::numeric(FLERR, arg[4], false, lmp); // store 1/2 factor with prefactor int count = 0; for (int i = ilo; i <= ihi; i++) { - k1[i] = 0.5*k1_one; - k2[i] = 0.5*k2_one; - k3[i] = 0.5*k3_one; - k4[i] = 0.5*k4_one; + k1[i] = 0.5 * k1_one; + k2[i] = 0.5 * k2_one; + k3[i] = 0.5 * k3_one; + k4[i] = 0.5 * k4_one; setflag[i] = 1; 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) { - fwrite(&k1[1],sizeof(double),atom->ndihedraltypes,fp); - fwrite(&k2[1],sizeof(double),atom->ndihedraltypes,fp); - fwrite(&k3[1],sizeof(double),atom->ndihedraltypes,fp); - fwrite(&k4[1],sizeof(double),atom->ndihedraltypes,fp); + fwrite(&k1[1], sizeof(double), atom->ndihedraltypes, fp); + fwrite(&k2[1], sizeof(double), atom->ndihedraltypes, fp); + fwrite(&k3[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(); if (comm->me == 0) { - 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,&k3[1],sizeof(double),atom->ndihedraltypes,fp,nullptr,error); - utils::sfread(FLERR,&k4[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, &k3[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(&k2[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(&k1[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(&k4[1], atom->ndihedraltypes, MPI_DOUBLE, 0, world); 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) { 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]); } -