diff --git a/src/CG-SDK/pair_lj_sdk_coul_long.cpp b/src/CG-SDK/pair_lj_sdk_coul_long.cpp index 0fd0d2d30a..e1d79789b8 100644 --- a/src/CG-SDK/pair_lj_sdk_coul_long.cpp +++ b/src/CG-SDK/pair_lj_sdk_coul_long.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -19,17 +18,17 @@ #include "pair_lj_sdk_coul_long.h" -#include -#include #include "atom.h" #include "comm.h" +#include "error.h" #include "force.h" #include "kspace.h" -#include "neighbor.h" -#include "neigh_list.h" #include "memory.h" -#include "error.h" +#include "neigh_list.h" +#include "neighbor.h" +#include +#include #define LMP_NEED_SDK_FIND_LJ_TYPE 1 #include "lj_sdk_common.h" @@ -37,13 +36,13 @@ using namespace LAMMPS_NS; using namespace LJSDKParms; -#define EWALD_F 1.12837917 -#define EWALD_P 0.3275911 -#define A1 0.254829592 -#define A2 -0.284496736 -#define A3 1.421413741 -#define A4 -1.453152027 -#define A5 1.061405429 +#define EWALD_F 1.12837917 +#define EWALD_P 0.3275911 +#define A1 0.254829592 +#define A2 -0.284496736 +#define A3 1.421413741 +#define A4 -1.453152027 +#define A5 1.061405429 /* ---------------------------------------------------------------------- */ @@ -86,19 +85,25 @@ PairLJSDKCoulLong::~PairLJSDKCoulLong() void PairLJSDKCoulLong::compute(int eflag, int vflag) { - ev_init(eflag,vflag); + ev_init(eflag, vflag); if (evflag) { if (eflag) { - if (force->newton_pair) eval<1,1,1>(); - else eval<1,1,0>(); + if (force->newton_pair) + eval<1, 1, 1>(); + else + eval<1, 1, 0>(); } else { - if (force->newton_pair) eval<1,0,1>(); - else eval<1,0,0>(); + if (force->newton_pair) + eval<1, 0, 1>(); + else + eval<1, 0, 0>(); } } else { - if (force->newton_pair) eval<0,0,1>(); - else eval<0,0,0>(); + if (force->newton_pair) + eval<0, 0, 1>(); + else + eval<0, 0, 0>(); } if (vflag_fdotr) virial_fdotr_compute(); @@ -106,29 +111,28 @@ void PairLJSDKCoulLong::compute(int eflag, int vflag) /* ---------------------------------------------------------------------- */ -template -void PairLJSDKCoulLong::eval() +template void PairLJSDKCoulLong::eval() { - int i,ii,j,jj,jtype,itable; - double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair; - double fraction,table; - double r,rsq,r2inv,forcecoul,forcelj,factor_coul,factor_lj; - double grij,expm2,prefactor,t,erfc; + int i, ii, j, jj, jtype, itable; + double qtmp, xtmp, ytmp, ztmp, delx, dely, delz, evdwl, ecoul, fpair; + double fraction, table; + double r, rsq, r2inv, forcecoul, forcelj, factor_coul, factor_lj; + double grij, expm2, prefactor, t, erfc; - const double * const * const x = atom->x; - double * const * const f = atom->f; - const double * const q = atom->q; - const int * const type = atom->type; + const double *const *const x = atom->x; + double *const *const f = atom->f; + const double *const q = atom->q; + const int *const type = atom->type; const int nlocal = atom->nlocal; - const double * const special_coul = force->special_coul; - const double * const special_lj = force->special_lj; + const double *const special_coul = force->special_coul; + const double *const special_lj = force->special_lj; const double qqrd2e = force->qqrd2e; - double fxtmp,fytmp,fztmp; + double fxtmp, fytmp, fztmp; const int inum = list->inum; - const int * const ilist = list->ilist; - const int * const numneigh = list->numneigh; - const int * const * const firstneigh = list->firstneigh; + const int *const ilist = list->ilist; + const int *const numneigh = list->numneigh; + const int *const *const firstneigh = list->firstneigh; // loop over neighbors of my atoms @@ -139,10 +143,10 @@ void PairLJSDKCoulLong::eval() xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; - fxtmp=fytmp=fztmp=0.0; + fxtmp = fytmp = fztmp = 0.0; const int itype = type[i]; - const int * const jlist = firstneigh[i]; + const int *const jlist = firstneigh[i]; const int jnum = numneigh[i]; for (jj = 0; jj < jnum; jj++) { @@ -156,26 +160,26 @@ void PairLJSDKCoulLong::eval() 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]; if (rsq < cutsq[itype][jtype]) { - r2inv = 1.0/rsq; + r2inv = 1.0 / rsq; const int ljt = lj_type[itype][jtype]; if (rsq < cut_coulsq) { if (!ncoultablebits || rsq <= tabinnersq) { r = sqrt(rsq); grij = g_ewald * r; - expm2 = exp(-grij*grij); - t = 1.0 / (1.0 + EWALD_P*grij); - erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2; - prefactor = qqrd2e * qtmp*q[j]/r; - forcecoul = prefactor * (erfc + EWALD_F*grij*expm2); - if (EFLAG) ecoul = prefactor*erfc; + expm2 = exp(-grij * grij); + t = 1.0 / (1.0 + EWALD_P * grij); + erfc = t * (A1 + t * (A2 + t * (A3 + t * (A4 + t * A5)))) * expm2; + prefactor = qqrd2e * qtmp * q[j] / r; + forcecoul = prefactor * (erfc + EWALD_F * grij * expm2); + if (EFLAG) ecoul = prefactor * erfc; if (factor_coul < 1.0) { - forcecoul -= (1.0-factor_coul)*prefactor; - if (EFLAG) ecoul -= (1.0-factor_coul)*prefactor; + forcecoul -= (1.0 - factor_coul) * prefactor; + if (EFLAG) ecoul -= (1.0 - factor_coul) * prefactor; } } else { union_int_float_t rsq_lookup; @@ -183,15 +187,14 @@ void PairLJSDKCoulLong::eval() itable = rsq_lookup.i & ncoulmask; itable >>= ncoulshiftbits; fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable]; - table = ftable[itable] + fraction*dftable[itable]; - forcecoul = qtmp*q[j] * table; - if (EFLAG) ecoul = qtmp*q[j] * - (etable[itable] + fraction*detable[itable]); + table = ftable[itable] + fraction * dftable[itable]; + forcecoul = qtmp * q[j] * table; + if (EFLAG) ecoul = qtmp * q[j] * (etable[itable] + fraction * detable[itable]); if (factor_coul < 1.0) { - table = ctable[itable] + fraction*dctable[itable]; - prefactor = qtmp*q[j] * table; - forcecoul -= (1.0-factor_coul)*prefactor; - if (EFLAG) ecoul -= (1.0-factor_coul)*prefactor; + table = ctable[itable] + fraction * dctable[itable]; + prefactor = qtmp * q[j] * table; + forcecoul -= (1.0 - factor_coul) * prefactor; + if (EFLAG) ecoul -= (1.0 - factor_coul) * prefactor; } } } @@ -199,30 +202,27 @@ void PairLJSDKCoulLong::eval() if (rsq < cut_ljsq[itype][jtype]) { if (ljt == LJ12_4) { - const double r4inv=r2inv*r2inv; - forcelj = r4inv*(lj1[itype][jtype]*r4inv*r4inv - - lj2[itype][jtype]); + const double r4inv = r2inv * r2inv; + forcelj = r4inv * (lj1[itype][jtype] * r4inv * r4inv - lj2[itype][jtype]); if (EFLAG) - evdwl = r4inv*(lj3[itype][jtype]*r4inv*r4inv - - lj4[itype][jtype]) - offset[itype][jtype]; + evdwl = r4inv * (lj3[itype][jtype] * r4inv * r4inv - lj4[itype][jtype]) - + offset[itype][jtype]; } else if (ljt == LJ9_6) { - const double r3inv = r2inv*sqrt(r2inv); - const double r6inv = r3inv*r3inv; - forcelj = r6inv*(lj1[itype][jtype]*r3inv - - lj2[itype][jtype]); + const double r3inv = r2inv * sqrt(r2inv); + const double r6inv = r3inv * r3inv; + forcelj = r6inv * (lj1[itype][jtype] * r3inv - lj2[itype][jtype]); if (EFLAG) - evdwl = r6inv*(lj3[itype][jtype]*r3inv - - lj4[itype][jtype]) - offset[itype][jtype]; + evdwl = + r6inv * (lj3[itype][jtype] * r3inv - lj4[itype][jtype]) - offset[itype][jtype]; } else if (ljt == LJ12_6) { - const double r6inv = r2inv*r2inv*r2inv; - forcelj = r6inv*(lj1[itype][jtype]*r6inv - - lj2[itype][jtype]); + const double r6inv = r2inv * r2inv * r2inv; + forcelj = r6inv * (lj1[itype][jtype] * r6inv - lj2[itype][jtype]); if (EFLAG) - evdwl = r6inv*(lj3[itype][jtype]*r6inv - - lj4[itype][jtype]) - offset[itype][jtype]; + evdwl = + r6inv * (lj3[itype][jtype] * r6inv - lj4[itype][jtype]) - offset[itype][jtype]; } forcelj *= factor_lj; if (EFLAG) evdwl *= factor_lj; @@ -230,17 +230,16 @@ void PairLJSDKCoulLong::eval() fpair = (forcecoul + forcelj) * r2inv; - fxtmp += delx*fpair; - fytmp += dely*fpair; - fztmp += delz*fpair; + fxtmp += delx * fpair; + fytmp += dely * fpair; + fztmp += 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 (EVFLAG) ev_tally(i,j,nlocal,NEWTON_PAIR, - evdwl,ecoul,fpair,delx,dely,delz); + if (EVFLAG) ev_tally(i, j, nlocal, NEWTON_PAIR, evdwl, ecoul, fpair, delx, dely, delz); } } f[i][0] += fxtmp; @@ -249,7 +248,6 @@ void PairLJSDKCoulLong::eval() } } - /* ---------------------------------------------------------------------- allocate all arrays ------------------------------------------------------------------------- */ @@ -257,33 +255,33 @@ void PairLJSDKCoulLong::eval() void PairLJSDKCoulLong::allocate() { allocated = 1; - int n = atom->ntypes; + int np1 = atom->ntypes + 1; - memory->create(setflag,n+1,n+1,"pair:setflag"); - memory->create(lj_type,n+1,n+1,"pair:lj_type"); - for (int i = 1; i <= n; i++) { - for (int j = i; j <= n; j++) { + memory->create(setflag, np1, np1, "pair:setflag"); + memory->create(lj_type, np1, np1, "pair:lj_type"); + for (int i = 1; i < np1; i++) { + for (int j = i; j < np1; j++) { setflag[i][j] = 0; lj_type[i][j] = LJ_NOT_SET; } } - memory->create(cutsq,n+1,n+1,"pair:cutsq"); + memory->create(cutsq, np1, np1, "pair:cutsq"); - memory->create(cut_lj,n+1,n+1,"pair:cut_lj"); - memory->create(cut_ljsq,n+1,n+1,"pair:cut_ljsq"); - memory->create(epsilon,n+1,n+1,"pair:epsilon"); - memory->create(sigma,n+1,n+1,"pair:sigma"); + memory->create(cut_lj, np1, np1, "pair:cut_lj"); + memory->create(cut_ljsq, np1, np1, "pair:cut_ljsq"); + memory->create(epsilon, np1, np1, "pair:epsilon"); + memory->create(sigma, np1, np1, "pair:sigma"); - memory->create(lj1,n+1,n+1,"pair:lj1"); - memory->create(lj2,n+1,n+1,"pair:lj2"); - memory->create(lj3,n+1,n+1,"pair:lj3"); - memory->create(lj4,n+1,n+1,"pair:lj4"); + memory->create(lj1, np1, np1, "pair:lj1"); + memory->create(lj2, np1, np1, "pair:lj2"); + memory->create(lj3, np1, np1, "pair:lj3"); + memory->create(lj4, np1, np1, "pair:lj4"); - memory->create(offset,n+1,n+1,"pair:offset"); + memory->create(offset, np1, np1, "pair:offset"); - memory->create(rminsq,n+1,n+1,"pair:rminsq"); - memory->create(emin,n+1,n+1,"pair:emin"); + memory->create(rminsq, np1, np1, "pair:rminsq"); + memory->create(emin, np1, np1, "pair:emin"); } /* ---------------------------------------------------------------------- @@ -292,16 +290,18 @@ void PairLJSDKCoulLong::allocate() void PairLJSDKCoulLong::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_lj_global = utils::numeric(FLERR,arg[0],false,lmp); - if (narg == 1) cut_coul = cut_lj_global; - else cut_coul = utils::numeric(FLERR,arg[1],false,lmp); + cut_lj_global = utils::numeric(FLERR, arg[0], false, lmp); + if (narg == 1) + cut_coul = cut_lj_global; + else + cut_coul = 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; j <= atom->ntypes; j++) if (setflag[i][j]) cut_lj[i][j] = cut_lj_global; @@ -314,27 +314,25 @@ void PairLJSDKCoulLong::settings(int narg, char **arg) void PairLJSDKCoulLong::coeff(int narg, char **arg) { - if (narg < 5 || narg > 6) - error->all(FLERR,"Incorrect args for pair coefficients"); + if (narg < 5 || narg > 6) 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); - int lj_type_one = find_lj_type(arg[2],lj_type_list); - if (lj_type_one == LJ_NOT_SET) - error->all(FLERR,"Cannot parse LJ type flag."); + int lj_type_one = find_lj_type(arg[2], lj_type_list); + if (lj_type_one == LJ_NOT_SET) error->all(FLERR, "Cannot parse LJ type flag."); - double epsilon_one = utils::numeric(FLERR,arg[3],false,lmp); - double sigma_one = utils::numeric(FLERR,arg[4],false,lmp); + double epsilon_one = utils::numeric(FLERR, arg[3], false, lmp); + double sigma_one = utils::numeric(FLERR, arg[4], false, lmp); double cut_lj_one = cut_lj_global; - if (narg == 6) cut_lj_one = utils::numeric(FLERR,arg[5],false,lmp); + if (narg == 6) cut_lj_one = utils::numeric(FLERR, arg[5], 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++) { lj_type[i][j] = lj_type_one; epsilon[i][j] = epsilon_one; sigma[i][j] = sigma_one; @@ -344,7 +342,7 @@ void PairLJSDKCoulLong::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"); } /* ---------------------------------------------------------------------- @@ -353,22 +351,20 @@ void PairLJSDKCoulLong::coeff(int narg, char **arg) void PairLJSDKCoulLong::init_style() { - if (!atom->q_flag) - error->all(FLERR,"Pair style lj/cut/coul/long requires atom attribute q"); + if (!atom->q_flag) error->all(FLERR, "Pair style lj/cut/coul/long requires atom attribute q"); - neighbor->request(this,instance_me); + neighbor->add_request(this); cut_coulsq = cut_coul * cut_coul; // insure use of KSpace long-range solver, set g_ewald - if (force->kspace == nullptr) - error->all(FLERR,"Pair style requires a KSpace style"); + if (force->kspace == nullptr) error->all(FLERR, "Pair style requires a KSpace style"); g_ewald = force->kspace->g_ewald; // setup force tables (no rRESPA support yet) - if (ncoultablebits) init_tables(cut_coul,nullptr); + if (ncoultablebits) init_tables(cut_coul, nullptr); } /* ---------------------------------------------------------------------- @@ -378,29 +374,28 @@ void PairLJSDKCoulLong::init_style() double PairLJSDKCoulLong::init_one(int i, int j) { if (setflag[i][j] == 0) - error->all(FLERR,"No mixing support for lj/sdk/coul/long. " + error->all(FLERR, + "No mixing support for lj/sdk/coul/long. " "Coefficients for all pairs need to be set explicitly."); const int ljt = lj_type[i][j]; - if (ljt == LJ_NOT_SET) - error->all(FLERR,"unrecognized LJ parameter flag"); + if (ljt == LJ_NOT_SET) error->all(FLERR, "unrecognized LJ parameter flag"); - double cut = MAX(cut_lj[i][j],cut_coul); + double cut = MAX(cut_lj[i][j], cut_coul); cut_ljsq[i][j] = cut_lj[i][j] * cut_lj[i][j]; - lj1[i][j] = lj_prefact[ljt] * lj_pow1[ljt] * epsilon[i][j] * - pow(sigma[i][j],lj_pow1[ljt]); - lj2[i][j] = lj_prefact[ljt] * lj_pow2[ljt] * epsilon[i][j] * - pow(sigma[i][j],lj_pow2[ljt]); - lj3[i][j] = lj_prefact[ljt] * epsilon[i][j] * pow(sigma[i][j],lj_pow1[ljt]); - lj4[i][j] = lj_prefact[ljt] * epsilon[i][j] * pow(sigma[i][j],lj_pow2[ljt]); + lj1[i][j] = lj_prefact[ljt] * lj_pow1[ljt] * epsilon[i][j] * pow(sigma[i][j], lj_pow1[ljt]); + lj2[i][j] = lj_prefact[ljt] * lj_pow2[ljt] * epsilon[i][j] * pow(sigma[i][j], lj_pow2[ljt]); + lj3[i][j] = lj_prefact[ljt] * epsilon[i][j] * pow(sigma[i][j], lj_pow1[ljt]); + lj4[i][j] = lj_prefact[ljt] * epsilon[i][j] * pow(sigma[i][j], lj_pow2[ljt]); if (offset_flag && (cut_lj[i][j] > 0.0)) { double ratio = sigma[i][j] / cut_lj[i][j]; - offset[i][j] = lj_prefact[ljt] * epsilon[i][j] * - (pow(ratio,lj_pow1[ljt]) - pow(ratio,lj_pow2[ljt])); - } else offset[i][j] = 0.0; + offset[i][j] = + lj_prefact[ljt] * epsilon[i][j] * (pow(ratio, lj_pow1[ljt]) - pow(ratio, lj_pow2[ljt])); + } else + offset[i][j] = 0.0; cut_ljsq[j][i] = cut_ljsq[i][j]; cut_lj[j][i] = cut_lj[i][j]; @@ -415,20 +410,19 @@ double PairLJSDKCoulLong::init_one(int i, int j) const double eps = epsilon[i][j]; const double sig = sigma[i][j]; - const double rmin = sig*exp(1.0/(lj_pow1[ljt]-lj_pow2[ljt]) - *log(lj_pow1[ljt]/lj_pow2[ljt]) ); - rminsq[j][i] = rminsq[i][j] = rmin*rmin; + const double rmin = + sig * exp(1.0 / (lj_pow1[ljt] - lj_pow2[ljt]) * log(lj_pow1[ljt] / lj_pow2[ljt])); + rminsq[j][i] = rminsq[i][j] = rmin * rmin; - const double ratio = sig/rmin; - const double emin_one = lj_prefact[ljt] * eps * (pow(ratio,lj_pow1[ljt]) - - pow(ratio,lj_pow2[ljt])); + const double ratio = sig / rmin; + const double emin_one = + lj_prefact[ljt] * eps * (pow(ratio, lj_pow1[ljt]) - pow(ratio, lj_pow2[ljt])); emin[j][i] = emin[i][j] = emin_one; // compute I,J contribution to long-range tail correction // count total # of atoms of type I and J via Allreduce - if (tail_flag) - error->all(FLERR,"Tail flag not supported by lj/sdk/coul/long pair style"); + if (tail_flag) error->all(FLERR, "Tail flag not supported by lj/sdk/coul/long pair style"); return cut; } @@ -441,15 +435,15 @@ void PairLJSDKCoulLong::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(&lj_type[i][j],sizeof(int),1,fp); - fwrite(&epsilon[i][j],sizeof(double),1,fp); - fwrite(&sigma[i][j],sizeof(double),1,fp); - fwrite(&cut_lj[i][j],sizeof(double),1,fp); + fwrite(&lj_type[i][j], sizeof(int), 1, fp); + fwrite(&epsilon[i][j], sizeof(double), 1, fp); + fwrite(&sigma[i][j], sizeof(double), 1, fp); + fwrite(&cut_lj[i][j], sizeof(double), 1, fp); } } } @@ -463,23 +457,23 @@ void PairLJSDKCoulLong::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,&lj_type[i][j],sizeof(int),1,fp,nullptr,error); - utils::sfread(FLERR,&epsilon[i][j],sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&sigma[i][j],sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&cut_lj[i][j],sizeof(double),1,fp,nullptr,error); + utils::sfread(FLERR, &lj_type[i][j], sizeof(int), 1, fp, nullptr, error); + utils::sfread(FLERR, &epsilon[i][j], sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &sigma[i][j], sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &cut_lj[i][j], sizeof(double), 1, fp, nullptr, error); } - MPI_Bcast(&lj_type[i][j],1,MPI_INT,0,world); - MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&cut_lj[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&lj_type[i][j], 1, MPI_INT, 0, world); + MPI_Bcast(&epsilon[i][j], 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&sigma[i][j], 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut_lj[i][j], 1, MPI_DOUBLE, 0, world); } } } @@ -490,13 +484,13 @@ void PairLJSDKCoulLong::read_restart(FILE *fp) void PairLJSDKCoulLong::write_restart_settings(FILE *fp) { - fwrite(&cut_lj_global,sizeof(double),1,fp); - fwrite(&cut_coul,sizeof(double),1,fp); - fwrite(&offset_flag,sizeof(int),1,fp); - fwrite(&mix_flag,sizeof(int),1,fp); - fwrite(&tail_flag,sizeof(int),1,fp); - fwrite(&ncoultablebits,sizeof(int),1,fp); - fwrite(&tabinner,sizeof(double),1,fp); + fwrite(&cut_lj_global, sizeof(double), 1, fp); + fwrite(&cut_coul, sizeof(double), 1, fp); + fwrite(&offset_flag, sizeof(int), 1, fp); + fwrite(&mix_flag, sizeof(int), 1, fp); + fwrite(&tail_flag, sizeof(int), 1, fp); + fwrite(&ncoultablebits, sizeof(int), 1, fp); + fwrite(&tabinner, sizeof(double), 1, fp); } /* ---------------------------------------------------------------------- @@ -506,21 +500,21 @@ void PairLJSDKCoulLong::write_restart_settings(FILE *fp) void PairLJSDKCoulLong::read_restart_settings(FILE *fp) { if (comm->me == 0) { - utils::sfread(FLERR,&cut_lj_global,sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&cut_coul,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,&tail_flag,sizeof(int),1,fp,nullptr,error); - utils::sfread(FLERR,&ncoultablebits,sizeof(int),1,fp,nullptr,error); - utils::sfread(FLERR,&tabinner,sizeof(double),1,fp,nullptr,error); + utils::sfread(FLERR, &cut_lj_global, sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &cut_coul, 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, &tail_flag, sizeof(int), 1, fp, nullptr, error); + utils::sfread(FLERR, &ncoultablebits, sizeof(int), 1, fp, nullptr, error); + utils::sfread(FLERR, &tabinner, sizeof(double), 1, fp, nullptr, error); } - MPI_Bcast(&cut_lj_global,1,MPI_DOUBLE,0,world); - MPI_Bcast(&cut_coul,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(&tail_flag,1,MPI_INT,0,world); - MPI_Bcast(&ncoultablebits,1,MPI_INT,0,world); - MPI_Bcast(&tabinner,1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut_lj_global, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut_coul, 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(&tail_flag, 1, MPI_INT, 0, world); + MPI_Bcast(&ncoultablebits, 1, MPI_INT, 0, world); + MPI_Bcast(&tabinner, 1, MPI_DOUBLE, 0, world); } /* ---------------------------------------------------------------------- @@ -529,7 +523,8 @@ void PairLJSDKCoulLong::read_restart_settings(FILE *fp) void PairLJSDKCoulLong::write_data(FILE *) { - error->one(FLERR, "Pair style lj/sdk/coul/* requires using " + error->one(FLERR, + "Pair style lj/sdk/coul/* requires using " "write_data with the 'pair ij' option"); } @@ -541,37 +536,35 @@ void PairLJSDKCoulLong::write_data_all(FILE *fp) { for (int i = 1; i <= atom->ntypes; i++) for (int j = i; j <= atom->ntypes; j++) - fprintf(fp,"%d %d %s %g %g %g\n",i,j,lj_type_list[lj_type[i][j]], - epsilon[i][j],sigma[i][j],cut_lj[i][j]); + fprintf(fp, "%d %d %s %g %g %g\n", i, j, lj_type_list[lj_type[i][j]], epsilon[i][j], + sigma[i][j], cut_lj[i][j]); } /* ---------------------------------------------------------------------- */ -double PairLJSDKCoulLong::single(int i, int j, int itype, int jtype, - double rsq, - double factor_coul, double factor_lj, - double &fforce) +double PairLJSDKCoulLong::single(int i, int j, int itype, int jtype, double rsq, double factor_coul, + double factor_lj, double &fforce) { - double r2inv,r,grij,expm2,t,erfc,prefactor; - double fraction,table,forcecoul,forcelj,phicoul,philj; + double r2inv, r, grij, expm2, t, erfc, prefactor; + double fraction, table, forcecoul, forcelj, phicoul, philj; int itable; forcecoul = forcelj = phicoul = philj = 0.0; - r2inv = 1.0/rsq; + r2inv = 1.0 / rsq; if (rsq < cut_coulsq) { if (!ncoultablebits || rsq <= tabinnersq) { r = sqrt(rsq); grij = g_ewald * r; - expm2 = exp(-grij*grij); - t = 1.0 / (1.0 + EWALD_P*grij); - erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2; - prefactor = force->qqrd2e * atom->q[i]*atom->q[j]/r; - forcecoul = prefactor * (erfc + EWALD_F*grij*expm2); - phicoul = prefactor*erfc; + expm2 = exp(-grij * grij); + t = 1.0 / (1.0 + EWALD_P * grij); + erfc = t * (A1 + t * (A2 + t * (A3 + t * (A4 + t * A5)))) * expm2; + prefactor = force->qqrd2e * atom->q[i] * atom->q[j] / r; + forcecoul = prefactor * (erfc + EWALD_F * grij * expm2); + phicoul = prefactor * erfc; if (factor_coul < 1.0) { - forcecoul -= (1.0-factor_coul)*prefactor; - phicoul -= (1.0-factor_coul)*prefactor; + forcecoul -= (1.0 - factor_coul) * prefactor; + phicoul -= (1.0 - factor_coul) * prefactor; } } else { union_int_float_t rsq_lookup_single; @@ -579,15 +572,15 @@ double PairLJSDKCoulLong::single(int i, int j, int itype, int jtype, itable = rsq_lookup_single.i & ncoulmask; itable >>= ncoulshiftbits; fraction = (rsq_lookup_single.f - rtable[itable]) * drtable[itable]; - table = ftable[itable] + fraction*dftable[itable]; - forcecoul = atom->q[i]*atom->q[j] * table; - table = etable[itable] + fraction*detable[itable]; - phicoul = atom->q[i]*atom->q[j] * table; + table = ftable[itable] + fraction * dftable[itable]; + forcecoul = atom->q[i] * atom->q[j] * table; + table = etable[itable] + fraction * detable[itable]; + phicoul = atom->q[i] * atom->q[j] * table; if (factor_coul < 1.0) { - table = ctable[itable] + fraction*dctable[itable]; - prefactor = atom->q[i]*atom->q[j] * table; - forcecoul -= (1.0-factor_coul)*prefactor; - phicoul -= (1.0-factor_coul)*prefactor; + table = ctable[itable] + fraction * dctable[itable]; + prefactor = atom->q[i] * atom->q[j] * table; + forcecoul -= (1.0 - factor_coul) * prefactor; + phicoul -= (1.0 - factor_coul) * prefactor; } } } @@ -596,27 +589,22 @@ double PairLJSDKCoulLong::single(int i, int j, int itype, int jtype, const int ljt = lj_type[itype][jtype]; if (ljt == LJ12_4) { - const double r4inv=r2inv*r2inv; - forcelj = r4inv*(lj1[itype][jtype]*r4inv*r4inv - - lj2[itype][jtype]); + const double r4inv = r2inv * r2inv; + forcelj = r4inv * (lj1[itype][jtype] * r4inv * r4inv - lj2[itype][jtype]); - philj = r4inv*(lj3[itype][jtype]*r4inv*r4inv - - lj4[itype][jtype]) - offset[itype][jtype]; + philj = + r4inv * (lj3[itype][jtype] * r4inv * r4inv - lj4[itype][jtype]) - offset[itype][jtype]; } else if (ljt == LJ9_6) { - const double r3inv = r2inv*sqrt(r2inv); - const double r6inv = r3inv*r3inv; - forcelj = r6inv*(lj1[itype][jtype]*r3inv - - lj2[itype][jtype]); - philj = r6inv*(lj3[itype][jtype]*r3inv - - lj4[itype][jtype]) - offset[itype][jtype]; + const double r3inv = r2inv * sqrt(r2inv); + const double r6inv = r3inv * r3inv; + forcelj = r6inv * (lj1[itype][jtype] * r3inv - lj2[itype][jtype]); + philj = r6inv * (lj3[itype][jtype] * r3inv - lj4[itype][jtype]) - offset[itype][jtype]; } else if (ljt == LJ12_6) { - const double r6inv = r2inv*r2inv*r2inv; - forcelj = r6inv*(lj1[itype][jtype]*r6inv - - lj2[itype][jtype]); - philj = r6inv*(lj3[itype][jtype]*r6inv - - lj4[itype][jtype]) - offset[itype][jtype]; + const double r6inv = r2inv * r2inv * r2inv; + forcelj = r6inv * (lj1[itype][jtype] * r6inv - lj2[itype][jtype]); + philj = r6inv * (lj3[itype][jtype] * r6inv - lj4[itype][jtype]) - offset[itype][jtype]; } forcelj *= factor_lj; philj *= factor_lj; @@ -632,18 +620,18 @@ double PairLJSDKCoulLong::single(int i, int j, int itype, int jtype, void *PairLJSDKCoulLong::extract(const char *str, int &dim) { dim = 2; - if (strcmp(str,"epsilon") == 0) return (void *) epsilon; - if (strcmp(str,"sigma") == 0) return (void *) sigma; - if (strcmp(str,"lj_type") == 0) return (void *) lj_type; - if (strcmp(str,"lj1") == 0) return (void *) lj1; - if (strcmp(str,"lj2") == 0) return (void *) lj2; - if (strcmp(str,"lj3") == 0) return (void *) lj3; - if (strcmp(str,"lj4") == 0) return (void *) lj4; - if (strcmp(str,"rminsq") == 0) return (void *) rminsq; - if (strcmp(str,"emin") == 0) return (void *) emin; + if (strcmp(str, "epsilon") == 0) return (void *) epsilon; + if (strcmp(str, "sigma") == 0) return (void *) sigma; + if (strcmp(str, "lj_type") == 0) return (void *) lj_type; + if (strcmp(str, "lj1") == 0) return (void *) lj1; + if (strcmp(str, "lj2") == 0) return (void *) lj2; + if (strcmp(str, "lj3") == 0) return (void *) lj3; + if (strcmp(str, "lj4") == 0) return (void *) lj4; + if (strcmp(str, "rminsq") == 0) return (void *) rminsq; + if (strcmp(str, "emin") == 0) return (void *) emin; dim = 0; - if (strcmp(str,"cut_coul") == 0) return (void *) &cut_coul; + if (strcmp(str, "cut_coul") == 0) return (void *) &cut_coul; return nullptr; } @@ -655,13 +643,13 @@ double PairLJSDKCoulLong::memory_usage() int n = atom->ntypes; // setflag/lj_type - bytes += (double)2 * (n+1)*(n+1)*sizeof(int); + bytes += (double) 2 * (n + 1) * (n + 1) * sizeof(int); // lj_cut/lj_cutsq/epsilon/sigma/offset/lj1/lj2/lj3/lj4/rminsq/emin - bytes += (double)11 * (n+1)*(n+1)*sizeof(double); + bytes += (double) 11 * (n + 1) * (n + 1) * sizeof(double); if (ncoultablebits) { - int ntable = 1<x; double **f = atom->f; @@ -104,34 +102,32 @@ void PairLJClass2::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]; if (rsq < cutsq[itype][jtype]) { - r2inv = 1.0/rsq; + r2inv = 1.0 / rsq; rinv = sqrt(r2inv); - r3inv = r2inv*rinv; - r6inv = r3inv*r3inv; - forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]); - fpair = factor_lj*forcelj*r2inv; + r3inv = r2inv * rinv; + r6inv = r3inv * r3inv; + forcelj = r6inv * (lj1[itype][jtype] * r3inv - lj2[itype][jtype]); + fpair = factor_lj * forcelj * r2inv; - 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) { - evdwl = r6inv*(lj3[itype][jtype]*r3inv-lj4[itype][jtype]) - - offset[itype][jtype]; + evdwl = r6inv * (lj3[itype][jtype] * r3inv - lj4[itype][jtype]) - offset[itype][jtype]; evdwl *= factor_lj; } - if (evflag) ev_tally(i,j,nlocal,newton_pair, - evdwl,0.0,fpair,delx,dely,delz); + if (evflag) ev_tally(i, j, nlocal, newton_pair, evdwl, 0.0, fpair, delx, dely, delz); } } } @@ -139,15 +135,14 @@ void PairLJClass2::compute(int eflag, int vflag) if (vflag_fdotr) virial_fdotr_compute(); } -/* ---------------------------------------------------------------------- -*/ +/* ---------------------------------------------------------------------- */ void PairLJClass2::compute_inner() { - int i,j,ii,jj,inum,jnum,itype,jtype; - double xtmp,ytmp,ztmp,delx,dely,delz,fpair; - double rsq,rinv,r2inv,r3inv,r6inv,forcelj,factor_lj,rsw; - int *ilist,*jlist,*numneigh,**firstneigh; + int i, j, ii, jj, inum, jnum, itype, jtype; + double xtmp, ytmp, ztmp, delx, dely, delz, fpair; + double rsq, rinv, r2inv, r3inv, r6inv, forcelj, factor_lj, rsw; + int *ilist, *jlist, *numneigh, **firstneigh; double **x = atom->x; double **f = atom->f; @@ -165,8 +160,8 @@ void PairLJClass2::compute_inner() double cut_out_off = cut_respa[1]; double cut_out_diff = cut_out_off - cut_out_on; - double cut_out_on_sq = cut_out_on*cut_out_on; - double cut_out_off_sq = cut_out_off*cut_out_off; + double cut_out_on_sq = cut_out_on * cut_out_on; + double cut_out_off_sq = cut_out_off * cut_out_off; // loop over neighbors of my atoms @@ -187,28 +182,28 @@ void PairLJClass2::compute_inner() 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; if (rsq < cut_out_off_sq) { - r2inv = 1.0/rsq; - rinv = sqrt(r2inv); - r3inv = r2inv*rinv; - r6inv = r3inv*r3inv; + r2inv = 1.0 / rsq; + rinv = sqrt(r2inv); + r3inv = r2inv * rinv; + r6inv = r3inv * r3inv; jtype = type[j]; - forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]); - fpair = factor_lj*forcelj*r2inv; + forcelj = r6inv * (lj1[itype][jtype] * r3inv - lj2[itype][jtype]); + fpair = factor_lj * forcelj * r2inv; if (rsq > cut_out_on_sq) { - rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff; - fpair *= 1.0 - rsw*rsw*(3.0 - 2.0*rsw); + rsw = (sqrt(rsq) - cut_out_on) / cut_out_diff; + fpair *= 1.0 - rsw * rsw * (3.0 - 2.0 * rsw); } - 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; } } } @@ -219,10 +214,10 @@ void PairLJClass2::compute_inner() void PairLJClass2::compute_middle() { - int i,j,ii,jj,inum,jnum,itype,jtype; - double xtmp,ytmp,ztmp,delx,dely,delz,fpair; - double rsq,rinv,r2inv,r3inv,r6inv,forcelj,factor_lj,rsw; - int *ilist,*jlist,*numneigh,**firstneigh; + int i, j, ii, jj, inum, jnum, itype, jtype; + double xtmp, ytmp, ztmp, delx, dely, delz, fpair; + double rsq, rinv, r2inv, r3inv, r6inv, forcelj, factor_lj, rsw; + int *ilist, *jlist, *numneigh, **firstneigh; double **x = atom->x; double **f = atom->f; @@ -243,10 +238,10 @@ void PairLJClass2::compute_middle() double cut_in_diff = cut_in_on - cut_in_off; double cut_out_diff = cut_out_off - cut_out_on; - double cut_in_off_sq = cut_in_off*cut_in_off; - double cut_in_on_sq = cut_in_on*cut_in_on; - double cut_out_on_sq = cut_out_on*cut_out_on; - double cut_out_off_sq = cut_out_off*cut_out_off; + double cut_in_off_sq = cut_in_off * cut_in_off; + double cut_in_on_sq = cut_in_on * cut_in_on; + double cut_out_on_sq = cut_out_on * cut_out_on; + double cut_out_off_sq = cut_out_off * cut_out_off; // loop over neighbors of my atoms @@ -267,32 +262,32 @@ void PairLJClass2::compute_middle() 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; if (rsq < cut_out_off_sq && rsq > cut_in_off_sq) { - r2inv = 1.0/rsq; - rinv = sqrt(r2inv); - r3inv = r2inv*rinv; - r6inv = r3inv*r3inv; + r2inv = 1.0 / rsq; + rinv = sqrt(r2inv); + r3inv = r2inv * rinv; + r6inv = r3inv * r3inv; jtype = type[j]; - forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]); - fpair = factor_lj*forcelj*r2inv; + forcelj = r6inv * (lj1[itype][jtype] * r3inv - lj2[itype][jtype]); + fpair = factor_lj * forcelj * r2inv; if (rsq < cut_in_on_sq) { - rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff; - fpair *= rsw*rsw*(3.0 - 2.0*rsw); + rsw = (sqrt(rsq) - cut_in_off) / cut_in_diff; + fpair *= rsw * rsw * (3.0 - 2.0 * rsw); } if (rsq > cut_out_on_sq) { - rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff; - fpair *= 1.0 + rsw*rsw*(2.0*rsw - 3.0); + rsw = (sqrt(rsq) - cut_out_on) / cut_out_diff; + fpair *= 1.0 + rsw * rsw * (2.0 * rsw - 3.0); } - 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; } } } @@ -303,13 +298,13 @@ void PairLJClass2::compute_middle() void PairLJClass2::compute_outer(int eflag, int vflag) { - int i,j,ii,jj,inum,jnum,itype,jtype; - double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair; - double rsq,rinv,r2inv,r3inv,r6inv,forcelj,factor_lj,rsw; - int *ilist,*jlist,*numneigh,**firstneigh; + int i, j, ii, jj, inum, jnum, itype, jtype; + double xtmp, ytmp, ztmp, delx, dely, delz, evdwl, fpair; + double rsq, rinv, r2inv, r3inv, r6inv, forcelj, factor_lj, rsw; + int *ilist, *jlist, *numneigh, **firstneigh; evdwl = 0.0; - ev_init(eflag,vflag); + ev_init(eflag, vflag); double **x = atom->x; double **f = atom->f; @@ -327,8 +322,8 @@ void PairLJClass2::compute_outer(int eflag, int vflag) double cut_in_on = cut_respa[3]; double cut_in_diff = cut_in_on - cut_in_off; - double cut_in_off_sq = cut_in_off*cut_in_off; - double cut_in_on_sq = cut_in_on*cut_in_on; + double cut_in_off_sq = cut_in_off * cut_in_off; + double cut_in_on_sq = cut_in_on * cut_in_on; // loop over neighbors of my atoms @@ -349,55 +344,54 @@ void PairLJClass2::compute_outer(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]; if (rsq < cutsq[itype][jtype]) { if (rsq > cut_in_off_sq) { - r2inv = 1.0/rsq; - rinv = sqrt(r2inv); - r3inv = r2inv*rinv; - r6inv = r3inv*r3inv; - forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]); - fpair = factor_lj*forcelj*r2inv; + r2inv = 1.0 / rsq; + rinv = sqrt(r2inv); + r3inv = r2inv * rinv; + r6inv = r3inv * r3inv; + forcelj = r6inv * (lj1[itype][jtype] * r3inv - lj2[itype][jtype]); + fpair = factor_lj * forcelj * r2inv; if (rsq < cut_in_on_sq) { - rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff; - fpair *= rsw*rsw*(3.0 - 2.0*rsw); + rsw = (sqrt(rsq) - cut_in_off) / cut_in_diff; + fpair *= rsw * rsw * (3.0 - 2.0 * rsw); } - 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) { - r2inv = 1.0/rsq; - rinv = sqrt(r2inv); - r3inv = r2inv*rinv; - r6inv = r3inv*r3inv; - evdwl = r6inv*(lj3[itype][jtype]*r3inv-lj4[itype][jtype]) - - offset[itype][jtype]; + r2inv = 1.0 / rsq; + rinv = sqrt(r2inv); + r3inv = r2inv * rinv; + r6inv = r3inv * r3inv; + evdwl = r6inv * (lj3[itype][jtype] * r3inv - lj4[itype][jtype]) - offset[itype][jtype]; evdwl *= factor_lj; } if (vflag) { if (rsq <= cut_in_off_sq) { - r2inv = 1.0/rsq; - rinv = sqrt(r2inv); - r3inv = r2inv*rinv; - r6inv = r3inv*r3inv; - forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]); - fpair = factor_lj*forcelj*r2inv; + r2inv = 1.0 / rsq; + rinv = sqrt(r2inv); + r3inv = r2inv * rinv; + r6inv = r3inv * r3inv; + forcelj = r6inv * (lj1[itype][jtype] * r3inv - lj2[itype][jtype]); + fpair = factor_lj * forcelj * r2inv; } else if (rsq < cut_in_on_sq) - fpair = factor_lj*forcelj*r2inv; + fpair = factor_lj * forcelj * r2inv; } - if (evflag) ev_tally(i,j,nlocal,newton_pair,evdwl,0.0,fpair,delx,dely,delz); + if (evflag) ev_tally(i, j, nlocal, newton_pair, evdwl, 0.0, fpair, delx, dely, delz); } } } @@ -411,21 +405,20 @@ void PairLJClass2::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(cutsq, n + 1, n + 1, "pair:cutsq"); - memory->create(cut,n+1,n+1,"pair:cut"); - memory->create(epsilon,n+1,n+1,"pair:epsilon"); - memory->create(sigma,n+1,n+1,"pair:sigma"); - memory->create(lj1,n+1,n+1,"pair:lj1"); - memory->create(lj2,n+1,n+1,"pair:lj2"); - memory->create(lj3,n+1,n+1,"pair:lj3"); - memory->create(lj4,n+1,n+1,"pair:lj4"); - memory->create(offset,n+1,n+1,"pair:offset"); + memory->create(cut, n + 1, n + 1, "pair:cut"); + memory->create(epsilon, n + 1, n + 1, "pair:epsilon"); + memory->create(sigma, n + 1, n + 1, "pair:sigma"); + memory->create(lj1, n + 1, n + 1, "pair:lj1"); + memory->create(lj2, n + 1, n + 1, "pair:lj2"); + memory->create(lj3, n + 1, n + 1, "pair:lj3"); + memory->create(lj4, n + 1, n + 1, "pair:lj4"); + memory->create(offset, n + 1, n + 1, "pair:offset"); } /* ---------------------------------------------------------------------- @@ -434,14 +427,14 @@ void PairLJClass2::allocate() void PairLJClass2::settings(int narg, char **arg) { - if (narg != 1) error->all(FLERR,"Illegal pair_style command"); + if (narg != 1) error->all(FLERR, "Illegal pair_style command"); - cut_global = utils::numeric(FLERR,arg[0],false,lmp); + cut_global = utils::numeric(FLERR, arg[0], 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; j <= atom->ntypes; j++) if (setflag[i][j]) cut[i][j] = cut_global; @@ -454,22 +447,22 @@ void PairLJClass2::settings(int narg, char **arg) void PairLJClass2::coeff(int narg, char **arg) { - if (narg < 4 || narg > 5) error->all(FLERR,"Incorrect args for pair coefficients"); + if (narg < 4 || narg > 5) 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 epsilon_one = utils::numeric(FLERR,arg[2],false,lmp); - double sigma_one = utils::numeric(FLERR,arg[3],false,lmp); + double epsilon_one = utils::numeric(FLERR, arg[2], false, lmp); + double sigma_one = utils::numeric(FLERR, arg[3], false, lmp); double cut_one = cut_global; - if (narg == 5) cut_one = utils::numeric(FLERR,arg[4],false,lmp); + if (narg == 5) cut_one = utils::numeric(FLERR, arg[4], 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++) { epsilon[i][j] = epsilon_one; sigma[i][j] = sigma_one; cut[i][j] = cut_one; @@ -478,7 +471,7 @@ void PairLJClass2::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"); } /* ---------------------------------------------------------------------- @@ -489,28 +482,22 @@ void PairLJClass2::init_style() { // request regular or rRESPA neighbor list - int irequest; - int respa = 0; + int list_style = NeighConst::REQ_DEFAULT; - if (update->whichflag == 1 && utils::strmatch(update->integrate_style,"^respa")) { - if (((Respa *) update->integrate)->level_inner >= 0) respa = 1; - if (((Respa *) update->integrate)->level_middle >= 0) respa = 2; + if (update->whichflag == 1 && utils::strmatch(update->integrate_style, "^respa")) { + auto respa = (Respa *) update->integrate; + if (respa->level_inner >= 0) list_style = NeighConst::REQ_RESPA_INOUT; + if (respa->level_middle >= 0) list_style = NeighConst::REQ_RESPA_ALL; } - - irequest = neighbor->request(this,instance_me); - - if (respa >= 1) { - neighbor->requests[irequest]->respaouter = 1; - neighbor->requests[irequest]->respainner = 1; - } - if (respa == 2) neighbor->requests[irequest]->respamiddle = 1; + neighbor->add_request(this, list_style); // set rRESPA cutoffs - if (utils::strmatch(update->integrate_style,"^respa") && + if (utils::strmatch(update->integrate_style, "^respa") && ((Respa *) update->integrate)->level_inner >= 0) cut_respa = ((Respa *) update->integrate)->cutoff; - else cut_respa = nullptr; + else + cut_respa = nullptr; } /* ---------------------------------------------------------------------- @@ -523,23 +510,22 @@ double PairLJClass2::init_one(int i, int j) // mix distance via user-defined rule if (setflag[i][j] == 0) { - epsilon[i][j] = 2.0 * sqrt(epsilon[i][i]*epsilon[j][j]) * - pow(sigma[i][i],3.0) * pow(sigma[j][j],3.0) / - (pow(sigma[i][i],6.0) + pow(sigma[j][j],6.0)); - sigma[i][j] = - pow((0.5 * (pow(sigma[i][i],6.0) + pow(sigma[j][j],6.0))),1.0/6.0); - cut[i][j] = mix_distance(cut[i][i],cut[j][j]); + epsilon[i][j] = 2.0 * sqrt(epsilon[i][i] * epsilon[j][j]) * pow(sigma[i][i], 3.0) * + pow(sigma[j][j], 3.0) / (pow(sigma[i][i], 6.0) + pow(sigma[j][j], 6.0)); + sigma[i][j] = pow((0.5 * (pow(sigma[i][i], 6.0) + pow(sigma[j][j], 6.0))), 1.0 / 6.0); + cut[i][j] = mix_distance(cut[i][i], cut[j][j]); } - lj1[i][j] = 18.0 * epsilon[i][j] * pow(sigma[i][j],9.0); - lj2[i][j] = 18.0 * epsilon[i][j] * pow(sigma[i][j],6.0); - lj3[i][j] = 2.0 * epsilon[i][j] * pow(sigma[i][j],9.0); - lj4[i][j] = 3.0 * epsilon[i][j] * pow(sigma[i][j],6.0); + lj1[i][j] = 18.0 * epsilon[i][j] * pow(sigma[i][j], 9.0); + lj2[i][j] = 18.0 * epsilon[i][j] * pow(sigma[i][j], 6.0); + lj3[i][j] = 2.0 * epsilon[i][j] * pow(sigma[i][j], 9.0); + lj4[i][j] = 3.0 * epsilon[i][j] * pow(sigma[i][j], 6.0); if (offset_flag && (cut[i][j] > 0.0)) { double ratio = sigma[i][j] / cut[i][j]; - offset[i][j] = epsilon[i][j] * (2.0*pow(ratio,9.0) - 3.0*pow(ratio,6.0)); - } else offset[i][j] = 0.0; + offset[i][j] = epsilon[i][j] * (2.0 * pow(ratio, 9.0) - 3.0 * pow(ratio, 6.0)); + } else + offset[i][j] = 0.0; lj1[j][i] = lj1[i][j]; lj2[j][i] = lj2[i][j]; @@ -550,7 +536,7 @@ double PairLJClass2::init_one(int i, int j) // check interior rRESPA cutoff if (cut_respa && cut[i][j] < cut_respa[3]) - error->all(FLERR,"Pair cutoff < Respa interior cutoff"); + error->all(FLERR, "Pair cutoff < Respa interior cutoff"); // compute I,J contribution to long-range tail correction // count total # of atoms of type I and J via Allreduce @@ -559,22 +545,21 @@ double PairLJClass2::init_one(int i, int j) int *type = atom->type; int nlocal = atom->nlocal; - double count[2],all[2]; + double count[2], all[2]; count[0] = count[1] = 0.0; for (int k = 0; k < nlocal; k++) { if (type[k] == i) count[0] += 1.0; if (type[k] == j) count[1] += 1.0; } - MPI_Allreduce(count,all,2,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(count, all, 2, MPI_DOUBLE, MPI_SUM, world); - double sig3 = sigma[i][j]*sigma[i][j]*sigma[i][j]; - double sig6 = sig3*sig3; - double rc3 = cut[i][j]*cut[i][j]*cut[i][j]; - double rc6 = rc3*rc3; - etail_ij = 2.0*MY_PI*all[0]*all[1]*epsilon[i][j] * - sig6 * (sig3 - 3.0*rc3) / (3.0*rc6); - ptail_ij = 2.0*MY_PI*all[0]*all[1]*epsilon[i][j] * - sig6 * (sig3 - 2.0*rc3) / rc6; + double sig3 = sigma[i][j] * sigma[i][j] * sigma[i][j]; + double sig6 = sig3 * sig3; + double rc3 = cut[i][j] * cut[i][j] * cut[i][j]; + double rc6 = rc3 * rc3; + etail_ij = + 2.0 * MY_PI * all[0] * all[1] * epsilon[i][j] * sig6 * (sig3 - 3.0 * rc3) / (3.0 * rc6); + ptail_ij = 2.0 * MY_PI * all[0] * all[1] * epsilon[i][j] * sig6 * (sig3 - 2.0 * rc3) / rc6; } return cut[i][j]; @@ -588,14 +573,14 @@ void PairLJClass2::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(&epsilon[i][j],sizeof(double),1,fp); - fwrite(&sigma[i][j],sizeof(double),1,fp); - fwrite(&cut[i][j],sizeof(double),1,fp); + fwrite(&epsilon[i][j], sizeof(double), 1, fp); + fwrite(&sigma[i][j], sizeof(double), 1, fp); + fwrite(&cut[i][j], sizeof(double), 1, fp); } } } @@ -609,21 +594,21 @@ void PairLJClass2::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,&epsilon[i][j],sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&sigma[i][j],sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&cut[i][j],sizeof(double),1,fp,nullptr,error); + utils::sfread(FLERR, &epsilon[i][j], sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &sigma[i][j], sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &cut[i][j], sizeof(double), 1, fp, nullptr, error); } - MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&epsilon[i][j], 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&sigma[i][j], 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut[i][j], 1, MPI_DOUBLE, 0, world); } } } @@ -634,10 +619,10 @@ void PairLJClass2::read_restart(FILE *fp) void PairLJClass2::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(&tail_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); + fwrite(&tail_flag, sizeof(int), 1, fp); } /* ---------------------------------------------------------------------- @@ -648,26 +633,24 @@ void PairLJClass2::read_restart_settings(FILE *fp) { int me = comm->me; if (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,&tail_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); + utils::sfread(FLERR, &tail_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(&tail_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); + MPI_Bcast(&tail_flag, 1, MPI_INT, 0, world); } - /* ---------------------------------------------------------------------- proc 0 writes to data file ------------------------------------------------------------------------- */ void PairLJClass2::write_data(FILE *fp) { - for (int i = 1; i <= atom->ntypes; i++) - fprintf(fp,"%d %g %g\n",i,epsilon[i][i],sigma[i][i]); + for (int i = 1; i <= atom->ntypes; i++) fprintf(fp, "%d %g %g\n", i, epsilon[i][i], sigma[i][i]); } /* ---------------------------------------------------------------------- @@ -678,27 +661,25 @@ void PairLJClass2::write_data_all(FILE *fp) { for (int i = 1; i <= atom->ntypes; i++) for (int j = i; j <= atom->ntypes; j++) - fprintf(fp,"%d %d %g %g %g\n",i,j,epsilon[i][j],sigma[i][j],cut[i][j]); + fprintf(fp, "%d %d %g %g %g\n", i, j, epsilon[i][j], sigma[i][j], cut[i][j]); } /* ---------------------------------------------------------------------- */ double PairLJClass2::single(int /*i*/, int /*j*/, int itype, int jtype, double rsq, - double /*factor_coul*/, double factor_lj, - double &fforce) + double /*factor_coul*/, double factor_lj, double &fforce) { - double r2inv,rinv,r3inv,r6inv,forcelj,philj; + double r2inv, rinv, r3inv, r6inv, forcelj, philj; - r2inv = 1.0/rsq; + r2inv = 1.0 / rsq; rinv = sqrt(r2inv); - r3inv = r2inv*rinv; - r6inv = r3inv*r3inv; - forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]); - fforce = factor_lj*forcelj*r2inv; + r3inv = r2inv * rinv; + r6inv = r3inv * r3inv; + forcelj = r6inv * (lj1[itype][jtype] * r3inv - lj2[itype][jtype]); + fforce = factor_lj * forcelj * r2inv; - philj = r6inv*(lj3[itype][jtype]*r3inv-lj4[itype][jtype]) - - offset[itype][jtype]; - return factor_lj*philj; + philj = r6inv * (lj3[itype][jtype] * r3inv - lj4[itype][jtype]) - offset[itype][jtype]; + return factor_lj * philj; } /* ---------------------------------------------------------------------- */ @@ -706,7 +687,7 @@ double PairLJClass2::single(int /*i*/, int /*j*/, int itype, int jtype, double r void *PairLJClass2::extract(const char *str, int &dim) { dim = 2; - if (strcmp(str,"epsilon") == 0) return (void *) epsilon; - if (strcmp(str,"sigma") == 0) return (void *) sigma; + if (strcmp(str, "epsilon") == 0) return (void *) epsilon; + if (strcmp(str, "sigma") == 0) return (void *) sigma; return nullptr; } diff --git a/src/CLASS2/pair_lj_class2_coul_cut.cpp b/src/CLASS2/pair_lj_class2_coul_cut.cpp index c4ac5ae7c8..179c91b48e 100644 --- a/src/CLASS2/pair_lj_class2_coul_cut.cpp +++ b/src/CLASS2/pair_lj_class2_coul_cut.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -14,8 +13,6 @@ #include "pair_lj_class2_coul_cut.h" -#include -#include #include "atom.h" #include "comm.h" #include "force.h" @@ -25,6 +22,8 @@ #include "memory.h" #include "error.h" +#include +#include using namespace LAMMPS_NS; using namespace MathConst; @@ -259,7 +258,7 @@ void PairLJClass2CoulCut::init_style() if (!atom->q_flag) error->all(FLERR,"Pair style lj/class2/coul/cut requires atom attribute q"); - neighbor->request(this,instance_me); + neighbor->add_request(this); } /* ---------------------------------------------------------------------- diff --git a/src/CLASS2/pair_lj_class2_coul_long.cpp b/src/CLASS2/pair_lj_class2_coul_long.cpp index 63265509c9..fb4e96d337 100644 --- a/src/CLASS2/pair_lj_class2_coul_long.cpp +++ b/src/CLASS2/pair_lj_class2_coul_long.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -14,32 +13,31 @@ #include "pair_lj_class2_coul_long.h" -#include -#include #include "atom.h" #include "comm.h" +#include "error.h" #include "force.h" #include "kspace.h" -#include "update.h" -#include "respa.h" -#include "neighbor.h" -#include "neigh_list.h" -#include "neigh_request.h" #include "math_const.h" #include "memory.h" -#include "error.h" +#include "neigh_list.h" +#include "neighbor.h" +#include "respa.h" +#include "update.h" +#include +#include using namespace LAMMPS_NS; using namespace MathConst; -#define EWALD_F 1.12837917 -#define EWALD_P 0.3275911 -#define A1 0.254829592 -#define A2 -0.284496736 -#define A3 1.421413741 -#define A4 -1.453152027 -#define A5 1.061405429 +#define EWALD_F 1.12837917 +#define EWALD_P 0.3275911 +#define A1 0.254829592 +#define A2 -0.284496736 +#define A3 1.421413741 +#define A4 -1.453152027 +#define A5 1.061405429 /* ---------------------------------------------------------------------- */ @@ -79,16 +77,16 @@ PairLJClass2CoulLong::~PairLJClass2CoulLong() void PairLJClass2CoulLong::compute(int eflag, int vflag) { - int i,j,ii,jj,inum,jnum,itable,itype,jtype; - double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair; - double fraction,table; - double rsq,r,rinv,r2inv,r3inv,r6inv,forcecoul,forcelj; - double grij,expm2,prefactor,t,erfc; - double factor_coul,factor_lj; - int *ilist,*jlist,*numneigh,**firstneigh; + int i, j, ii, jj, inum, jnum, itable, itype, jtype; + double qtmp, xtmp, ytmp, ztmp, delx, dely, delz, evdwl, ecoul, fpair; + double fraction, table; + double rsq, r, rinv, r2inv, r3inv, r6inv, forcecoul, forcelj; + double grij, expm2, prefactor, t, erfc; + double factor_coul, factor_lj; + int *ilist, *jlist, *numneigh, **firstneigh; evdwl = ecoul = 0.0; - ev_init(eflag,vflag); + ev_init(eflag, vflag); double **x = atom->x; double **f = atom->f; @@ -126,75 +124,77 @@ void PairLJClass2CoulLong::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]; if (rsq < cutsq[itype][jtype]) { - r2inv = 1.0/rsq; + r2inv = 1.0 / rsq; if (rsq < cut_coulsq) { if (!ncoultablebits || rsq <= tabinnersq) { r = sqrt(rsq); grij = g_ewald * r; - expm2 = exp(-grij*grij); - t = 1.0 / (1.0 + EWALD_P*grij); - erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2; - prefactor = qqrd2e * qtmp*q[j]/r; - forcecoul = prefactor * (erfc + EWALD_F*grij*expm2); - if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; + expm2 = exp(-grij * grij); + t = 1.0 / (1.0 + EWALD_P * grij); + erfc = t * (A1 + t * (A2 + t * (A3 + t * (A4 + t * A5)))) * expm2; + prefactor = qqrd2e * qtmp * q[j] / r; + forcecoul = prefactor * (erfc + EWALD_F * grij * expm2); + if (factor_coul < 1.0) forcecoul -= (1.0 - factor_coul) * prefactor; } else { union_int_float_t rsq_lookup; rsq_lookup.f = rsq; itable = rsq_lookup.i & ncoulmask; itable >>= ncoulshiftbits; fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable]; - table = ftable[itable] + fraction*dftable[itable]; - forcecoul = qtmp*q[j] * table; + table = ftable[itable] + fraction * dftable[itable]; + forcecoul = qtmp * q[j] * table; if (factor_coul < 1.0) { - table = ctable[itable] + fraction*dctable[itable]; - prefactor = qtmp*q[j] * table; - forcecoul -= (1.0-factor_coul)*prefactor; + table = ctable[itable] + fraction * dctable[itable]; + prefactor = qtmp * q[j] * table; + forcecoul -= (1.0 - factor_coul) * prefactor; } } - } else forcecoul = 0.0; + } else + forcecoul = 0.0; if (rsq < cut_ljsq[itype][jtype]) { rinv = sqrt(r2inv); - r3inv = r2inv*rinv; - r6inv = r3inv*r3inv; - forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]); - } else forcelj = 0.0; + r3inv = r2inv * rinv; + r6inv = r3inv * r3inv; + forcelj = r6inv * (lj1[itype][jtype] * r3inv - lj2[itype][jtype]); + } else + forcelj = 0.0; - fpair = (forcecoul + factor_lj*forcelj) * r2inv; + fpair = (forcecoul + factor_lj * forcelj) * r2inv; - 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 (rsq < cut_coulsq) { if (!ncoultablebits || rsq <= tabinnersq) - ecoul = prefactor*erfc; + ecoul = prefactor * erfc; else { - table = etable[itable] + fraction*detable[itable]; - ecoul = qtmp*q[j] * table; + table = etable[itable] + fraction * detable[itable]; + ecoul = qtmp * q[j] * table; } - if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor; - } else ecoul = 0.0; + if (factor_coul < 1.0) ecoul -= (1.0 - factor_coul) * prefactor; + } else + ecoul = 0.0; if (rsq < cut_ljsq[itype][jtype]) { - evdwl = r6inv*(lj3[itype][jtype]*r3inv-lj4[itype][jtype]) - - offset[itype][jtype]; + evdwl = r6inv * (lj3[itype][jtype] * r3inv - lj4[itype][jtype]) - offset[itype][jtype]; evdwl *= factor_lj; - } else evdwl = 0.0; + } else + evdwl = 0.0; } - if (evflag) ev_tally(i,j,nlocal,newton_pair, - evdwl,ecoul,fpair,delx,dely,delz); + if (evflag) ev_tally(i, j, nlocal, newton_pair, evdwl, ecoul, fpair, delx, dely, delz); } } } @@ -206,11 +206,11 @@ void PairLJClass2CoulLong::compute(int eflag, int vflag) void PairLJClass2CoulLong::compute_inner() { - int i,j,ii,jj,inum,jnum,itype,jtype; - double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,fpair; - double rsq,rinv,r2inv,r3inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj; + int i, j, ii, jj, inum, jnum, itype, jtype; + double qtmp, xtmp, ytmp, ztmp, delx, dely, delz, fpair; + double rsq, rinv, r2inv, r3inv, r6inv, forcecoul, forcelj, factor_coul, factor_lj; double rsw; - int *ilist,*jlist,*numneigh,**firstneigh; + int *ilist, *jlist, *numneigh, **firstneigh; double **x = atom->x; double **f = atom->f; @@ -231,8 +231,8 @@ void PairLJClass2CoulLong::compute_inner() double cut_out_off = cut_respa[1]; double cut_out_diff = cut_out_off - cut_out_on; - double cut_out_on_sq = cut_out_on*cut_out_on; - double cut_out_off_sq = cut_out_off*cut_out_off; + double cut_out_on_sq = cut_out_on * cut_out_on; + double cut_out_off_sq = cut_out_off * cut_out_off; // loop over neighbors of my atoms @@ -255,34 +255,35 @@ void PairLJClass2CoulLong::compute_inner() 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; if (rsq < cut_out_off_sq) { - r2inv = 1.0/rsq; - forcecoul = qqrd2e * qtmp*q[j]*sqrt(r2inv); - if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*forcecoul; + r2inv = 1.0 / rsq; + forcecoul = qqrd2e * qtmp * q[j] * sqrt(r2inv); + if (factor_coul < 1.0) forcecoul -= (1.0 - factor_coul) * forcecoul; jtype = type[j]; if (rsq < cut_ljsq[itype][jtype]) { - rinv = sqrt(r2inv); - r3inv = r2inv*rinv; - r6inv = r3inv*r3inv; - forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]); - } else forcelj = 0.0; + rinv = sqrt(r2inv); + r3inv = r2inv * rinv; + r6inv = r3inv * r3inv; + forcelj = r6inv * (lj1[itype][jtype] * r3inv - lj2[itype][jtype]); + } else + forcelj = 0.0; - fpair = (forcecoul + factor_lj*forcelj) * r2inv; + fpair = (forcecoul + factor_lj * forcelj) * r2inv; if (rsq > cut_out_on_sq) { - rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff; - fpair *= 1.0 + rsw*rsw*(2.0*rsw-3.0); + rsw = (sqrt(rsq) - cut_out_on) / cut_out_diff; + fpair *= 1.0 + rsw * rsw * (2.0 * rsw - 3.0); } - 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; } } } @@ -293,11 +294,11 @@ void PairLJClass2CoulLong::compute_inner() void PairLJClass2CoulLong::compute_middle() { - int i,j,ii,jj,inum,jnum,itype,jtype; - double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,fpair; - double rsq,rinv,r2inv,r3inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj; + int i, j, ii, jj, inum, jnum, itype, jtype; + double qtmp, xtmp, ytmp, ztmp, delx, dely, delz, fpair; + double rsq, rinv, r2inv, r3inv, r6inv, forcecoul, forcelj, factor_coul, factor_lj; double rsw; - int *ilist,*jlist,*numneigh,**firstneigh; + int *ilist, *jlist, *numneigh, **firstneigh; double **x = atom->x; double **f = atom->f; @@ -321,10 +322,10 @@ void PairLJClass2CoulLong::compute_middle() double cut_in_diff = cut_in_on - cut_in_off; double cut_out_diff = cut_out_off - cut_out_on; - double cut_in_off_sq = cut_in_off*cut_in_off; - double cut_in_on_sq = cut_in_on*cut_in_on; - double cut_out_on_sq = cut_out_on*cut_out_on; - double cut_out_off_sq = cut_out_off*cut_out_off; + double cut_in_off_sq = cut_in_off * cut_in_off; + double cut_in_on_sq = cut_in_on * cut_in_on; + double cut_out_on_sq = cut_out_on * cut_out_on; + double cut_out_off_sq = cut_out_off * cut_out_off; // loop over neighbors of my atoms @@ -347,38 +348,39 @@ void PairLJClass2CoulLong::compute_middle() 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; if (rsq < cut_out_off_sq && rsq > cut_in_off_sq) { - r2inv = 1.0/rsq; - forcecoul = qqrd2e * qtmp*q[j]*sqrt(r2inv); - if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*forcecoul; + r2inv = 1.0 / rsq; + forcecoul = qqrd2e * qtmp * q[j] * sqrt(r2inv); + if (factor_coul < 1.0) forcecoul -= (1.0 - factor_coul) * forcecoul; jtype = type[j]; if (rsq < cut_ljsq[itype][jtype]) { - rinv = sqrt(r2inv); - r3inv = r2inv*rinv; - r6inv = r3inv*r3inv; - forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]); - } else forcelj = 0.0; + rinv = sqrt(r2inv); + r3inv = r2inv * rinv; + r6inv = r3inv * r3inv; + forcelj = r6inv * (lj1[itype][jtype] * r3inv - lj2[itype][jtype]); + } else + forcelj = 0.0; - fpair = (forcecoul + factor_lj*forcelj) * r2inv; + fpair = (forcecoul + factor_lj * forcelj) * r2inv; if (rsq < cut_in_on_sq) { - rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff; - fpair *= rsw*rsw*(3.0 - 2.0*rsw); + rsw = (sqrt(rsq) - cut_in_off) / cut_in_diff; + fpair *= rsw * rsw * (3.0 - 2.0 * rsw); } if (rsq > cut_out_on_sq) { - rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff; - fpair *= 1.0 + rsw*rsw*(2.0*rsw - 3.0); + rsw = (sqrt(rsq) - cut_out_on) / cut_out_diff; + fpair *= 1.0 + rsw * rsw * (2.0 * rsw - 3.0); } - 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; } } } @@ -389,17 +391,17 @@ void PairLJClass2CoulLong::compute_middle() void PairLJClass2CoulLong::compute_outer(int eflag, int vflag) { - int i,j,ii,jj,inum,jnum,itype,jtype,itable; - double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair; - double fraction,table; - double r,rinv,r2inv,r3inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj; - double grij,expm2,prefactor,t,erfc; + int i, j, ii, jj, inum, jnum, itype, jtype, itable; + double qtmp, xtmp, ytmp, ztmp, delx, dely, delz, evdwl, ecoul, fpair; + double fraction, table; + double r, rinv, r2inv, r3inv, r6inv, forcecoul, forcelj, factor_coul, factor_lj; + double grij, expm2, prefactor, t, erfc; double rsw; - int *ilist,*jlist,*numneigh,**firstneigh; + int *ilist, *jlist, *numneigh, **firstneigh; double rsq; evdwl = ecoul = 0.0; - ev_init(eflag,vflag); + ev_init(eflag, vflag); double **x = atom->x; double **f = atom->f; @@ -420,8 +422,8 @@ void PairLJClass2CoulLong::compute_outer(int eflag, int vflag) double cut_in_on = cut_respa[3]; double cut_in_diff = cut_in_on - cut_in_off; - double cut_in_off_sq = cut_in_off*cut_in_off; - double cut_in_on_sq = cut_in_on*cut_in_on; + double cut_in_off_sq = cut_in_off * cut_in_off; + double cut_in_on_sq = cut_in_on * cut_in_on; // loop over neighbors of my atoms @@ -444,32 +446,30 @@ void PairLJClass2CoulLong::compute_outer(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]; if (rsq < cutsq[itype][jtype]) { - r2inv = 1.0/rsq; + r2inv = 1.0 / rsq; if (rsq < cut_coulsq) { if (!ncoultablebits || rsq <= tabinnersq) { r = sqrt(rsq); grij = g_ewald * r; - expm2 = exp(-grij*grij); - t = 1.0 / (1.0 + EWALD_P*grij); - erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2; - prefactor = qqrd2e * qtmp*q[j]/r; - forcecoul = prefactor * (erfc + EWALD_F*grij*expm2 - 1.0); + expm2 = exp(-grij * grij); + t = 1.0 / (1.0 + EWALD_P * grij); + erfc = t * (A1 + t * (A2 + t * (A3 + t * (A4 + t * A5)))) * expm2; + prefactor = qqrd2e * qtmp * q[j] / r; + forcecoul = prefactor * (erfc + EWALD_F * grij * expm2 - 1.0); if (rsq > cut_in_off_sq) { if (rsq < cut_in_on_sq) { - rsw = (r - cut_in_off)/cut_in_diff; - forcecoul += prefactor*rsw*rsw*(3.0 - 2.0*rsw); + rsw = (r - cut_in_off) / cut_in_diff; + forcecoul += prefactor * rsw * rsw * (3.0 - 2.0 * rsw); if (factor_coul < 1.0) - forcecoul -= - (1.0-factor_coul)*prefactor*rsw*rsw*(3.0 - 2.0*rsw); + forcecoul -= (1.0 - factor_coul) * prefactor * rsw * rsw * (3.0 - 2.0 * rsw); } else { forcecoul += prefactor; - if (factor_coul < 1.0) - forcecoul -= (1.0-factor_coul)*prefactor; + if (factor_coul < 1.0) forcecoul -= (1.0 - factor_coul) * prefactor; } } } else { @@ -478,96 +478,99 @@ void PairLJClass2CoulLong::compute_outer(int eflag, int vflag) itable = rsq_lookup.i & ncoulmask; itable >>= ncoulshiftbits; fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable]; - table = ftable[itable] + fraction*dftable[itable]; - forcecoul = qtmp*q[j] * table; + table = ftable[itable] + fraction * dftable[itable]; + forcecoul = qtmp * q[j] * table; if (factor_coul < 1.0) { - table = ctable[itable] + fraction*dctable[itable]; - prefactor = qtmp*q[j] * table; - forcecoul -= (1.0-factor_coul)*prefactor; + table = ctable[itable] + fraction * dctable[itable]; + prefactor = qtmp * q[j] * table; + forcecoul -= (1.0 - factor_coul) * prefactor; } } - } else forcecoul = 0.0; + } else + forcecoul = 0.0; if (rsq < cut_ljsq[itype][jtype] && rsq > cut_in_off_sq) { - rinv = sqrt(r2inv); - r3inv = r2inv*rinv; - r6inv = r3inv*r3inv; - forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]); + rinv = sqrt(r2inv); + r3inv = r2inv * rinv; + r6inv = r3inv * r3inv; + forcelj = r6inv * (lj1[itype][jtype] * r3inv - lj2[itype][jtype]); if (rsq < cut_in_on_sq) { - rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff; - forcelj *= rsw*rsw*(3.0 - 2.0*rsw); + rsw = (sqrt(rsq) - cut_in_off) / cut_in_diff; + forcelj *= rsw * rsw * (3.0 - 2.0 * rsw); } - } else forcelj = 0.0; + } else + forcelj = 0.0; fpair = (forcecoul + forcelj) * r2inv; - 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 (rsq < cut_coulsq) { if (!ncoultablebits || rsq <= tabinnersq) { - ecoul = prefactor*erfc; - if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor; + ecoul = prefactor * erfc; + if (factor_coul < 1.0) ecoul -= (1.0 - factor_coul) * prefactor; } else { - table = etable[itable] + fraction*detable[itable]; - ecoul = qtmp*q[j] * table; + table = etable[itable] + fraction * detable[itable]; + ecoul = qtmp * q[j] * table; if (factor_coul < 1.0) { - table = ptable[itable] + fraction*dptable[itable]; - prefactor = qtmp*q[j] * table; - ecoul -= (1.0-factor_coul)*prefactor; + table = ptable[itable] + fraction * dptable[itable]; + prefactor = qtmp * q[j] * table; + ecoul -= (1.0 - factor_coul) * prefactor; } } - } else ecoul = 0.0; + } else + ecoul = 0.0; if (rsq < cut_ljsq[itype][jtype]) { - rinv = sqrt(r2inv); - r3inv = r2inv*rinv; - r6inv = r3inv*r3inv; - evdwl = r6inv*(lj3[itype][jtype]*r3inv-lj4[itype][jtype]) - - offset[itype][jtype]; + rinv = sqrt(r2inv); + r3inv = r2inv * rinv; + r6inv = r3inv * r3inv; + evdwl = r6inv * (lj3[itype][jtype] * r3inv - lj4[itype][jtype]) - offset[itype][jtype]; evdwl *= factor_lj; - } else evdwl = 0.0; + } else + evdwl = 0.0; } if (vflag) { if (rsq < cut_coulsq) { if (!ncoultablebits || rsq <= tabinnersq) { - forcecoul = prefactor * (erfc + EWALD_F*grij*expm2); - if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; + forcecoul = prefactor * (erfc + EWALD_F * grij * expm2); + if (factor_coul < 1.0) forcecoul -= (1.0 - factor_coul) * prefactor; } else { - table = vtable[itable] + fraction*dvtable[itable]; - forcecoul = qtmp*q[j] * table; + table = vtable[itable] + fraction * dvtable[itable]; + forcecoul = qtmp * q[j] * table; if (factor_coul < 1.0) { - table = ptable[itable] + fraction*dptable[itable]; - prefactor = qtmp*q[j] * table; - forcecoul -= (1.0-factor_coul)*prefactor; + table = ptable[itable] + fraction * dptable[itable]; + prefactor = qtmp * q[j] * table; + forcecoul -= (1.0 - factor_coul) * prefactor; } } - } else forcecoul = 0.0; + } else + forcecoul = 0.0; if (rsq <= cut_in_off_sq) { rinv = sqrt(r2inv); - r3inv = r2inv*rinv; - r6inv = r3inv*r3inv; - forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]); + r3inv = r2inv * rinv; + r6inv = r3inv * r3inv; + forcelj = r6inv * (lj1[itype][jtype] * r3inv - lj2[itype][jtype]); } else if (rsq <= cut_in_on_sq) { - rinv = sqrt(r2inv); - r3inv = r2inv*rinv; - r6inv = r3inv*r3inv; - forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]); + rinv = sqrt(r2inv); + r3inv = r2inv * rinv; + r6inv = r3inv * r3inv; + forcelj = r6inv * (lj1[itype][jtype] * r3inv - lj2[itype][jtype]); } - fpair = (forcecoul + factor_lj*forcelj) * r2inv; + fpair = (forcecoul + factor_lj * forcelj) * r2inv; } - if (evflag) ev_tally(i,j,nlocal,newton_pair, - evdwl,ecoul,fpair,delx,dely,delz); + if (evflag) ev_tally(i, j, nlocal, newton_pair, evdwl, ecoul, fpair, delx, dely, delz); } } } @@ -582,22 +585,21 @@ void PairLJClass2CoulLong::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(cutsq, n + 1, n + 1, "pair:cutsq"); - memory->create(cut_lj,n+1,n+1,"pair:cut_lj"); - memory->create(cut_ljsq,n+1,n+1,"pair:cut_ljsq"); - memory->create(epsilon,n+1,n+1,"pair:epsilon"); - memory->create(sigma,n+1,n+1,"pair:sigma"); - memory->create(lj1,n+1,n+1,"pair:lj1"); - memory->create(lj2,n+1,n+1,"pair:lj2"); - memory->create(lj3,n+1,n+1,"pair:lj3"); - memory->create(lj4,n+1,n+1,"pair:lj4"); - memory->create(offset,n+1,n+1,"pair:offset"); + memory->create(cut_lj, n + 1, n + 1, "pair:cut_lj"); + memory->create(cut_ljsq, n + 1, n + 1, "pair:cut_ljsq"); + memory->create(epsilon, n + 1, n + 1, "pair:epsilon"); + memory->create(sigma, n + 1, n + 1, "pair:sigma"); + memory->create(lj1, n + 1, n + 1, "pair:lj1"); + memory->create(lj2, n + 1, n + 1, "pair:lj2"); + memory->create(lj3, n + 1, n + 1, "pair:lj3"); + memory->create(lj4, n + 1, n + 1, "pair:lj4"); + memory->create(offset, n + 1, n + 1, "pair:offset"); } /* ---------------------------------------------------------------------- @@ -606,16 +608,18 @@ void PairLJClass2CoulLong::allocate() void PairLJClass2CoulLong::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_lj_global = utils::numeric(FLERR,arg[0],false,lmp); - if (narg == 1) cut_coul = cut_lj_global; - else cut_coul = utils::numeric(FLERR,arg[1],false,lmp); + cut_lj_global = utils::numeric(FLERR, arg[0], false, lmp); + if (narg == 1) + cut_coul = cut_lj_global; + else + cut_coul = 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; j <= atom->ntypes; j++) if (setflag[i][j]) cut_lj[i][j] = cut_lj_global; @@ -628,24 +632,23 @@ void PairLJClass2CoulLong::settings(int narg, char **arg) void PairLJClass2CoulLong::coeff(int narg, char **arg) { - if (narg < 4 || narg > 5) - error->all(FLERR,"Incorrect args for pair coefficients"); + if (narg < 4 || narg > 5) 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 epsilon_one = utils::numeric(FLERR,arg[2],false,lmp); - double sigma_one = utils::numeric(FLERR,arg[3],false,lmp); + double epsilon_one = utils::numeric(FLERR, arg[2], false, lmp); + double sigma_one = utils::numeric(FLERR, arg[3], false, lmp); - double cut_lj_one = cut_lj_global; - if (narg == 5) cut_lj_one = utils::numeric(FLERR,arg[4],false,lmp); + double cut_lj_one = cut_lj_global; + if (narg == 5) cut_lj_one = utils::numeric(FLERR, arg[4], 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++) { epsilon[i][j] = epsilon_one; sigma[i][j] = sigma_one; cut_lj[i][j] = cut_lj_one; @@ -654,7 +657,7 @@ void PairLJClass2CoulLong::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"); } /* ---------------------------------------------------------------------- @@ -663,45 +666,36 @@ void PairLJClass2CoulLong::coeff(int narg, char **arg) void PairLJClass2CoulLong::init_style() { - if (!atom->q_flag) - error->all(FLERR, - "Pair style lj/class2/coul/long requires atom attribute q"); + if (!atom->q_flag) error->all(FLERR, "Pair style lj/class2/coul/long requires atom attribute q"); // request regular or rRESPA neighbor list - int irequest; - int respa = 0; + int list_style = NeighConst::REQ_DEFAULT; - if (update->whichflag == 1 && utils::strmatch(update->integrate_style,"^respa")) { - if (((Respa *) update->integrate)->level_inner >= 0) respa = 1; - if (((Respa *) update->integrate)->level_middle >= 0) respa = 2; + if (update->whichflag == 1 && utils::strmatch(update->integrate_style, "^respa")) { + auto respa = (Respa *) update->integrate; + if (respa->level_inner >= 0) list_style = NeighConst::REQ_RESPA_INOUT; + if (respa->level_middle >= 0) list_style = NeighConst::REQ_RESPA_ALL; } - - irequest = neighbor->request(this,instance_me); - - if (respa >= 1) { - neighbor->requests[irequest]->respaouter = 1; - neighbor->requests[irequest]->respainner = 1; - } - if (respa == 2) neighbor->requests[irequest]->respamiddle = 1; + neighbor->add_request(this, list_style); cut_coulsq = cut_coul * cut_coul; // set rRESPA cutoffs - if (utils::strmatch(update->integrate_style,"^respa") && + if (utils::strmatch(update->integrate_style, "^respa") && ((Respa *) update->integrate)->level_inner >= 0) cut_respa = ((Respa *) update->integrate)->cutoff; - else cut_respa = nullptr; + else + cut_respa = nullptr; // insure use of KSpace long-range solver, set g_ewald - if (force->kspace == nullptr) - error->all(FLERR,"Pair style requires a KSpace style"); + if (force->kspace == nullptr) error->all(FLERR, "Pair style requires a KSpace style"); g_ewald = force->kspace->g_ewald; // setup force tables - if (ncoultablebits) init_tables(cut_coul,cut_respa); + if (ncoultablebits) init_tables(cut_coul, cut_respa); } /* ---------------------------------------------------------------------- @@ -714,26 +708,25 @@ double PairLJClass2CoulLong::init_one(int i, int j) // mix distance via user-defined rule if (setflag[i][j] == 0) { - epsilon[i][j] = 2.0 * sqrt(epsilon[i][i]*epsilon[j][j]) * - pow(sigma[i][i],3.0) * pow(sigma[j][j],3.0) / - (pow(sigma[i][i],6.0) + pow(sigma[j][j],6.0)); - sigma[i][j] = - pow((0.5 * (pow(sigma[i][i],6.0) + pow(sigma[j][j],6.0))),1.0/6.0); - cut_lj[i][j] = mix_distance(cut_lj[i][i],cut_lj[j][j]); + epsilon[i][j] = 2.0 * sqrt(epsilon[i][i] * epsilon[j][j]) * pow(sigma[i][i], 3.0) * + pow(sigma[j][j], 3.0) / (pow(sigma[i][i], 6.0) + pow(sigma[j][j], 6.0)); + sigma[i][j] = pow((0.5 * (pow(sigma[i][i], 6.0) + pow(sigma[j][j], 6.0))), 1.0 / 6.0); + cut_lj[i][j] = mix_distance(cut_lj[i][i], cut_lj[j][j]); } - double cut = MAX(cut_lj[i][j],cut_coul); + double cut = MAX(cut_lj[i][j], cut_coul); cut_ljsq[i][j] = cut_lj[i][j] * cut_lj[i][j]; - lj1[i][j] = 18.0 * epsilon[i][j] * pow(sigma[i][j],9.0); - lj2[i][j] = 18.0 * epsilon[i][j] * pow(sigma[i][j],6.0); - lj3[i][j] = 2.0 * epsilon[i][j] * pow(sigma[i][j],9.0); - lj4[i][j] = 3.0 * epsilon[i][j] * pow(sigma[i][j],6.0); + lj1[i][j] = 18.0 * epsilon[i][j] * pow(sigma[i][j], 9.0); + lj2[i][j] = 18.0 * epsilon[i][j] * pow(sigma[i][j], 6.0); + lj3[i][j] = 2.0 * epsilon[i][j] * pow(sigma[i][j], 9.0); + lj4[i][j] = 3.0 * epsilon[i][j] * pow(sigma[i][j], 6.0); if (offset_flag && (cut_lj[i][j] > 0.0)) { double ratio = sigma[i][j] / cut_lj[i][j]; - offset[i][j] = epsilon[i][j] * (2.0*pow(ratio,9.0) - 3.0*pow(ratio,6.0)); - } else offset[i][j] = 0.0; + offset[i][j] = epsilon[i][j] * (2.0 * pow(ratio, 9.0) - 3.0 * pow(ratio, 6.0)); + } else + offset[i][j] = 0.0; cut_ljsq[j][i] = cut_ljsq[i][j]; lj1[j][i] = lj1[i][j]; @@ -744,8 +737,8 @@ double PairLJClass2CoulLong::init_one(int i, int j) // check interior rRESPA cutoff - if (cut_respa && MIN(cut_lj[i][j],cut_coul) < cut_respa[3]) - error->all(FLERR,"Pair cutoff < Respa interior cutoff"); + if (cut_respa && MIN(cut_lj[i][j], cut_coul) < cut_respa[3]) + error->all(FLERR, "Pair cutoff < Respa interior cutoff"); // compute I,J contribution to long-range tail correction // count total # of atoms of type I and J via Allreduce @@ -754,22 +747,21 @@ double PairLJClass2CoulLong::init_one(int i, int j) int *type = atom->type; int nlocal = atom->nlocal; - double count[2],all[2]; + double count[2], all[2]; count[0] = count[1] = 0.0; for (int k = 0; k < nlocal; k++) { if (type[k] == i) count[0] += 1.0; if (type[k] == j) count[1] += 1.0; } - MPI_Allreduce(count,all,2,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(count, all, 2, MPI_DOUBLE, MPI_SUM, world); - double sig3 = sigma[i][j]*sigma[i][j]*sigma[i][j]; - double sig6 = sig3*sig3; - double rc3 = cut_lj[i][j]*cut_lj[i][j]*cut_lj[i][j]; - double rc6 = rc3*rc3; - etail_ij = 2.0*MY_PI*all[0]*all[1]*epsilon[i][j] * - sig6 * (sig3 - 3.0*rc3) / (3.0*rc6); - ptail_ij = 2.0*MY_PI*all[0]*all[1]*epsilon[i][j] * - sig6 * (sig3 - 2.0*rc3) / rc6; + double sig3 = sigma[i][j] * sigma[i][j] * sigma[i][j]; + double sig6 = sig3 * sig3; + double rc3 = cut_lj[i][j] * cut_lj[i][j] * cut_lj[i][j]; + double rc6 = rc3 * rc3; + etail_ij = + 2.0 * MY_PI * all[0] * all[1] * epsilon[i][j] * sig6 * (sig3 - 3.0 * rc3) / (3.0 * rc6); + ptail_ij = 2.0 * MY_PI * all[0] * all[1] * epsilon[i][j] * sig6 * (sig3 - 2.0 * rc3) / rc6; } return cut; @@ -783,14 +775,14 @@ void PairLJClass2CoulLong::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(&epsilon[i][j],sizeof(double),1,fp); - fwrite(&sigma[i][j],sizeof(double),1,fp); - fwrite(&cut_lj[i][j],sizeof(double),1,fp); + fwrite(&epsilon[i][j], sizeof(double), 1, fp); + fwrite(&sigma[i][j], sizeof(double), 1, fp); + fwrite(&cut_lj[i][j], sizeof(double), 1, fp); } } } @@ -804,21 +796,21 @@ void PairLJClass2CoulLong::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,&epsilon[i][j],sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&sigma[i][j],sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&cut_lj[i][j],sizeof(double),1,fp,nullptr,error); + utils::sfread(FLERR, &epsilon[i][j], sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &sigma[i][j], sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &cut_lj[i][j], sizeof(double), 1, fp, nullptr, error); } - MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&cut_lj[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&epsilon[i][j], 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&sigma[i][j], 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut_lj[i][j], 1, MPI_DOUBLE, 0, world); } } } @@ -829,13 +821,13 @@ void PairLJClass2CoulLong::read_restart(FILE *fp) void PairLJClass2CoulLong::write_restart_settings(FILE *fp) { - fwrite(&cut_lj_global,sizeof(double),1,fp); - fwrite(&cut_coul,sizeof(double),1,fp); - fwrite(&offset_flag,sizeof(int),1,fp); - fwrite(&mix_flag,sizeof(int),1,fp); - fwrite(&tail_flag,sizeof(int),1,fp); - fwrite(&ncoultablebits,sizeof(int),1,fp); - fwrite(&tabinner,sizeof(double),1,fp); + fwrite(&cut_lj_global, sizeof(double), 1, fp); + fwrite(&cut_coul, sizeof(double), 1, fp); + fwrite(&offset_flag, sizeof(int), 1, fp); + fwrite(&mix_flag, sizeof(int), 1, fp); + fwrite(&tail_flag, sizeof(int), 1, fp); + fwrite(&ncoultablebits, sizeof(int), 1, fp); + fwrite(&tabinner, sizeof(double), 1, fp); } /* ---------------------------------------------------------------------- @@ -845,21 +837,21 @@ void PairLJClass2CoulLong::write_restart_settings(FILE *fp) void PairLJClass2CoulLong::read_restart_settings(FILE *fp) { if (comm->me == 0) { - utils::sfread(FLERR,&cut_lj_global,sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&cut_coul,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,&tail_flag,sizeof(int),1,fp,nullptr,error); - utils::sfread(FLERR,&ncoultablebits,sizeof(int),1,fp,nullptr,error); - utils::sfread(FLERR,&tabinner,sizeof(double),1,fp,nullptr,error); + utils::sfread(FLERR, &cut_lj_global, sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &cut_coul, 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, &tail_flag, sizeof(int), 1, fp, nullptr, error); + utils::sfread(FLERR, &ncoultablebits, sizeof(int), 1, fp, nullptr, error); + utils::sfread(FLERR, &tabinner, sizeof(double), 1, fp, nullptr, error); } - MPI_Bcast(&cut_lj_global,1,MPI_DOUBLE,0,world); - MPI_Bcast(&cut_coul,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(&tail_flag,1,MPI_INT,0,world); - MPI_Bcast(&ncoultablebits,1,MPI_INT,0,world); - MPI_Bcast(&tabinner,1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut_lj_global, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut_coul, 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(&tail_flag, 1, MPI_INT, 0, world); + MPI_Bcast(&ncoultablebits, 1, MPI_INT, 0, world); + MPI_Bcast(&tabinner, 1, MPI_DOUBLE, 0, world); } /* ---------------------------------------------------------------------- @@ -868,8 +860,7 @@ void PairLJClass2CoulLong::read_restart_settings(FILE *fp) void PairLJClass2CoulLong::write_data(FILE *fp) { - for (int i = 1; i <= atom->ntypes; i++) - fprintf(fp,"%d %g %g\n",i,epsilon[i][i],sigma[i][i]); + for (int i = 1; i <= atom->ntypes; i++) fprintf(fp, "%d %g %g\n", i, epsilon[i][i], sigma[i][i]); } /* ---------------------------------------------------------------------- @@ -880,69 +871,68 @@ void PairLJClass2CoulLong::write_data_all(FILE *fp) { for (int i = 1; i <= atom->ntypes; i++) for (int j = i; j <= atom->ntypes; j++) - fprintf(fp,"%d %d %g %g %g\n",i,j,epsilon[i][j],sigma[i][j],cut_lj[i][j]); + fprintf(fp, "%d %d %g %g %g\n", i, j, epsilon[i][j], sigma[i][j], cut_lj[i][j]); } /* ---------------------------------------------------------------------- */ -double PairLJClass2CoulLong::single(int i, int j, int itype, int jtype, - double rsq, - double factor_coul, double factor_lj, - double &fforce) +double PairLJClass2CoulLong::single(int i, int j, int itype, int jtype, double rsq, + double factor_coul, double factor_lj, double &fforce) { - double r2inv,r,rinv,r3inv,r6inv,grij,expm2,t,erfc,prefactor; - double fraction,table,forcecoul,forcelj,phicoul,philj; + double r2inv, r, rinv, r3inv, r6inv, grij, expm2, t, erfc, prefactor; + double fraction, table, forcecoul, forcelj, phicoul, philj; int itable; - r2inv = 1.0/rsq; + r2inv = 1.0 / rsq; if (rsq < cut_coulsq) { if (!ncoultablebits || rsq <= tabinnersq) { r = sqrt(rsq); grij = g_ewald * r; - expm2 = exp(-grij*grij); - t = 1.0 / (1.0 + EWALD_P*grij); - erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2; - prefactor = force->qqrd2e * atom->q[i]*atom->q[j]/r; - forcecoul = prefactor * (erfc + EWALD_F*grij*expm2); - if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; + expm2 = exp(-grij * grij); + t = 1.0 / (1.0 + EWALD_P * grij); + erfc = t * (A1 + t * (A2 + t * (A3 + t * (A4 + t * A5)))) * expm2; + prefactor = force->qqrd2e * atom->q[i] * atom->q[j] / r; + forcecoul = prefactor * (erfc + EWALD_F * grij * expm2); + if (factor_coul < 1.0) forcecoul -= (1.0 - factor_coul) * prefactor; } else { union_int_float_t rsq_lookup; rsq_lookup.f = rsq; itable = rsq_lookup.i & ncoulmask; itable >>= ncoulshiftbits; fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable]; - table = ftable[itable] + fraction*dftable[itable]; - forcecoul = atom->q[i]*atom->q[j] * table; + table = ftable[itable] + fraction * dftable[itable]; + forcecoul = atom->q[i] * atom->q[j] * table; if (factor_coul < 1.0) { - table = ctable[itable] + fraction*dctable[itable]; - prefactor = atom->q[i]*atom->q[j] * table; - forcecoul -= (1.0-factor_coul)*prefactor; + table = ctable[itable] + fraction * dctable[itable]; + prefactor = atom->q[i] * atom->q[j] * table; + forcecoul -= (1.0 - factor_coul) * prefactor; } } - } else forcecoul = 0.0; + } else + forcecoul = 0.0; if (rsq < cut_ljsq[itype][jtype]) { rinv = sqrt(r2inv); - r3inv = r2inv*rinv; - r6inv = r3inv*r3inv; - forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]); - } else forcelj = 0.0; - fforce = (forcecoul + factor_lj*forcelj) * r2inv; + r3inv = r2inv * rinv; + r6inv = r3inv * r3inv; + forcelj = r6inv * (lj1[itype][jtype] * r3inv - lj2[itype][jtype]); + } else + forcelj = 0.0; + fforce = (forcecoul + factor_lj * forcelj) * r2inv; double eng = 0.0; if (rsq < cut_coulsq) { if (!ncoultablebits || rsq <= tabinnersq) - phicoul = prefactor*erfc; + phicoul = prefactor * erfc; else { - table = etable[itable] + fraction*detable[itable]; - phicoul = atom->q[i]*atom->q[j] * table; + table = etable[itable] + fraction * detable[itable]; + phicoul = atom->q[i] * atom->q[j] * table; } - if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor; + if (factor_coul < 1.0) phicoul -= (1.0 - factor_coul) * prefactor; eng += phicoul; } if (rsq < cut_ljsq[itype][jtype]) { - philj = r6inv*(lj3[itype][jtype]*r3inv-lj4[itype][jtype]) - - offset[itype][jtype]; - eng += factor_lj*philj; + philj = r6inv * (lj3[itype][jtype] * r3inv - lj4[itype][jtype]) - offset[itype][jtype]; + eng += factor_lj * philj; } return eng; @@ -953,9 +943,9 @@ double PairLJClass2CoulLong::single(int i, int j, int itype, int jtype, void *PairLJClass2CoulLong::extract(const char *str, int &dim) { dim = 0; - if (strcmp(str,"cut_coul") == 0) return (void *) &cut_coul; + if (strcmp(str, "cut_coul") == 0) return (void *) &cut_coul; dim = 2; - if (strcmp(str,"epsilon") == 0) return (void *) epsilon; - if (strcmp(str,"sigma") == 0) return (void *) sigma; + if (strcmp(str, "epsilon") == 0) return (void *) epsilon; + if (strcmp(str, "sigma") == 0) return (void *) sigma; return nullptr; } diff --git a/src/COLLOID/pair_brownian.cpp b/src/COLLOID/pair_brownian.cpp index 0640ae3e73..5782586654 100644 --- a/src/COLLOID/pair_brownian.cpp +++ b/src/COLLOID/pair_brownian.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -18,26 +17,26 @@ #include "pair_brownian.h" -#include -#include #include "atom.h" #include "comm.h" -#include "force.h" -#include "neighbor.h" -#include "neigh_list.h" #include "domain.h" -#include "update.h" -#include "modify.h" +#include "error.h" #include "fix.h" #include "fix_wall.h" +#include "force.h" #include "input.h" -#include "variable.h" -#include "random_mars.h" #include "math_const.h" #include "math_special.h" #include "memory.h" -#include "error.h" +#include "modify.h" +#include "neigh_list.h" +#include "neighbor.h" +#include "random_mars.h" +#include "update.h" +#include "variable.h" +#include +#include using namespace LAMMPS_NS; using namespace MathConst; @@ -45,7 +44,7 @@ using namespace MathSpecial; // same as fix_wall.cpp -enum{EDGE,CONSTANT,VARIABLE}; +enum { EDGE, CONSTANT, VARIABLE }; /* ---------------------------------------------------------------------- */ @@ -73,12 +72,12 @@ PairBrownian::~PairBrownian() void PairBrownian::compute(int eflag, int vflag) { - int i,j,ii,jj,inum,jnum,itype,jtype; - double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz,tx,ty,tz; - double rsq,r,h_sep,radi; - int *ilist,*jlist,*numneigh,**firstneigh; + int i, j, ii, jj, inum, jnum, itype, jtype; + double xtmp, ytmp, ztmp, delx, dely, delz, fx, fy, fz, tx, ty, tz; + double rsq, r, h_sep, radi; + int *ilist, *jlist, *numneigh, **firstneigh; - ev_init(eflag,vflag); + ev_init(eflag, vflag); double **x = atom->x; double **f = atom->f; @@ -91,19 +90,18 @@ void PairBrownian::compute(int eflag, int vflag) double vxmu2f = force->vxmu2f; double randr; double prethermostat; - double xl[3],a_sq,a_sh,a_pu,Fbmag; - double p1[3],p2[3],p3[3]; + double xl[3], a_sq, a_sh, a_pu, Fbmag; + double p1[3], p2[3], p3[3]; int overlaps = 0; // This section of code adjusts R0/RT0/RS0 if necessary due to changes // in the volume fraction as a result of fix deform or moving walls double dims[3], wallcoord; - if (flagVF) // Flag for volume fraction corrections - if (flagdeform || flagwall == 2) { // Possible changes in volume fraction + if (flagVF) // Flag for volume fraction corrections + if (flagdeform || flagwall == 2) { // Possible changes in volume fraction if (flagdeform && !flagwall) - for (j = 0; j < 3; j++) - dims[j] = domain->prd[j]; + for (j = 0; j < 3; j++) dims[j] = domain->prd[j]; else if (flagwall == 2 || (flagdeform && flagwall == 1)) { double wallhi[3], walllo[3]; for (int j = 0; j < 3; j++) { @@ -115,31 +113,32 @@ void PairBrownian::compute(int eflag, int vflag) int side = wallfix->wallwhich[m] % 2; if (wallfix->xstyle[m] == VARIABLE) { wallcoord = input->variable->compute_equal(wallfix->xindex[m]); - } - else wallcoord = wallfix->coord0[m]; - if (side == 0) walllo[dim] = wallcoord; - else wallhi[dim] = wallcoord; + } else + wallcoord = wallfix->coord0[m]; + if (side == 0) + walllo[dim] = wallcoord; + else + wallhi[dim] = wallcoord; } - for (int j = 0; j < 3; j++) - dims[j] = wallhi[j] - walllo[j]; + for (int j = 0; j < 3; j++) dims[j] = wallhi[j] - walllo[j]; } - double vol_T = dims[0]*dims[1]*dims[2]; - double vol_f = vol_P/vol_T; + double vol_T = dims[0] * dims[1] * dims[2]; + double vol_f = vol_P / vol_T; if (flaglog == 0) { - R0 = 6*MY_PI*mu*rad*(1.0 + 2.16*vol_f); - RT0 = 8*MY_PI*mu*cube(rad); + R0 = 6 * MY_PI * mu * rad * (1.0 + 2.16 * vol_f); + RT0 = 8 * MY_PI * mu * cube(rad); //RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3)*(1.0 + 3.33*vol_f + 2.80*vol_f*vol_f); } else { - R0 = 6*MY_PI*mu*rad*(1.0 + 2.725*vol_f - 6.583*vol_f*vol_f); - RT0 = 8*MY_PI*mu*cube(rad)*(1.0 + 0.749*vol_f - 2.469*vol_f*vol_f); + R0 = 6 * MY_PI * mu * rad * (1.0 + 2.725 * vol_f - 6.583 * vol_f * vol_f); + RT0 = 8 * MY_PI * mu * cube(rad) * (1.0 + 0.749 * vol_f - 2.469 * vol_f * vol_f); //RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3)*(1.0 + 3.64*vol_f - 6.95*vol_f*vol_f); } } // scale factor for Brownian moments - prethermostat = sqrt(24.0*force->boltz*t_target/update->dt); - prethermostat *= sqrt(force->vxmu2f/force->ftm2v/force->mvv2e); + prethermostat = sqrt(24.0 * force->boltz * t_target / update->dt); + prethermostat *= sqrt(force->vxmu2f / force->ftm2v / force->mvv2e); inum = list->inum; ilist = list->ilist; @@ -159,13 +158,13 @@ void PairBrownian::compute(int eflag, int vflag) // FLD contribution to force and torque due to isotropic terms if (flagfld) { - f[i][0] += prethermostat*sqrt(R0)*(random->uniform()-0.5); - f[i][1] += prethermostat*sqrt(R0)*(random->uniform()-0.5); - f[i][2] += prethermostat*sqrt(R0)*(random->uniform()-0.5); + f[i][0] += prethermostat * sqrt(R0) * (random->uniform() - 0.5); + f[i][1] += prethermostat * sqrt(R0) * (random->uniform() - 0.5); + f[i][2] += prethermostat * sqrt(R0) * (random->uniform() - 0.5); if (flaglog) { - torque[i][0] += prethermostat*sqrt(RT0)*(random->uniform()-0.5); - torque[i][1] += prethermostat*sqrt(RT0)*(random->uniform()-0.5); - torque[i][2] += prethermostat*sqrt(RT0)*(random->uniform()-0.5); + torque[i][0] += prethermostat * sqrt(RT0) * (random->uniform() - 0.5); + torque[i][1] += prethermostat * sqrt(RT0) * (random->uniform() - 0.5); + torque[i][2] += prethermostat * sqrt(RT0) * (random->uniform() - 0.5); } } @@ -178,7 +177,7 @@ void PairBrownian::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]; if (rsq < cutsq[itype][jtype]) { @@ -186,7 +185,7 @@ void PairBrownian::compute(int eflag, int vflag) // scalar resistances a_sq and a_sh - h_sep = r - 2.0*radi; + h_sep = r - 2.0 * radi; // check for overlaps @@ -194,35 +193,34 @@ void PairBrownian::compute(int eflag, int vflag) // if less than minimum gap, use minimum gap instead - if (r < cut_inner[itype][jtype]) - h_sep = cut_inner[itype][jtype] - 2.0*radi; + if (r < cut_inner[itype][jtype]) h_sep = cut_inner[itype][jtype] - 2.0 * radi; // scale h_sep by radi - h_sep = h_sep/radi; + h_sep = h_sep / radi; // scalar resistances if (flaglog) { - a_sq = 6.0*MY_PI*mu*radi*(1.0/4.0/h_sep + 9.0/40.0*log(1.0/h_sep)); - a_sh = 6.0*MY_PI*mu*radi*(1.0/6.0*log(1.0/h_sep)); - a_pu = 8.0*MY_PI*mu*cube(radi)*(3.0/160.0*log(1.0/h_sep)); + a_sq = 6.0 * MY_PI * mu * radi * (1.0 / 4.0 / h_sep + 9.0 / 40.0 * log(1.0 / h_sep)); + a_sh = 6.0 * MY_PI * mu * radi * (1.0 / 6.0 * log(1.0 / h_sep)); + a_pu = 8.0 * MY_PI * mu * cube(radi) * (3.0 / 160.0 * log(1.0 / h_sep)); } else - a_sq = 6.0*MY_PI*mu*radi*(1.0/4.0/h_sep); + a_sq = 6.0 * MY_PI * mu * radi * (1.0 / 4.0 / h_sep); // generate the Pairwise Brownian Force: a_sq - Fbmag = prethermostat*sqrt(a_sq); + Fbmag = prethermostat * sqrt(a_sq); // generate a random number - randr = random->uniform()-0.5; + randr = random->uniform() - 0.5; // contribution due to Brownian motion - fx = Fbmag*randr*delx/r; - fy = Fbmag*randr*dely/r; - fz = Fbmag*randr*delz/r; + fx = Fbmag * randr * delx / r; + fy = Fbmag * randr * dely / r; + fz = Fbmag * randr * delz / r; // add terms due to a_sh @@ -230,31 +228,33 @@ void PairBrownian::compute(int eflag, int vflag) // generate two orthogonal vectors to the line of centers - p1[0] = delx/r; p1[1] = dely/r; p1[2] = delz/r; - set_3_orthogonal_vectors(p1,p2,p3); + p1[0] = delx / r; + p1[1] = dely / r; + p1[2] = delz / r; + set_3_orthogonal_vectors(p1, p2, p3); // magnitude - Fbmag = prethermostat*sqrt(a_sh); + Fbmag = prethermostat * sqrt(a_sh); // force in each of the two directions - randr = random->uniform()-0.5; - fx += Fbmag*randr*p2[0]; - fy += Fbmag*randr*p2[1]; - fz += Fbmag*randr*p2[2]; + randr = random->uniform() - 0.5; + fx += Fbmag * randr * p2[0]; + fy += Fbmag * randr * p2[1]; + fz += Fbmag * randr * p2[2]; - randr = random->uniform()-0.5; - fx += Fbmag*randr*p3[0]; - fy += Fbmag*randr*p3[1]; - fz += Fbmag*randr*p3[2]; + randr = random->uniform() - 0.5; + fx += Fbmag * randr * p3[0]; + fy += Fbmag * randr * p3[1]; + fz += Fbmag * randr * p3[2]; } // scale forces to appropriate units - fx = vxmu2f*fx; - fy = vxmu2f*fy; - fz = vxmu2f*fz; + fx = vxmu2f * fx; + fy = vxmu2f * fy; + fz = vxmu2f * fz; // sum to total force @@ -279,15 +279,15 @@ void PairBrownian::compute(int eflag, int vflag) // location of the point of closest approach on I from its center - xl[0] = -delx/r*radi; - xl[1] = -dely/r*radi; - xl[2] = -delz/r*radi; + xl[0] = -delx / r * radi; + xl[1] = -dely / r * radi; + xl[2] = -delz / r * radi; // torque = xl_cross_F - tx = xl[1]*fz - xl[2]*fy; - ty = xl[2]*fx - xl[0]*fz; - tz = xl[0]*fy - xl[1]*fx; + tx = xl[1] * fz - xl[2] * fy; + ty = xl[2] * fx - xl[0] * fz; + tz = xl[0] * fy - xl[1] * fx; // torque is same on both particles @@ -303,19 +303,19 @@ void PairBrownian::compute(int eflag, int vflag) // torque due to a_pu - Fbmag = prethermostat*sqrt(a_pu); + Fbmag = prethermostat * sqrt(a_pu); // force in each direction - randr = random->uniform()-0.5; - tx = Fbmag*randr*p2[0]; - ty = Fbmag*randr*p2[1]; - tz = Fbmag*randr*p2[2]; + randr = random->uniform() - 0.5; + tx = Fbmag * randr * p2[0]; + ty = Fbmag * randr * p2[1]; + tz = Fbmag * randr * p2[2]; - randr = random->uniform()-0.5; - tx += Fbmag*randr*p3[0]; - ty += Fbmag*randr*p3[1]; - tz += Fbmag*randr*p3[2]; + randr = random->uniform() - 0.5; + tx += Fbmag * randr * p3[0]; + ty += Fbmag * randr * p3[1]; + tz += Fbmag * randr * p3[2]; // torque has opposite sign on two particles @@ -330,15 +330,14 @@ void PairBrownian::compute(int eflag, int vflag) } } - if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair, - 0.0,0.0,-fx,-fy,-fz,delx,dely,delz); + if (evflag) + ev_tally_xyz(i, j, nlocal, newton_pair, 0.0, 0.0, -fx, -fy, -fz, delx, dely, delz); } } } int print_overlaps = 0; - if (print_overlaps && overlaps) - printf("Number of overlaps=%d\n",overlaps); + if (print_overlaps && overlaps) printf("Number of overlaps=%d\n", overlaps); if (vflag_fdotr) virial_fdotr_compute(); } @@ -350,17 +349,16 @@ void PairBrownian::compute(int eflag, int vflag) void PairBrownian::allocate() { allocated = 1; - int n = atom->ntypes; + int np1 = atom->ntypes + 1; - 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; + memory->create(setflag, np1, np1, "pair:setflag"); + for (int i = 1; i < np1; i++) + for (int j = i; j < np1; j++) setflag[i][j] = 0; - memory->create(cutsq,n+1,n+1,"pair:cutsq"); + memory->create(cutsq, np1, np1, "pair:cutsq"); - memory->create(cut,n+1,n+1,"pair:cut"); - memory->create(cut_inner,n+1,n+1,"pair:cut_inner"); + memory->create(cut, np1, np1, "pair:cut"); + memory->create(cut_inner, np1, np1, "pair:cut_inner"); } /* ---------------------------------------------------------------------- @@ -369,24 +367,25 @@ void PairBrownian::allocate() void PairBrownian::settings(int narg, char **arg) { - if (narg != 7 && narg != 9) error->all(FLERR,"Illegal pair_style command"); + if (narg != 7 && narg != 9) error->all(FLERR, "Illegal pair_style command"); - mu = utils::numeric(FLERR,arg[0],false,lmp); - flaglog = utils::inumeric(FLERR,arg[1],false,lmp); - flagfld = utils::inumeric(FLERR,arg[2],false,lmp); - cut_inner_global = utils::numeric(FLERR,arg[3],false,lmp); - cut_global = utils::numeric(FLERR,arg[4],false,lmp); - t_target = utils::numeric(FLERR,arg[5],false,lmp); - seed = utils::inumeric(FLERR,arg[6],false,lmp); + mu = utils::numeric(FLERR, arg[0], false, lmp); + flaglog = utils::inumeric(FLERR, arg[1], false, lmp); + flagfld = utils::inumeric(FLERR, arg[2], false, lmp); + cut_inner_global = utils::numeric(FLERR, arg[3], false, lmp); + cut_global = utils::numeric(FLERR, arg[4], false, lmp); + t_target = utils::numeric(FLERR, arg[5], false, lmp); + seed = utils::inumeric(FLERR, arg[6], false, lmp); flagHI = flagVF = 1; if (narg == 9) { - flagHI = utils::inumeric(FLERR,arg[7],false,lmp); - flagVF = utils::inumeric(FLERR,arg[8],false,lmp); + flagHI = utils::inumeric(FLERR, arg[7], false, lmp); + flagVF = utils::inumeric(FLERR, arg[8], false, lmp); } if (flaglog == 1 && flagHI == 0) { - error->warning(FLERR,"Cannot include log terms without 1/r terms; " + error->warning(FLERR, + "Cannot include log terms without 1/r terms; " "setting flagHI to 1"); flagHI = 1; } @@ -394,7 +393,7 @@ void PairBrownian::settings(int narg, char **arg) // initialize Marsaglia RNG with processor-unique seed delete random; - random = new RanMars(lmp,seed + comm->me); + random = new RanMars(lmp, seed + comm->me); // reset cutoffs that have been explicitly set @@ -414,33 +413,32 @@ void PairBrownian::settings(int narg, char **arg) void PairBrownian::coeff(int narg, char **arg) { - if (narg != 2 && narg != 4) - error->all(FLERR,"Incorrect args for pair coefficients"); + if (narg != 2 && 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 cut_inner_one = cut_inner_global; double cut_one = cut_global; if (narg == 4) { - cut_inner_one = utils::numeric(FLERR,arg[2],false,lmp); - cut_one = utils::numeric(FLERR,arg[3],false,lmp); + cut_inner_one = utils::numeric(FLERR, arg[2], false, lmp); + 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++) { cut_inner[i][j] = cut_inner_one; cut[i][j] = cut_one; setflag[i][j] = 1; count++; } - if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); + if (count == 0) error->all(FLERR, "Incorrect args for pair coefficients"); } /* ---------------------------------------------------------------------- @@ -449,18 +447,15 @@ void PairBrownian::coeff(int narg, char **arg) void PairBrownian::init_style() { - if (!atom->sphere_flag) - error->all(FLERR,"Pair brownian requires atom style sphere"); + if (!atom->sphere_flag) error->all(FLERR, "Pair brownian requires atom style sphere"); // if newton off, forces between atoms ij will be double computed // using different random numbers if (force->newton_pair == 0 && comm->me == 0) - error->warning(FLERR, - "Pair brownian needs newton pair on for " - "momentum conservation"); + error->warning(FLERR, "Pair brownian needs newton pair on for momentum conservation"); - neighbor->request(this,instance_me); + neighbor->add_request(this); // insure all particles are finite-size // for pair hybrid, should limit test to types using the pair style @@ -469,17 +464,15 @@ void PairBrownian::init_style() int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) - if (radius[i] == 0.0) - error->one(FLERR,"Pair brownian requires extended particles"); + if (radius[i] == 0.0) error->one(FLERR, "Pair brownian requires extended particles"); // require monodisperse system with same radii for all types double radtype; for (int i = 1; i <= atom->ntypes; i++) { - if (!atom->radius_consistency(i,radtype)) - error->all(FLERR,"Pair brownian requires monodisperse particles"); - if (i > 1 && radtype != rad) - error->all(FLERR,"Pair brownian requires monodisperse particles"); + if (!atom->radius_consistency(i, radtype)) + error->all(FLERR, "Pair brownian requires monodisperse particles"); + if (i > 1 && radtype != rad) error->all(FLERR, "Pair brownian requires monodisperse particles"); rad = radtype; } @@ -496,15 +489,13 @@ void PairBrownian::init_style() flagdeform = flagwall = 0; for (int i = 0; i < modify->nfix; i++) { - if (strcmp(modify->fix[i]->style,"deform") == 0) + if (strcmp(modify->fix[i]->style, "deform") == 0) flagdeform = 1; - else if (strstr(modify->fix[i]->style,"wall") != nullptr) { - if (flagwall) - error->all(FLERR, - "Cannot use multiple fix wall commands with pair brownian"); - flagwall = 1; // Walls exist + else if (strstr(modify->fix[i]->style, "wall") != nullptr) { + if (flagwall) error->all(FLERR, "Cannot use multiple fix wall commands with pair brownian"); + flagwall = 1; // Walls exist wallfix = (FixWall *) modify->fix[i]; - if (wallfix->xflag) flagwall = 2; // Moving walls exist + if (wallfix->xflag) flagwall = 2; // Moving walls exist } } @@ -512,7 +503,8 @@ void PairBrownian::init_style() // vol_T = total volumeshearing = flagdeform = flagwall = 0; double vol_T, wallcoord; - if (!flagwall) vol_T = domain->xprd*domain->yprd*domain->zprd; + if (!flagwall) + vol_T = domain->xprd * domain->yprd * domain->zprd; else { double wallhi[3], walllo[3]; for (int j = 0; j < 3; j++) { @@ -528,31 +520,33 @@ void PairBrownian::init_style() wallcoord = input->variable->compute_equal(wallfix->xindex[m]); } - else wallcoord = wallfix->coord0[m]; + else + wallcoord = wallfix->coord0[m]; - if (side == 0) walllo[dim] = wallcoord; - else wallhi[dim] = wallcoord; + if (side == 0) + walllo[dim] = wallcoord; + else + wallhi[dim] = wallcoord; } - vol_T = (wallhi[0] - walllo[0]) * (wallhi[1] - walllo[1]) * - (wallhi[2] - walllo[2]); + vol_T = (wallhi[0] - walllo[0]) * (wallhi[1] - walllo[1]) * (wallhi[2] - walllo[2]); } // vol_P = volume of particles, assuming mono-dispersity // vol_f = volume fraction - vol_P = atom->natoms*(4.0/3.0)*MY_PI*cube(rad); + vol_P = atom->natoms * (4.0 / 3.0) * MY_PI * cube(rad); - double vol_f = vol_P/vol_T; + double vol_f = vol_P / vol_T; // set isotropic constants if (!flagVF) vol_f = 0; if (flaglog == 0) { - R0 = 6*MY_PI*mu*rad*(1.0 + 2.16*vol_f); - RT0 = 8*MY_PI*mu*cube(rad); // not actually needed + R0 = 6 * MY_PI * mu * rad * (1.0 + 2.16 * vol_f); + RT0 = 8 * MY_PI * mu * cube(rad); // not actually needed } else { - R0 = 6*MY_PI*mu*rad*(1.0 + 2.725*vol_f - 6.583*vol_f*vol_f); - RT0 = 8*MY_PI*mu*cube(rad)*(1.0 + 0.749*vol_f - 2.469*vol_f*vol_f); + R0 = 6 * MY_PI * mu * rad * (1.0 + 2.725 * vol_f - 6.583 * vol_f * vol_f); + RT0 = 8 * MY_PI * mu * cube(rad) * (1.0 + 0.749 * vol_f - 2.469 * vol_f * vol_f); } } @@ -563,8 +557,8 @@ void PairBrownian::init_style() double PairBrownian::init_one(int i, int j) { if (setflag[i][j] == 0) { - cut_inner[i][j] = mix_distance(cut_inner[i][i],cut_inner[j][j]); - cut[i][j] = mix_distance(cut[i][i],cut[j][j]); + cut_inner[i][j] = mix_distance(cut_inner[i][i], cut_inner[j][j]); + cut[i][j] = mix_distance(cut[i][i], cut[j][j]); } cut_inner[j][i] = cut_inner[i][j]; @@ -580,13 +574,13 @@ void PairBrownian::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(&cut_inner[i][j],sizeof(double),1,fp); - fwrite(&cut[i][j],sizeof(double),1,fp); + fwrite(&cut_inner[i][j], sizeof(double), 1, fp); + fwrite(&cut[i][j], sizeof(double), 1, fp); } } } @@ -600,19 +594,19 @@ void PairBrownian::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,&cut_inner[i][j],sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&cut[i][j],sizeof(double),1,fp,nullptr,error); + utils::sfread(FLERR, &cut_inner[i][j], sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &cut[i][j], sizeof(double), 1, fp, nullptr, error); } - MPI_Bcast(&cut_inner[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut_inner[i][j], 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut[i][j], 1, MPI_DOUBLE, 0, world); } } } @@ -623,17 +617,17 @@ void PairBrownian::read_restart(FILE *fp) void PairBrownian::write_restart_settings(FILE *fp) { - fwrite(&mu,sizeof(double),1,fp); - fwrite(&flaglog,sizeof(int),1,fp); - fwrite(&flagfld,sizeof(int),1,fp); - fwrite(&cut_inner_global,sizeof(double),1,fp); - fwrite(&cut_global,sizeof(double),1,fp); - fwrite(&t_target,sizeof(double),1,fp); - fwrite(&seed,sizeof(int),1,fp); - fwrite(&offset_flag,sizeof(int),1,fp); - fwrite(&mix_flag,sizeof(int),1,fp); - fwrite(&flagHI,sizeof(int),1,fp); - fwrite(&flagVF,sizeof(int),1,fp); + fwrite(&mu, sizeof(double), 1, fp); + fwrite(&flaglog, sizeof(int), 1, fp); + fwrite(&flagfld, sizeof(int), 1, fp); + fwrite(&cut_inner_global, sizeof(double), 1, fp); + fwrite(&cut_global, sizeof(double), 1, fp); + fwrite(&t_target, sizeof(double), 1, fp); + fwrite(&seed, sizeof(int), 1, fp); + fwrite(&offset_flag, sizeof(int), 1, fp); + fwrite(&mix_flag, sizeof(int), 1, fp); + fwrite(&flagHI, sizeof(int), 1, fp); + fwrite(&flagVF, sizeof(int), 1, fp); } /* ---------------------------------------------------------------------- @@ -644,57 +638,56 @@ void PairBrownian::read_restart_settings(FILE *fp) { int me = comm->me; if (me == 0) { - utils::sfread(FLERR,&mu,sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&flaglog,sizeof(int),1,fp,nullptr,error); - utils::sfread(FLERR,&flagfld,sizeof(int),1,fp,nullptr,error); - utils::sfread(FLERR,&cut_inner_global,sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&cut_global,sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&t_target, sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&seed, sizeof(int),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,&flagHI,sizeof(int),1,fp,nullptr,error); - utils::sfread(FLERR,&flagVF,sizeof(int),1,fp,nullptr,error); + utils::sfread(FLERR, &mu, sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &flaglog, sizeof(int), 1, fp, nullptr, error); + utils::sfread(FLERR, &flagfld, sizeof(int), 1, fp, nullptr, error); + utils::sfread(FLERR, &cut_inner_global, sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &cut_global, sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &t_target, sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &seed, sizeof(int), 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, &flagHI, sizeof(int), 1, fp, nullptr, error); + utils::sfread(FLERR, &flagVF, sizeof(int), 1, fp, nullptr, error); } - MPI_Bcast(&mu,1,MPI_DOUBLE,0,world); - MPI_Bcast(&flaglog,1,MPI_INT,0,world); - MPI_Bcast(&flagfld,1,MPI_INT,0,world); - MPI_Bcast(&cut_inner_global,1,MPI_DOUBLE,0,world); - MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world); - MPI_Bcast(&t_target,1,MPI_DOUBLE,0,world); - MPI_Bcast(&seed,1,MPI_INT,0,world); - MPI_Bcast(&offset_flag,1,MPI_INT,0,world); - MPI_Bcast(&mix_flag,1,MPI_INT,0,world); - MPI_Bcast(&flagHI,1,MPI_INT,0,world); - MPI_Bcast(&flagVF,1,MPI_INT,0,world); + MPI_Bcast(&mu, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&flaglog, 1, MPI_INT, 0, world); + MPI_Bcast(&flagfld, 1, MPI_INT, 0, world); + MPI_Bcast(&cut_inner_global, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut_global, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&t_target, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&seed, 1, MPI_INT, 0, world); + MPI_Bcast(&offset_flag, 1, MPI_INT, 0, world); + MPI_Bcast(&mix_flag, 1, MPI_INT, 0, world); + MPI_Bcast(&flagHI, 1, MPI_INT, 0, world); + MPI_Bcast(&flagVF, 1, MPI_INT, 0, world); // additional setup based on restart parameters delete random; - random = new RanMars(lmp,seed + comm->me); + random = new RanMars(lmp, seed + comm->me); } /* ----------------------------------------------------------------------*/ -void PairBrownian::set_3_orthogonal_vectors(double p1[3], - double p2[3], double p3[3]) +void PairBrownian::set_3_orthogonal_vectors(double p1[3], double p2[3], double p3[3]) { double norm; - int ix,iy,iz; + int ix, iy, iz; // find the index of maximum magnitude and store it in iz if (fabs(p1[0]) > fabs(p1[1])) { - iz=0; - ix=1; - iy=2; + iz = 0; + ix = 1; + iy = 2; } else { - iz=1; - ix=2; - iy=0; + iz = 1; + ix = 2; + iy = 0; } - if (iz==0) { + if (iz == 0) { if (fabs(p1[0]) < fabs(p1[2])) { iz = 2; ix = 0; @@ -710,21 +703,21 @@ void PairBrownian::set_3_orthogonal_vectors(double p1[3], // set p2 arbitrarily such that it's orthogonal to p1 - p2[ix]=1.0; - p2[iy]=1.0; - p2[iz] = -(p1[ix]*p2[ix] + p1[iy]*p2[iy])/p1[iz]; + p2[ix] = 1.0; + p2[iy] = 1.0; + p2[iz] = -(p1[ix] * p2[ix] + p1[iy] * p2[iy]) / p1[iz]; // normalize p2 - norm = sqrt(p2[0]*p2[0] + p2[1]*p2[1] + p2[2]*p2[2]); + norm = sqrt(p2[0] * p2[0] + p2[1] * p2[1] + p2[2] * p2[2]); - p2[0] = p2[0]/norm; - p2[1] = p2[1]/norm; - p2[2] = p2[2]/norm; + p2[0] = p2[0] / norm; + p2[1] = p2[1] / norm; + p2[2] = p2[2] / norm; // Set p3 by taking the cross product p3=p2xp1 - p3[0] = p1[1]*p2[2] - p1[2]*p2[1]; - p3[1] = p1[2]*p2[0] - p1[0]*p2[2]; - p3[2] = p1[0]*p2[1] - p1[1]*p2[0]; + p3[0] = p1[1] * p2[2] - p1[2] * p2[1]; + p3[1] = p1[2] * p2[0] - p1[0] * p2[2]; + p3[2] = p1[0] * p2[1] - p1[1] * p2[0]; } diff --git a/src/COLLOID/pair_brownian_poly.cpp b/src/COLLOID/pair_brownian_poly.cpp index 0c3b9f1ed0..2edbbadc0f 100644 --- a/src/COLLOID/pair_brownian_poly.cpp +++ b/src/COLLOID/pair_brownian_poly.cpp @@ -19,24 +19,24 @@ #include "pair_brownian_poly.h" -#include -#include #include "atom.h" -#include "force.h" -#include "neighbor.h" -#include "neigh_list.h" -#include "neigh_request.h" #include "domain.h" -#include "update.h" -#include "modify.h" +#include "error.h" #include "fix.h" #include "fix_wall.h" +#include "force.h" #include "input.h" -#include "variable.h" -#include "random_mars.h" #include "math_const.h" #include "math_special.h" -#include "error.h" +#include "modify.h" +#include "neigh_list.h" +#include "neighbor.h" +#include "random_mars.h" +#include "update.h" +#include "variable.h" + +#include +#include using namespace LAMMPS_NS; using namespace MathConst; @@ -340,9 +340,7 @@ void PairBrownianPoly::init_style() if (radius[i] == 0.0) error->one(FLERR,"Pair brownian/poly requires extended particles"); - int irequest = neighbor->request(this,instance_me); - neighbor->requests[irequest]->half = 0; - neighbor->requests[irequest]->full = 1; + neighbor->add_request(this, NeighConst::REQ_FULL); // set the isotropic constants that depend on the volume fraction // vol_T = total volume diff --git a/src/COLLOID/pair_colloid.cpp b/src/COLLOID/pair_colloid.cpp index 2c9887e8a5..891aae9233 100644 --- a/src/COLLOID/pair_colloid.cpp +++ b/src/COLLOID/pair_colloid.cpp @@ -18,7 +18,6 @@ #include "pair_colloid.h" -#include #include "atom.h" #include "comm.h" #include "force.h" @@ -27,6 +26,7 @@ #include "memory.h" #include "error.h" +#include using namespace LAMMPS_NS; using namespace MathSpecial; diff --git a/src/COLLOID/pair_lubricate.cpp b/src/COLLOID/pair_lubricate.cpp index 526f12e54c..c2117cf1ad 100644 --- a/src/COLLOID/pair_lubricate.cpp +++ b/src/COLLOID/pair_lubricate.cpp @@ -19,24 +19,24 @@ #include "pair_lubricate.h" -#include -#include #include "atom.h" #include "comm.h" -#include "force.h" -#include "neighbor.h" -#include "neigh_list.h" #include "domain.h" -#include "modify.h" +#include "error.h" #include "fix.h" #include "fix_deform.h" #include "fix_wall.h" +#include "force.h" #include "input.h" -#include "variable.h" #include "math_const.h" #include "memory.h" -#include "error.h" +#include "modify.h" +#include "neigh_list.h" +#include "neighbor.h" +#include "variable.h" +#include +#include using namespace LAMMPS_NS; using namespace MathConst; @@ -536,7 +536,7 @@ void PairLubricate::init_style() if (comm->ghost_velocity == 0) error->all(FLERR,"Pair lubricate requires ghost atoms store velocity"); - neighbor->request(this,instance_me); + neighbor->add_request(this); // require that atom radii are identical within each type // require monodisperse system with same radii for all types diff --git a/src/COLLOID/pair_lubricateU.cpp b/src/COLLOID/pair_lubricateU.cpp index 499e998b85..3d4fbcbc7d 100644 --- a/src/COLLOID/pair_lubricateU.cpp +++ b/src/COLLOID/pair_lubricateU.cpp @@ -18,24 +18,24 @@ #include "pair_lubricateU.h" -#include -#include #include "atom.h" #include "comm.h" -#include "force.h" -#include "neighbor.h" -#include "neigh_list.h" #include "domain.h" -#include "update.h" -#include "math_const.h" -#include "modify.h" +#include "error.h" #include "fix.h" #include "fix_wall.h" +#include "force.h" #include "input.h" -#include "variable.h" +#include "math_const.h" #include "memory.h" -#include "error.h" +#include "modify.h" +#include "neigh_list.h" +#include "neighbor.h" +#include "update.h" +#include "variable.h" +#include +#include using namespace LAMMPS_NS; using namespace MathConst; @@ -1769,7 +1769,7 @@ void PairLubricateU::init_style() if (comm->ghost_velocity == 0) error->all(FLERR,"Pair lubricateU requires ghost atoms store velocity"); - neighbor->request(this,instance_me); + neighbor->add_request(this); // require that atom radii are identical within each type // require monodisperse system with same radii for all types diff --git a/src/COLLOID/pair_lubricateU_poly.cpp b/src/COLLOID/pair_lubricateU_poly.cpp index a778cf6d15..0aed0df97f 100644 --- a/src/COLLOID/pair_lubricateU_poly.cpp +++ b/src/COLLOID/pair_lubricateU_poly.cpp @@ -20,23 +20,23 @@ #include "pair_lubricateU_poly.h" -#include -#include #include "atom.h" #include "comm.h" -#include "force.h" -#include "neighbor.h" -#include "neigh_list.h" -#include "neigh_request.h" #include "domain.h" -#include "modify.h" +#include "error.h" #include "fix.h" #include "fix_wall.h" +#include "force.h" #include "input.h" -#include "variable.h" #include "math_const.h" #include "memory.h" -#include "error.h" +#include "modify.h" +#include "neigh_list.h" +#include "neighbor.h" +#include "variable.h" + +#include +#include using namespace LAMMPS_NS; using namespace MathConst; @@ -1235,7 +1235,5 @@ void PairLubricateUPoly::init_style() RS0 = 20.0/3.0*MY_PI*mu*(1.0 + 3.64*vol_f - 6.95*vol_f*vol_f); } - int irequest = neighbor->request(this,instance_me); - neighbor->requests[irequest]->half = 0; - neighbor->requests[irequest]->full = 1; + neighbor->add_request(this, NeighConst::REQ_FULL); } diff --git a/src/COLLOID/pair_lubricate_poly.cpp b/src/COLLOID/pair_lubricate_poly.cpp index 2d921debfe..97aa6bf6ff 100644 --- a/src/COLLOID/pair_lubricate_poly.cpp +++ b/src/COLLOID/pair_lubricate_poly.cpp @@ -20,23 +20,23 @@ #include "pair_lubricate_poly.h" -#include -#include #include "atom.h" #include "comm.h" -#include "force.h" -#include "neighbor.h" -#include "neigh_list.h" -#include "neigh_request.h" #include "domain.h" -#include "modify.h" +#include "error.h" #include "fix.h" #include "fix_deform.h" #include "fix_wall.h" +#include "force.h" #include "input.h" -#include "variable.h" #include "math_const.h" -#include "error.h" +#include "modify.h" +#include "neigh_list.h" +#include "neighbor.h" +#include "variable.h" + +#include +#include using namespace LAMMPS_NS; using namespace MathConst; @@ -445,9 +445,7 @@ void PairLubricatePoly::init_style() if (radius[i] == 0.0) error->one(FLERR,"Pair lubricate/poly requires extended particles"); - int irequest = neighbor->request(this,instance_me); - neighbor->requests[irequest]->half = 0; - neighbor->requests[irequest]->full = 1; + neighbor->add_request(this, NeighConst::REQ_FULL); // set the isotropic constants that depend on the volume fraction // vol_T = total volume diff --git a/src/COLLOID/pair_yukawa_colloid.cpp b/src/COLLOID/pair_yukawa_colloid.cpp index 323390169a..67b518fe70 100644 --- a/src/COLLOID/pair_yukawa_colloid.cpp +++ b/src/COLLOID/pair_yukawa_colloid.cpp @@ -124,14 +124,13 @@ void PairYukawaColloid::init_style() if (!atom->sphere_flag) error->all(FLERR,"Pair yukawa/colloid requires atom style sphere"); - neighbor->request(this,instance_me); + neighbor->add_request(this); // require that atom radii are identical within each type for (int i = 1; i <= atom->ntypes; i++) if (!atom->radius_consistency(i,rad[i])) - error->all(FLERR,"Pair yukawa/colloid requires atoms with same type " - "have same radius"); + error->all(FLERR,"Pair yukawa/colloid requires atoms with same type have same radius"); } /* ---------------------------------------------------------------------- diff --git a/src/DIELECTRIC/fix_polarize_functional.cpp b/src/DIELECTRIC/fix_polarize_functional.cpp index d631447d89..15140e7d58 100644 --- a/src/DIELECTRIC/fix_polarize_functional.cpp +++ b/src/DIELECTRIC/fix_polarize_functional.cpp @@ -40,7 +40,6 @@ #include "memory.h" #include "msm_dielectric.h" #include "neigh_list.h" -#include "neigh_request.h" #include "neighbor.h" #include "pair_coul_cut_dielectric.h" #include "pair_coul_long_dielectric.h" @@ -265,12 +264,7 @@ void FixPolarizeFunctional::init() // need a full neighbor list w/ Newton off and ghost neighbors // built whenever re-neighboring occurs - int irequest = neighbor->request(this, instance_me); - neighbor->requests[irequest]->pair = 0; - neighbor->requests[irequest]->fix = 1; - neighbor->requests[irequest]->half = 0; - neighbor->requests[irequest]->full = 1; - neighbor->requests[irequest]->occasional = 0; + neighbor->add_request(this, NeighConst::REQ_FULL | NeighConst::REQ_OCCASIONAL); if (force->kspace) g_ewald = force->kspace->g_ewald; diff --git a/src/DIELECTRIC/pair_coul_cut_dielectric.cpp b/src/DIELECTRIC/pair_coul_cut_dielectric.cpp index 9596ce5be2..7e3602967f 100644 --- a/src/DIELECTRIC/pair_coul_cut_dielectric.cpp +++ b/src/DIELECTRIC/pair_coul_cut_dielectric.cpp @@ -24,7 +24,6 @@ #include "math_const.h" #include "memory.h" #include "neigh_list.h" -#include "neigh_request.h" #include "neighbor.h" #include @@ -166,9 +165,7 @@ void PairCoulCutDielectric::init_style() avec = (AtomVecDielectric *) atom->style_match("dielectric"); if (!avec) error->all(FLERR, "Pair coul/cut/dielectric requires atom style dielectric"); - int irequest = neighbor->request(this, instance_me); - neighbor->requests[irequest]->half = 0; - neighbor->requests[irequest]->full = 1; + neighbor->add_request(this, NeighConst::REQ_FULL); } /* ---------------------------------------------------------------------- */ diff --git a/src/DIELECTRIC/pair_coul_long_dielectric.cpp b/src/DIELECTRIC/pair_coul_long_dielectric.cpp index 4bf7d48c2b..66a4473bfc 100644 --- a/src/DIELECTRIC/pair_coul_long_dielectric.cpp +++ b/src/DIELECTRIC/pair_coul_long_dielectric.cpp @@ -25,7 +25,6 @@ #include "math_const.h" #include "memory.h" #include "neigh_list.h" -#include "neigh_request.h" #include "neighbor.h" #include @@ -211,9 +210,7 @@ void PairCoulLongDielectric::init_style() avec = (AtomVecDielectric *) atom->style_match("dielectric"); if (!avec) error->all(FLERR, "Pair coul/long/dielectric requires atom style dielectric"); - int irequest = neighbor->request(this, instance_me); - neighbor->requests[irequest]->half = 0; - neighbor->requests[irequest]->full = 1; + neighbor->add_request(this, NeighConst::REQ_FULL); cut_coulsq = cut_coul * cut_coul; diff --git a/src/DIELECTRIC/pair_lj_cut_coul_cut_dielectric.cpp b/src/DIELECTRIC/pair_lj_cut_coul_cut_dielectric.cpp index 3e840459b7..94ecdac578 100644 --- a/src/DIELECTRIC/pair_lj_cut_coul_cut_dielectric.cpp +++ b/src/DIELECTRIC/pair_lj_cut_coul_cut_dielectric.cpp @@ -24,7 +24,6 @@ #include "math_const.h" #include "memory.h" #include "neigh_list.h" -#include "neigh_request.h" #include "neighbor.h" #include @@ -194,9 +193,7 @@ void PairLJCutCoulCutDielectric::init_style() avec = (AtomVecDielectric *) atom->style_match("dielectric"); if (!avec) error->all(FLERR, "Pair lj/cut/coul/cut/dielectric requires atom style dielectric"); - int irequest = neighbor->request(this, instance_me); - neighbor->requests[irequest]->half = 0; - neighbor->requests[irequest]->full = 1; + neighbor->add_request(this, NeighConst::REQ_FULL); } /* ---------------------------------------------------------------------- */ diff --git a/src/DIELECTRIC/pair_lj_cut_coul_debye_dielectric.cpp b/src/DIELECTRIC/pair_lj_cut_coul_debye_dielectric.cpp index 5220be7b62..8858a444a8 100644 --- a/src/DIELECTRIC/pair_lj_cut_coul_debye_dielectric.cpp +++ b/src/DIELECTRIC/pair_lj_cut_coul_debye_dielectric.cpp @@ -24,7 +24,6 @@ #include "math_const.h" #include "memory.h" #include "neigh_list.h" -#include "neigh_request.h" #include "neighbor.h" #include @@ -198,9 +197,7 @@ void PairLJCutCoulDebyeDielectric::init_style() avec = (AtomVecDielectric *) atom->style_match("dielectric"); if (!avec) error->all(FLERR, "Pair lj/cut/coul/debye/dielectric requires atom style dielectric"); - int irequest = neighbor->request(this, instance_me); - neighbor->requests[irequest]->half = 0; - neighbor->requests[irequest]->full = 1; + neighbor->add_request(this, NeighConst::REQ_FULL); } /* ---------------------------------------------------------------------- */ diff --git a/src/DIELECTRIC/pair_lj_cut_coul_long_dielectric.cpp b/src/DIELECTRIC/pair_lj_cut_coul_long_dielectric.cpp index b256dd5a6c..22a725d045 100644 --- a/src/DIELECTRIC/pair_lj_cut_coul_long_dielectric.cpp +++ b/src/DIELECTRIC/pair_lj_cut_coul_long_dielectric.cpp @@ -25,7 +25,6 @@ #include "math_const.h" #include "memory.h" #include "neigh_list.h" -#include "neigh_request.h" #include "neighbor.h" #include @@ -248,9 +247,7 @@ void PairLJCutCoulLongDielectric::init_style() avec = (AtomVecDielectric *) atom->style_match("dielectric"); if (!avec) error->all(FLERR, "Pair lj/cut/coul/long/dielectric requires atom style dielectric"); - int irequest = neighbor->request(this, instance_me); - neighbor->requests[irequest]->half = 0; - neighbor->requests[irequest]->full = 1; + neighbor->add_request(this, NeighConst::REQ_FULL); cut_coulsq = cut_coul * cut_coul; diff --git a/src/DIELECTRIC/pair_lj_cut_coul_msm_dielectric.cpp b/src/DIELECTRIC/pair_lj_cut_coul_msm_dielectric.cpp index b729abef8b..751f4d509a 100644 --- a/src/DIELECTRIC/pair_lj_cut_coul_msm_dielectric.cpp +++ b/src/DIELECTRIC/pair_lj_cut_coul_msm_dielectric.cpp @@ -25,7 +25,6 @@ #include "math_const.h" #include "memory.h" #include "neigh_list.h" -#include "neigh_request.h" #include "neighbor.h" #include @@ -356,9 +355,7 @@ void PairLJCutCoulMSMDielectric::init_style() avec = (AtomVecDielectric *) atom->style_match("dielectric"); if (!avec) error->all(FLERR, "Pair lj/cut/coul/msm/dielectric requires atom style dielectric"); - int irequest = neighbor->request(this, instance_me); - neighbor->requests[irequest]->half = 0; - neighbor->requests[irequest]->full = 1; + neighbor->add_request(this, NeighConst::REQ_FULL); cut_coulsq = cut_coul * cut_coul; diff --git a/src/DIELECTRIC/pair_lj_long_coul_long_dielectric.cpp b/src/DIELECTRIC/pair_lj_long_coul_long_dielectric.cpp index 25a533632b..3fd86587e9 100644 --- a/src/DIELECTRIC/pair_lj_long_coul_long_dielectric.cpp +++ b/src/DIELECTRIC/pair_lj_long_coul_long_dielectric.cpp @@ -25,7 +25,6 @@ #include "math_extra.h" #include "memory.h" #include "neigh_list.h" -#include "neigh_request.h" #include "neighbor.h" #include @@ -75,9 +74,7 @@ void PairLJLongCoulLongDielectric::init_style() avec = (AtomVecDielectric *) atom->style_match("dielectric"); if (!avec) error->all(FLERR, "Pair lj/long/coul/long/dielectric requires atom style dielectric"); - int irequest = neighbor->request(this, instance_me); - neighbor->requests[irequest]->half = 0; - neighbor->requests[irequest]->full = 1; + neighbor->add_request(this, NeighConst::REQ_FULL); } /* ---------------------------------------------------------------------- diff --git a/src/DIPOLE/pair_lj_cut_dipole_cut.cpp b/src/DIPOLE/pair_lj_cut_dipole_cut.cpp index 0f8a7317c6..6f0ed58520 100644 --- a/src/DIPOLE/pair_lj_cut_dipole_cut.cpp +++ b/src/DIPOLE/pair_lj_cut_dipole_cut.cpp @@ -14,17 +14,17 @@ #include "pair_lj_cut_dipole_cut.h" -#include -#include #include "atom.h" -#include "neighbor.h" -#include "neigh_list.h" #include "comm.h" +#include "error.h" #include "force.h" #include "memory.h" -#include "error.h" +#include "neigh_list.h" +#include "neighbor.h" #include "update.h" +#include +#include using namespace LAMMPS_NS; @@ -362,7 +362,7 @@ void PairLJCutDipoleCut::init_style() if (!atom->q_flag || !atom->mu_flag || !atom->torque_flag) error->all(FLERR,"Pair dipole/cut requires atom attributes q, mu, torque"); - neighbor->request(this,instance_me); + neighbor->add_request(this); } /* ---------------------------------------------------------------------- diff --git a/src/DIPOLE/pair_lj_cut_dipole_long.cpp b/src/DIPOLE/pair_lj_cut_dipole_long.cpp index b13e48b808..6d6ddb748a 100644 --- a/src/DIPOLE/pair_lj_cut_dipole_long.cpp +++ b/src/DIPOLE/pair_lj_cut_dipole_long.cpp @@ -14,19 +14,19 @@ #include "pair_lj_cut_dipole_long.h" -#include -#include #include "atom.h" #include "comm.h" -#include "neighbor.h" -#include "neigh_list.h" +#include "error.h" #include "force.h" #include "kspace.h" #include "math_const.h" #include "memory.h" -#include "error.h" +#include "neigh_list.h" +#include "neighbor.h" #include "update.h" +#include +#include using namespace LAMMPS_NS; using namespace MathConst; @@ -454,7 +454,7 @@ void PairLJCutDipoleLong::init_style() cut_coulsq = cut_coul * cut_coul; - neighbor->request(this,instance_me); + neighbor->add_request(this); } /* ---------------------------------------------------------------------- diff --git a/src/DIPOLE/pair_lj_long_dipole_long.cpp b/src/DIPOLE/pair_lj_long_dipole_long.cpp index 0d94d7116d..76dc487528 100644 --- a/src/DIPOLE/pair_lj_long_dipole_long.cpp +++ b/src/DIPOLE/pair_lj_long_dipole_long.cpp @@ -235,7 +235,7 @@ void PairLJLongDipoleLong::init_style() if (!atom->mu_flag || !atom->torque_flag) error->all(FLERR,"Pair lj/long/dipole/long requires atom attributes mu, torque"); - neighbor->request(this,instance_me); + neighbor->add_request(this); cut_coulsq = cut_coul * cut_coul; diff --git a/src/DIPOLE/pair_lj_sf_dipole_sf.cpp b/src/DIPOLE/pair_lj_sf_dipole_sf.cpp index 7d88513910..f32f50461c 100644 --- a/src/DIPOLE/pair_lj_sf_dipole_sf.cpp +++ b/src/DIPOLE/pair_lj_sf_dipole_sf.cpp @@ -37,12 +37,6 @@ static int warn_single = 0; /* ---------------------------------------------------------------------- */ -PairLJSFDipoleSF::PairLJSFDipoleSF(LAMMPS *lmp) : Pair(lmp) -{ -} - -/* ---------------------------------------------------------------------- */ - PairLJSFDipoleSF::~PairLJSFDipoleSF() { if (allocated) { @@ -418,7 +412,7 @@ void PairLJSFDipoleSF::init_style() if (!atom->q_flag || !atom->mu_flag || !atom->torque_flag) error->all(FLERR,"Pair dipole/sf requires atom attributes q, mu, torque"); - neighbor->request(this,instance_me); + neighbor->add_request(this); } /* ---------------------------------------------------------------------- diff --git a/src/DIPOLE/pair_lj_sf_dipole_sf.h b/src/DIPOLE/pair_lj_sf_dipole_sf.h index 22b9320d8a..809d352472 100644 --- a/src/DIPOLE/pair_lj_sf_dipole_sf.h +++ b/src/DIPOLE/pair_lj_sf_dipole_sf.h @@ -26,7 +26,7 @@ namespace LAMMPS_NS { class PairLJSFDipoleSF : public Pair { public: - PairLJSFDipoleSF(class LAMMPS *); + PairLJSFDipoleSF(class LAMMPS *_lmp) : Pair(_lmp) {}; ~PairLJSFDipoleSF() override; void compute(int, int) override; void settings(int, char **) override; diff --git a/src/DPD-BASIC/pair_dpd.cpp b/src/DPD-BASIC/pair_dpd.cpp index dec016706e..caa8161573 100644 --- a/src/DPD-BASIC/pair_dpd.cpp +++ b/src/DPD-BASIC/pair_dpd.cpp @@ -18,17 +18,17 @@ #include "pair_dpd.h" -#include #include "atom.h" #include "comm.h" -#include "update.h" -#include "force.h" -#include "neighbor.h" -#include "neigh_list.h" -#include "random_mars.h" -#include "memory.h" #include "error.h" +#include "force.h" +#include "memory.h" +#include "neigh_list.h" +#include "neighbor.h" +#include "random_mars.h" +#include "update.h" +#include using namespace LAMMPS_NS; @@ -257,10 +257,10 @@ void PairDPD::init_style() // if newton off, forces between atoms ij will be double computed // using different random numbers - if (force->newton_pair == 0 && comm->me == 0) error->warning(FLERR, - "Pair dpd needs newton pair on for momentum conservation"); + if (force->newton_pair == 0 && comm->me == 0) + error->warning(FLERR, "Pair dpd needs newton pair on for momentum conservation"); - neighbor->request(this,instance_me); + neighbor->add_request(this); } /* ---------------------------------------------------------------------- diff --git a/src/DPD-BASIC/pair_dpd_ext.cpp b/src/DPD-BASIC/pair_dpd_ext.cpp index 82e0cde0d6..5e45790913 100644 --- a/src/DPD-BASIC/pair_dpd_ext.cpp +++ b/src/DPD-BASIC/pair_dpd_ext.cpp @@ -326,10 +326,10 @@ void PairDPDExt::init_style() // if newton off, forces between atoms ij will be double computed // using different random numbers - if (force->newton_pair == 0 && comm->me == 0) error->warning(FLERR, - "Pair dpd needs newton pair on for momentum conservation"); + if (force->newton_pair == 0 && comm->me == 0) + error->warning(FLERR, "Pair dpd needs newton pair on for momentum conservation"); - neighbor->request(this,instance_me); + neighbor->add_request(this); } /* ---------------------------------------------------------------------- diff --git a/src/DPD-BASIC/pair_dpd_tstat.cpp b/src/DPD-BASIC/pair_dpd_tstat.cpp index 71c9151856..502241db73 100644 --- a/src/DPD-BASIC/pair_dpd_tstat.cpp +++ b/src/DPD-BASIC/pair_dpd_tstat.cpp @@ -14,15 +14,15 @@ #include "pair_dpd_tstat.h" -#include #include "atom.h" -#include "update.h" +#include "comm.h" +#include "error.h" #include "force.h" #include "neigh_list.h" -#include "comm.h" #include "random_mars.h" -#include "error.h" +#include "update.h" +#include using namespace LAMMPS_NS; diff --git a/src/DPD-MESO/pair_edpd.cpp b/src/DPD-MESO/pair_edpd.cpp index 836ea27e8b..ddbbd05085 100644 --- a/src/DPD-MESO/pair_edpd.cpp +++ b/src/DPD-MESO/pair_edpd.cpp @@ -375,10 +375,10 @@ void PairEDPD::init_style() // if newton off, forces between atoms ij will be double computed // using different random numbers - if (force->newton_pair == 0 && comm->me == 0) error->warning(FLERR, - "Pair tdpd needs newton pair on for momentum conservation"); + if (force->newton_pair == 0 && comm->me == 0) + error->warning(FLERR, "Pair tdpd needs newton pair on for momentum conservation"); - neighbor->request(this,instance_me); + neighbor->add_request(this); } /* ---------------------------------------------------------------------- diff --git a/src/DPD-MESO/pair_mdpd.cpp b/src/DPD-MESO/pair_mdpd.cpp index 1b490aa190..053d322f00 100644 --- a/src/DPD-MESO/pair_mdpd.cpp +++ b/src/DPD-MESO/pair_mdpd.cpp @@ -287,10 +287,10 @@ void PairMDPD::init_style() // if newton off, forces between atoms ij will be double computed // using different random numbers - if (force->newton_pair == 0 && comm->me == 0) error->warning(FLERR, - "Pair mdpd needs newton pair on for momentum conservation"); + if (force->newton_pair == 0 && comm->me == 0) + error->warning(FLERR, "Pair mdpd needs newton pair on for momentum conservation"); - neighbor->request(this,instance_me); + neighbor->add_request(this); } /* ---------------------------------------------------------------------- diff --git a/src/DPD-MESO/pair_tdpd.cpp b/src/DPD-MESO/pair_tdpd.cpp index 70168c0e2a..35f440426e 100644 --- a/src/DPD-MESO/pair_tdpd.cpp +++ b/src/DPD-MESO/pair_tdpd.cpp @@ -328,10 +328,9 @@ void PairTDPD::init_style() // using different random numbers if (force->newton_pair == 0 && comm->me == 0) - error->warning(FLERR,"Pair tdpd needs newton pair on " - "for momentum conservation"); + error->warning(FLERR,"Pair tdpd needs newton pair on for momentum conservation"); - neighbor->request(this,instance_me); + neighbor->add_request(this); } /* ----------------------------------------------------------------------