From ebe5e6c02492e8ea3eae4607e6c6096b5a3da462 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 29 Jul 2021 14:44:15 -0400 Subject: [PATCH] reformat --- src/INTERLAYER/pair_coul_shield.cpp | 223 ++++++++++++++-------------- 1 file changed, 115 insertions(+), 108 deletions(-) diff --git a/src/INTERLAYER/pair_coul_shield.cpp b/src/INTERLAYER/pair_coul_shield.cpp index 2cda85ead6..dbb12abf2e 100644 --- a/src/INTERLAYER/pair_coul_shield.cpp +++ b/src/INTERLAYER/pair_coul_shield.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -20,22 +19,24 @@ ------------------------------------------------------------------------- */ #include "pair_coul_shield.h" -#include + #include "atom.h" #include "comm.h" -#include "force.h" -#include "neighbor.h" -#include "neigh_list.h" -#include "memory.h" -#include "math_special.h" #include "error.h" +#include "force.h" +#include "math_special.h" +#include "memory.h" +#include "neigh_list.h" +#include "neighbor.h" +#include using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -PairCoulShield::PairCoulShield(LAMMPS *lmp) : Pair(lmp) { +PairCoulShield::PairCoulShield(LAMMPS *lmp) : Pair(lmp) +{ tap_flag = 1; } @@ -57,13 +58,13 @@ PairCoulShield::~PairCoulShield() void PairCoulShield::compute(int eflag, int vflag) { - int i,j,ii,jj,inum,jnum,itype,jtype; - double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,ecoul,fpair,Tap,dTap; - double rsq,r,r3,rarg,th,depsdr,epsr,forcecoul,factor_coul,Vc,fvc; - int *ilist,*jlist,*numneigh,**firstneigh; + int i, j, ii, jj, inum, jnum, itype, jtype; + double qtmp, xtmp, ytmp, ztmp, delx, dely, delz, ecoul, fpair, Tap, dTap; + double rsq, r, r3, rarg, th, depsdr, epsr, forcecoul, factor_coul, Vc, fvc; + int *ilist, *jlist, *numneigh, **firstneigh; ecoul = 0.0; - ev_init(eflag,vflag); + ev_init(eflag, vflag); double **x = atom->x; double **f = atom->f; @@ -100,47 +101,51 @@ void PairCoulShield::compute(int eflag, int vflag) delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; + rsq = delx * delx + dely * dely + delz * delz; jtype = type[j]; // only include the interaction between different layers if (rsq < cutsq[itype][jtype] && atom->molecule[i] != atom->molecule[j]) { r = sqrt(rsq); - r3 = rsq*r; - rarg = 1.0/sigmae[itype][jtype]; + r3 = rsq * r; + rarg = 1.0 / sigmae[itype][jtype]; th = r3 + MathSpecial::cube(rarg); - epsr = 1.0/pow(th,0.333333333333333333333333); + epsr = 1.0 / pow(th, 1.0 / 3.0); depsdr = MathSpecial::square(epsr); depsdr *= depsdr; - Vc = qqrd2e*qtmp*q[j]*epsr; + Vc = qqrd2e * qtmp * q[j] * epsr; // turn on/off taper function if (tap_flag) { - Tap = calc_Tap(r,cut[itype][jtype]); - dTap = calc_dTap(r,cut[itype][jtype]); - } else {Tap = 1.0; dTap = 0.0;} + Tap = calc_Tap(r, cut[itype][jtype]); + dTap = calc_dTap(r, cut[itype][jtype]); + } else { + Tap = 1.0; + dTap = 0.0; + } - forcecoul = qqrd2e*qtmp*q[j]*r*depsdr; - fvc = forcecoul*Tap - Vc*dTap/r; - fpair = factor_coul*fvc; + forcecoul = qqrd2e * qtmp * q[j] * r * depsdr; + fvc = forcecoul * Tap - Vc * dTap / r; + fpair = factor_coul * fvc; - f[i][0] += delx*fpair; - f[i][1] += dely*fpair; - f[i][2] += delz*fpair; + f[i][0] += delx * fpair; + f[i][1] += dely * fpair; + f[i][2] += delz * fpair; if (newton_pair || j < nlocal) { - f[j][0] -= delx*fpair; - f[j][1] -= dely*fpair; - f[j][2] -= delz*fpair; + f[j][0] -= delx * fpair; + f[j][1] -= dely * fpair; + f[j][2] -= delz * fpair; } if (eflag) { - if (tap_flag) ecoul = Vc*Tap; - else ecoul = Vc - offset[itype][jtype]; + if (tap_flag) + ecoul = Vc * Tap; + else + ecoul = Vc - offset[itype][jtype]; ecoul *= factor_coul; } - if (evflag) ev_tally(i,j,nlocal,newton_pair,0.0, - ecoul,fpair,delx,dely,delz); + if (evflag) ev_tally(i, j, nlocal, newton_pair, 0.0, ecoul, fpair, delx, dely, delz); } } } @@ -157,15 +162,14 @@ void PairCoulShield::allocate() allocated = 1; int n = atom->ntypes; - memory->create(setflag,n+1,n+1,"pair:setflag"); + memory->create(setflag, n + 1, n + 1, "pair:setflag"); for (int i = 1; i <= n; i++) - for (int j = i; j <= n; j++) - setflag[i][j] = 0; + for (int j = i; j <= n; j++) setflag[i][j] = 0; - memory->create(cutsq,n+1,n+1,"pair:cutsq"); - memory->create(cut,n+1,n+1,"pair:cut"); - memory->create(sigmae,n+1,n+1,"pair:sigmae"); - memory->create(offset,n+1,n+1,"pair:offset"); + memory->create(cutsq, n + 1, n + 1, "pair:cutsq"); + memory->create(cut, n + 1, n + 1, "pair:cut"); + memory->create(sigmae, n + 1, n + 1, "pair:sigmae"); + memory->create(offset, n + 1, n + 1, "pair:offset"); } /* ---------------------------------------------------------------------- @@ -174,17 +178,17 @@ void PairCoulShield::allocate() void PairCoulShield::settings(int narg, char **arg) { - if (narg < 1 || narg > 2) error->all(FLERR,"Illegal pair_style command"); + if (narg < 1 || narg > 2) error->all(FLERR, "Illegal pair_style command"); - cut_global = utils::numeric(FLERR,arg[0],false,lmp); - if (narg == 2) tap_flag = utils::numeric(FLERR,arg[1],false,lmp); + cut_global = utils::numeric(FLERR, arg[0], false, lmp); + if (narg == 2) tap_flag = utils::numeric(FLERR, arg[1], false, lmp); // reset cutoffs that have been explicitly set if (allocated) { - int i,j; + int i, j; for (i = 1; i <= atom->ntypes; i++) - for (j = i+1; j <= atom->ntypes; j++) + for (j = i + 1; j <= atom->ntypes; j++) if (setflag[i][j]) cut[i][j] = cut_global; } } @@ -195,21 +199,21 @@ void PairCoulShield::settings(int narg, char **arg) void PairCoulShield::coeff(int narg, char **arg) { - if (narg < 3 || narg > 4) error->all(FLERR,"Incorrect args for pair coefficients"); + if (narg < 3 || narg > 4) error->all(FLERR, "Incorrect args for pair coefficients"); if (!allocated) allocate(); - int ilo,ihi,jlo,jhi; - utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error); - utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error); + int ilo, ihi, jlo, jhi; + utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error); + utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error); - double sigmae_one = utils::numeric(FLERR,arg[2],false,lmp); + double sigmae_one = utils::numeric(FLERR, arg[2], false, lmp); double cut_one = cut_global; - if (narg == 4) cut_one = utils::numeric(FLERR,arg[3],false,lmp); + if (narg == 4) cut_one = utils::numeric(FLERR, arg[3], false, lmp); int count = 0; for (int i = ilo; i <= ihi; i++) { - for (int j = MAX(jlo,i); j <= jhi; j++) { + for (int j = MAX(jlo, i); j <= jhi; j++) { sigmae[i][j] = sigmae_one; cut[i][j] = cut_one; setflag[i][j] = 1; @@ -217,22 +221,20 @@ void PairCoulShield::coeff(int narg, char **arg) } } - if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); + if (count == 0) error->all(FLERR, "Incorrect args for pair coefficients"); } - /* ---------------------------------------------------------------------- init specific to this pair style ------------------------------------------------------------------------- */ void PairCoulShield::init_style() { - if (!atom->q_flag) - error->all(FLERR,"Pair style coul/shield requires atom attribute q"); + if (!atom->q_flag) error->all(FLERR, "Pair style coul/shield requires atom attribute q"); if (!atom->molecule_flag) - error->all(FLERR,"Pair style coul/shield requires atom attribute molecule"); + error->all(FLERR, "Pair style coul/shield requires atom attribute molecule"); - neighbor->request(this,instance_me); + neighbor->request(this, instance_me); } /* ---------------------------------------------------------------------- @@ -242,22 +244,23 @@ void PairCoulShield::init_style() double PairCoulShield::init_one(int i, int j) { if (setflag[i][j] == 0) { - error->all(FLERR,"for pair style coul/shield, parameters need to be set explicitly for all pairs."); + error->all(FLERR, + "for pair style coul/shield, parameters need to be set explicitly for all pairs."); } double *q = atom->q; double qqrd2e = force->qqrd2e; - double r,r3,rarg,th,epsr; + double r, r3, rarg, th, epsr; if (offset_flag) { - r = cut[i][j]; - r3 = r*r*r; - rarg = 1.0/sigmae[i][j]; - th = r3 + MathSpecial::cube(rarg); - epsr = 1.0/pow(th,0.333333333333333333); - offset[i][j] = qqrd2e*q[i]*q[j]*epsr; - } else offset[i][j] = 0.0; - + r = cut[i][j]; + r3 = r * r * r; + rarg = 1.0 / sigmae[i][j]; + th = r3 + MathSpecial::cube(rarg); + epsr = 1.0 / pow(th, 1.0/3.0); + offset[i][j] = qqrd2e * q[i] * q[j] * epsr; + } else + offset[i][j] = 0.0; sigmae[j][i] = sigmae[i][j]; offset[j][i] = offset[i][j]; @@ -274,13 +277,13 @@ void PairCoulShield::write_restart(FILE *fp) { write_restart_settings(fp); - int i,j; + int i, j; for (i = 1; i <= atom->ntypes; i++) for (j = i; j <= atom->ntypes; j++) { - fwrite(&setflag[i][j],sizeof(int),1,fp); + fwrite(&setflag[i][j], sizeof(int), 1, fp); if (setflag[i][j]) { - fwrite(&sigmae[i][j],sizeof(double),1,fp); - fwrite(&cut[i][j],sizeof(double),1,fp); + fwrite(&sigmae[i][j], sizeof(double), 1, fp); + fwrite(&cut[i][j], sizeof(double), 1, fp); } } } @@ -294,19 +297,19 @@ void PairCoulShield::read_restart(FILE *fp) read_restart_settings(fp); allocate(); - int i,j; + int i, j; int me = comm->me; for (i = 1; i <= atom->ntypes; i++) for (j = i; j <= atom->ntypes; j++) { - if (me == 0) utils::sfread(FLERR,&setflag[i][j],sizeof(int),1,fp,nullptr,error); - MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); + if (me == 0) utils::sfread(FLERR, &setflag[i][j], sizeof(int), 1, fp, nullptr, error); + MPI_Bcast(&setflag[i][j], 1, MPI_INT, 0, world); if (setflag[i][j]) { if (me == 0) { - utils::sfread(FLERR,&sigmae[i][j],sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&cut[i][j],sizeof(double),1,fp,nullptr,error); + utils::sfread(FLERR, &sigmae[i][j], sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &cut[i][j], sizeof(double), 1, fp, nullptr, error); } - MPI_Bcast(&sigmae[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&sigmae[i][j], 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut[i][j], 1, MPI_DOUBLE, 0, world); } } } @@ -317,9 +320,9 @@ void PairCoulShield::read_restart(FILE *fp) void PairCoulShield::write_restart_settings(FILE *fp) { - fwrite(&cut_global,sizeof(double),1,fp); - fwrite(&offset_flag,sizeof(int),1,fp); - fwrite(&mix_flag,sizeof(int),1,fp); + fwrite(&cut_global, sizeof(double), 1, fp); + fwrite(&offset_flag, sizeof(int), 1, fp); + fwrite(&mix_flag, sizeof(int), 1, fp); } /* ---------------------------------------------------------------------- @@ -329,23 +332,22 @@ void PairCoulShield::write_restart_settings(FILE *fp) void PairCoulShield::read_restart_settings(FILE *fp) { if (comm->me == 0) { - utils::sfread(FLERR,&cut_global,sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&offset_flag,sizeof(int),1,fp,nullptr,error); - utils::sfread(FLERR,&mix_flag,sizeof(int),1,fp,nullptr,error); + utils::sfread(FLERR, &cut_global, sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &offset_flag, sizeof(int), 1, fp, nullptr, error); + utils::sfread(FLERR, &mix_flag, sizeof(int), 1, fp, nullptr, error); } - MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world); - MPI_Bcast(&offset_flag,1,MPI_INT,0,world); - MPI_Bcast(&mix_flag,1,MPI_INT,0,world); + MPI_Bcast(&cut_global, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&offset_flag, 1, MPI_INT, 0, world); + MPI_Bcast(&mix_flag, 1, MPI_INT, 0, world); } /* ---------------------------------------------------------------------- */ -double PairCoulShield::single(int i, int j, int itype, int jtype, - double rsq, double factor_coul, double /*factor_lj*/, - double &fforce) +double PairCoulShield::single(int i, int j, int itype, int jtype, double rsq, double factor_coul, + double /*factor_lj*/, double &fforce) { - double r, rarg,Vc,fvc,forcecoul,phishieldec; - double r3,th,epsr,depsdr,Tap,dTap; + double r, rarg, Vc, fvc, forcecoul, phishieldec; + double r3, th, epsr, depsdr, Tap, dTap; double *q = atom->q; double qqrd2e = force->qqrd2e; @@ -357,25 +359,30 @@ double PairCoulShield::single(int i, int j, int itype, int jtype, } r = sqrt(rsq); - r3 = rsq*r; - rarg = 1.0/sigmae[itype][jtype]; + r3 = rsq * r; + rarg = 1.0 / sigmae[itype][jtype]; th = r3 + MathSpecial::cube(rarg); - epsr = 1.0/pow(th,0.333333333333333333); - depsdr = epsr*epsr; + epsr = 1.0 / pow(th, 1.0/3.0); + depsdr = epsr * epsr; depsdr *= depsdr; - Vc = qqrd2e*q[i]*q[j]*epsr; + Vc = qqrd2e * q[i] * q[j] * epsr; // turn on/off taper function if (tap_flag) { - Tap = calc_Tap(r,cut[itype][jtype]); - dTap = calc_dTap(r,cut[itype][jtype]); - } else {Tap = 1.0; dTap = 0.0;} + Tap = calc_Tap(r, cut[itype][jtype]); + dTap = calc_dTap(r, cut[itype][jtype]); + } else { + Tap = 1.0; + dTap = 0.0; + } - forcecoul = qqrd2e*q[i]*q[j]*r*depsdr; - fvc = forcecoul*Tap - Vc*dTap/r; - fforce = factor_coul*fvc; + forcecoul = qqrd2e * q[i] * q[j] * r * depsdr; + fvc = forcecoul * Tap - Vc * dTap / r; + fforce = factor_coul * fvc; - if (tap_flag) phishieldec = Vc*Tap; - else phishieldec = Vc - offset[itype][jtype]; - return factor_coul*phishieldec; + if (tap_flag) + phishieldec = Vc * Tap; + else + phishieldec = Vc - offset[itype][jtype]; + return factor_coul * phishieldec; }