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

reformat

Browse Source
This commit is contained in:
Axel Kohlmeyer
2021-07-29 14:44:15 -04:00
parent 5b6b6b8767
commit ebe5e6c024
1 changed files with 115 additions and 108 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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