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

convert more neighbor list requests to the new API

Browse Source
This commit is contained in:
Axel Kohlmeyer
2022-03-08 11:38:41 -05:00
parent deb694745b
commit dfbeedfded
31 changed files with 1056 additions and 1146 deletions
Download Patch File Download Diff File Expand all files Collapse all files

464
src/CG-SDK/pair_lj_sdk_coul_long.cpp
View File

@ -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 <cmath>
#include <cstring>
#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 <cmath>
#include <cstring>
#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 <int EVFLAG, int EFLAG, int NEWTON_PAIR>
void PairLJSDKCoulLong::eval()
template <int EVFLAG, int EFLAG, int NEWTON_PAIR> 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<<ncoultablebits;
bytes += (double)8 * ntable*sizeof(double);
int ntable = 1 << ncoultablebits;
bytes += (double) 8 * ntable * sizeof(double);
}
return bytes;

405
src/CLASS2/pair_lj_class2.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -19,7 +18,6 @@
#include "math_const.h"
#include "memory.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "neighbor.h"
#include "respa.h"
#include "update.h"
@ -65,13 +63,13 @@ PairLJClass2::~PairLJClass2()
void PairLJClass2::compute(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;
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;
int *ilist, *jlist, *numneigh, **firstneigh;
evdwl = 0.0;
ev_init(eflag,vflag);
ev_init(eflag, vflag);
double **x = atom->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;
}

7
src/CLASS2/pair_lj_class2_coul_cut.cpp
View File

@ -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 <cmath>
#include <cstring>
#include "atom.h"
#include "comm.h"
#include "force.h"
@ -25,6 +22,8 @@
#include "memory.h"
#include "error.h"
#include <cmath>
#include <cstring>
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);
}
/* ----------------------------------------------------------------------

652
src/CLASS2/pair_lj_class2_coul_long.cpp
View File

File diff suppressed because it is too large Load Diff

427
src/COLLOID/pair_brownian.cpp
View File

@ -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 <cmath>
#include <cstring>
#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 <cmath>
#include <cstring>
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];
}

26
src/COLLOID/pair_brownian_poly.cpp
View File

@ -19,24 +19,24 @@
#include "pair_brownian_poly.h"
#include <cmath>
#include <cstring>
#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 <cmath>
#include <cstring>
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

2
src/COLLOID/pair_colloid.cpp
View File

@ -18,7 +18,6 @@
#include "pair_colloid.h"
#include <cmath>
#include "atom.h"
#include "comm.h"
#include "force.h"
@ -27,6 +26,7 @@
#include "memory.h"
#include "error.h"
#include <cmath>
using namespace LAMMPS_NS;
using namespace MathSpecial;

18
src/COLLOID/pair_lubricate.cpp
View File

@ -19,24 +19,24 @@
#include "pair_lubricate.h"
#include <cmath>
#include <cstring>
#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 <cmath>
#include <cstring>
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

22
src/COLLOID/pair_lubricateU.cpp
View File

@ -18,24 +18,24 @@
#include "pair_lubricateU.h"
#include <cmath>
#include <cstring>
#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 <cmath>
#include <cstring>
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

22
src/COLLOID/pair_lubricateU_poly.cpp
View File

@ -20,23 +20,23 @@
#include "pair_lubricateU_poly.h"
#include <cmath>
#include <cstring>
#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 <cmath>
#include <cstring>
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);
}

22
src/COLLOID/pair_lubricate_poly.cpp
View File

@ -20,23 +20,23 @@
#include "pair_lubricate_poly.h"
#include <cmath>
#include <cstring>
#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 <cmath>
#include <cstring>
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

5
src/COLLOID/pair_yukawa_colloid.cpp
View File

@ -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");
}
/* ----------------------------------------------------------------------

8
src/DIELECTRIC/fix_polarize_functional.cpp
View File

@ -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;

5
src/DIELECTRIC/pair_coul_cut_dielectric.cpp
View File

@ -24,7 +24,6 @@
#include "math_const.h"
#include "memory.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "neighbor.h"
#include <cmath>
@ -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);
}
/* ---------------------------------------------------------------------- */

5
src/DIELECTRIC/pair_coul_long_dielectric.cpp
View File

@ -25,7 +25,6 @@
#include "math_const.h"
#include "memory.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "neighbor.h"
#include <cmath>
@ -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;

5
src/DIELECTRIC/pair_lj_cut_coul_cut_dielectric.cpp
View File

@ -24,7 +24,6 @@
#include "math_const.h"
#include "memory.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "neighbor.h"
#include <cmath>
@ -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);
}
/* ---------------------------------------------------------------------- */

5
src/DIELECTRIC/pair_lj_cut_coul_debye_dielectric.cpp
View File

@ -24,7 +24,6 @@
#include "math_const.h"
#include "memory.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "neighbor.h"
#include <cmath>
@ -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);
}
/* ---------------------------------------------------------------------- */

5
src/DIELECTRIC/pair_lj_cut_coul_long_dielectric.cpp
View File

@ -25,7 +25,6 @@
#include "math_const.h"
#include "memory.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "neighbor.h"
#include <cmath>
@ -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;

5
src/DIELECTRIC/pair_lj_cut_coul_msm_dielectric.cpp
View File

@ -25,7 +25,6 @@
#include "math_const.h"
#include "memory.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "neighbor.h"
#include <cmath>
@ -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;

5
src/DIELECTRIC/pair_lj_long_coul_long_dielectric.cpp
View File

@ -25,7 +25,6 @@
#include "math_extra.h"
#include "memory.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "neighbor.h"
#include <cmath>
@ -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);
}
/* ----------------------------------------------------------------------

12
src/DIPOLE/pair_lj_cut_dipole_cut.cpp
View File

@ -14,17 +14,17 @@
#include "pair_lj_cut_dipole_cut.h"
#include <cmath>
#include <cstring>
#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 <cmath>
#include <cstring>
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);
}
/* ----------------------------------------------------------------------

12
src/DIPOLE/pair_lj_cut_dipole_long.cpp
View File

@ -14,19 +14,19 @@
#include "pair_lj_cut_dipole_long.h"
#include <cmath>
#include <cstring>
#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 <cmath>
#include <cstring>
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);
}
/* ----------------------------------------------------------------------

2
src/DIPOLE/pair_lj_long_dipole_long.cpp
View File

@ -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;

8
src/DIPOLE/pair_lj_sf_dipole_sf.cpp
View File

@ -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);
}
/* ----------------------------------------------------------------------

2
src/DIPOLE/pair_lj_sf_dipole_sf.h
View File

@ -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;

20
src/DPD-BASIC/pair_dpd.cpp
View File

@ -18,17 +18,17 @@
#include "pair_dpd.h"
#include <cmath>
#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 <cmath>
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);
}
/* ----------------------------------------------------------------------

6
src/DPD-BASIC/pair_dpd_ext.cpp
View File

@ -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);
}
/* ----------------------------------------------------------------------

8
src/DPD-BASIC/pair_dpd_tstat.cpp
View File

@ -14,15 +14,15 @@
#include "pair_dpd_tstat.h"
#include <cmath>
#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 <cmath>
using namespace LAMMPS_NS;

6
src/DPD-MESO/pair_edpd.cpp
View File

@ -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);
}
/* ----------------------------------------------------------------------

6
src/DPD-MESO/pair_mdpd.cpp
View File

@ -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);
}
/* ----------------------------------------------------------------------

5
src/DPD-MESO/pair_tdpd.cpp
View File

@ -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);
}
/* ----------------------------------------------------------------------