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

apply clang-format

Browse Source
This commit is contained in:
Axel Kohlmeyer
2022-02-16 06:14:19 -05:00
parent 32cde04000
commit 08968cbdbf
7 changed files with 498 additions and 484 deletions
Download Patch File Download Diff File Expand all files Collapse all files

10
doc/src/pair_ilp_tmd.rst
View File

@ -30,7 +30,7 @@ Description
"""""""""""
The *ilp/tmd* style computes the registry-dependent interlayer
potential (ILP) potential for itransition metal dichalcogenide (TMD)
potential (ILP) potential for itransition metal dichalcogenide (TMD)
as described in :ref:`(Ouyang3) <Ouyang3>`.
.. math::
@ -57,11 +57,11 @@ calculating the normals.
.. note::
Since each MX2 (M = Mo, W and X = S, Se Te) layer contains two
Since each MX2 (M = Mo, W and X = S, Se Te) layer contains two
sublayers of X atoms and one sublayer of M atoms, the definition of the
normal vectors used for graphene and h-BN is no longer valid for TMDs.
In :ref:`(Ouyang3) <Ouyang3>`, a new definition is proposed, where for
each atom `i`, its six nearest neighboring atoms belonging to the same
In :ref:`(Ouyang3) <Ouyang3>`, a new definition is proposed, where for
each atom `i`, its six nearest neighboring atoms belonging to the same
sublayer are chosen to define the normal vector `{\bf n}_i`.
The parameter file (e.g. TMD.ILP), is intended for use with *metal*
@ -75,7 +75,7 @@ list for calculating the normals for each atom pair.
The parameters presented in the parameter file (e.g. BNCH.ILP),
are fitted with taper function by setting the cutoff equal to 16.0
Angstrom. Using different cutoff or taper function should be careful.
These parameters provide a good description in both short- and long-range
These parameters provide a good description in both short- and long-range
interaction regimes. This feature is essential for simulations in high pressure
regime (i.e., the interlayer distance is smaller than the equilibrium
distance). The benchmark tests and comparison of these parameters can

31
src/INTERLAYER/pair_ilp_graphene_hbn.cpp
View File

@ -58,8 +58,7 @@ static const char cite_ilp[] =
/* ---------------------------------------------------------------------- */
PairILPGrapheneHBN::PairILPGrapheneHBN(LAMMPS *lmp)
: Pair(lmp), variant(ILP_GrhBN)
PairILPGrapheneHBN::PairILPGrapheneHBN(LAMMPS *lmp) : Pair(lmp), variant(ILP_GrhBN)
{
restartinfo = 0;
one_coeff = 1;
@ -212,20 +211,20 @@ void PairILPGrapheneHBN::read_file(char *filename)
if (comm->me == 0) {
std::string potential_name;
switch (variant) {
case ILP_GrhBN:
potential_name = "ilp/graphene/hbn";
break;
case ILP_GrhBN:
potential_name = "ilp/graphene/hbn";
break;
case ILP_TMD:
potential_name = "ilp/tmd";
break;
case ILP_TMD:
potential_name = "ilp/tmd";
break;
case SAIP_METAL:
potential_name = "saip/metal";
break;
case SAIP_METAL:
potential_name = "saip/metal";
break;
default:
error->one(FLERR,"Unknown ILP style variant {}",variant);
default:
error->one(FLERR, "Unknown ILP style variant {}", variant);
}
PotentialFileReader reader(lmp, filename, potential_name, unit_convert_flag);
char *line;
@ -342,10 +341,8 @@ void PairILPGrapheneHBN::read_file(char *filename)
void PairILPGrapheneHBN::init_style()
{
if (force->newton_pair == 0)
error->all(FLERR, "Pair style ilp/* requires newton pair on");
if (!atom->molecule_flag)
error->all(FLERR, "Pair style ilp/* requires atom attribute molecule");
if (force->newton_pair == 0) error->all(FLERR, "Pair style ilp/* requires newton pair on");
if (!atom->molecule_flag) error->all(FLERR, "Pair style ilp/* requires atom attribute molecule");
// need a full neighbor list, including neighbors of ghosts

4
src/INTERLAYER/pair_ilp_graphene_hbn.h
View File

@ -68,8 +68,8 @@ class PairILPGrapheneHBN : public Pair {
double ***dnormdri;
double ****dnormal;
// adds for ilp/tmd
int Nnei; // max # of nearest neighbors for one atom
// adds for ilp/tmd
int Nnei; // max # of nearest neighbors for one atom
double **dnn;
double **vect;
double **pvet;

779
src/INTERLAYER/pair_ilp_tmd.cpp
View File

File diff suppressed because it is too large Load Diff

1
src/INTERLAYER/pair_ilp_tmd.h
View File

@ -43,7 +43,6 @@ class PairILPTMD : public PairILPGrapheneHBN {
return k % range;
}
/**************************************************************/
};
} // namespace LAMMPS_NS

12
src/INTERLAYER/pair_kolmogorov_crespi_z.cpp
View File

@ -28,8 +28,8 @@
#include "error.h"
#include "force.h"
#include "memory.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neighbor.h"
#include "potential_file_reader.h"
#include <cmath>
@ -160,15 +160,15 @@ void PairKolmogorovCrespiZ::compute(int eflag, int vflag)
if (evflag) {
ev_tally(i, j, nlocal, newton_pair, evdwl, 0.0, fpair, delx, dely, delz);
if (vflag_either) {
double fi[3],fj[3];
double fi[3], fj[3];
fi[0] = delx * fpair1;
fi[1] = dely * fpair1;
fi[2] = 0;
fj[0] = -delx * fpair1;
fj[1] = -dely * fpair1;
fj[2] = 0;
v_tally2_newton(i,fi,x[i]);
v_tally2_newton(j,fj,x[j]);
v_tally2_newton(i, fi, x[i]);
v_tally2_newton(j, fj, x[j]);
}
}
}
@ -246,9 +246,9 @@ void PairKolmogorovCrespiZ::coeff(int narg, char **arg)
void PairKolmogorovCrespiZ::init_style()
{
if (force->newton_pair == 0)
error->all(FLERR,"Pair style kolmogorov/crespi/z requires newton pair on");
error->all(FLERR, "Pair style kolmogorov/crespi/z requires newton pair on");
neighbor->request(this,instance_me);
neighbor->request(this, instance_me);
}
/* ----------------------------------------------------------------------

145
src/INTERLAYER/pair_saip_metal.cpp
View File

@ -52,7 +52,8 @@ static const char cite_saip[] =
/* ---------------------------------------------------------------------- */
PairSAIPMETAL::PairSAIPMETAL(LAMMPS *lmp) : PairILPGrapheneHBN(lmp) {
PairSAIPMETAL::PairSAIPMETAL(LAMMPS *lmp) : PairILPGrapheneHBN(lmp)
{
variant = SAIP_METAL;
if (lmp->citeme) lmp->citeme->add(cite_saip);
}
@ -71,19 +72,18 @@ void PairSAIPMETAL::settings(int narg, char **arg)
if (narg == 2) tap_flag = utils::numeric(FLERR, arg[1], false, lmp);
}
/* ----------------------------------------------------------------------
Repulsive forces and energy
------------------------------------------------------------------------- */
void PairSAIPMETAL::calc_FRep(int eflag, int /* vflag */)
{
int i,j,ii,jj,inum,jnum,itype,jtype,k,kk;
double prodnorm1,fkcx,fkcy,fkcz,filp;
double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair,fpair1;
double rsq,r,Rcut,rhosq1,exp0,exp1,Tap,dTap,Vilp;
double frho1,Erep,fsum,rdsq1;
int *ilist,*jlist,*numneigh,**firstneigh;
int i, j, ii, jj, inum, jnum, itype, jtype, k, kk;
double prodnorm1, fkcx, fkcy, fkcz, filp;
double xtmp, ytmp, ztmp, delx, dely, delz, evdwl, fpair, fpair1;
double rsq, r, Rcut, rhosq1, exp0, exp1, Tap, dTap, Vilp;
double frho1, Erep, fsum, rdsq1;
int *ilist, *jlist, *numneigh, **firstneigh;
int *ILP_neighs_i;
evdwl = 0.0;
@ -123,35 +123,38 @@ void PairSAIPMETAL::calc_FRep(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;
// only include the interaction between different layers
if (rsq < cutsq[itype][jtype] && atom->molecule[i] != atom->molecule[j]) {
int iparam_ij = elem2param[map[itype]][map[jtype]];
Param& p = params[iparam_ij];
Param &p = params[iparam_ij];
r = sqrt(rsq);
// turn on/off taper function
if (tap_flag) {
Rcut = sqrt(cutsq[itype][jtype]);
Tap = calc_Tap(r,Rcut);
dTap = calc_dTap(r,Rcut);
} else {Tap = 1.0; dTap = 0.0;}
Tap = calc_Tap(r, Rcut);
dTap = calc_dTap(r, Rcut);
} else {
Tap = 1.0;
dTap = 0.0;
}
// for atoms in bulk materials
if (strcmp(elements[map[itype]],"C") != 0 && strcmp(elements[map[itype]],"H") != 0 &&
strcmp(elements[map[itype]],"B") != 0 && strcmp(elements[map[itype]],"N") != 0) {
if (strcmp(elements[map[itype]], "C") != 0 && strcmp(elements[map[itype]], "H") != 0 &&
strcmp(elements[map[itype]], "B") != 0 && strcmp(elements[map[itype]], "N") != 0) {
// Set ni along rij
exp0 = exp(-p.lambda*(r-p.z0));
frho1 = 1.0*p.C;
Erep = 0.5*p.epsilon + frho1;
Vilp = exp0*Erep;
exp0 = exp(-p.lambda * (r - p.z0));
frho1 = 1.0 * p.C;
Erep = 0.5 * p.epsilon + frho1;
Vilp = exp0 * Erep;
// derivatives
fpair = p.lambda*exp0/r*Erep;
filp = fpair*Tap - Vilp*dTap/r;
fkcx = delx*filp;
fkcy = dely*filp;
fkcz = delz*filp;
fpair = p.lambda * exp0 / r * Erep;
filp = fpair * Tap - Vilp * dTap / r;
fkcx = delx * filp;
fkcy = dely * filp;
fkcz = delz * filp;
f[i][0] += fkcx;
f[i][1] += fkcy;
@ -159,46 +162,48 @@ void PairSAIPMETAL::calc_FRep(int eflag, int /* vflag */)
f[j][0] -= fkcx;
f[j][1] -= fkcy;
f[j][2] -= fkcz;
if (eflag) pvector[1] += evdwl = Tap*Vilp;
if (evflag) ev_tally(i,j,nlocal,newton_pair,evdwl,0.0,filp,delx,dely,delz);
}
else { // for atoms in 2D materials
if (eflag) pvector[1] += evdwl = Tap * Vilp;
if (evflag) ev_tally(i, j, nlocal, newton_pair, evdwl, 0.0, filp, delx, dely, delz);
} else { // for atoms in 2D materials
// Calculate the transverse distance
prodnorm1 = normal[i][0]*delx + normal[i][1]*dely + normal[i][2]*delz;
rhosq1 = rsq - prodnorm1*prodnorm1; // rho_ij
rdsq1 = rhosq1*p.delta2inv; // (rho_ij/delta)^2
prodnorm1 = normal[i][0] * delx + normal[i][1] * dely + normal[i][2] * delz;
rhosq1 = rsq - prodnorm1 * prodnorm1; // rho_ij
rdsq1 = rhosq1 * p.delta2inv; // (rho_ij/delta)^2
// store exponents
exp0 = exp(-p.lambda*(r-p.z0));
exp0 = exp(-p.lambda * (r - p.z0));
exp1 = exp(-rdsq1);
frho1 = exp1*p.C;
Erep = 0.5*p.epsilon + frho1;
Vilp = exp0*Erep;
frho1 = exp1 * p.C;
Erep = 0.5 * p.epsilon + frho1;
Vilp = exp0 * Erep;
// derivatives
fpair = p.lambda*exp0/r*Erep;
fpair1 = 2.0*exp0*frho1*p.delta2inv;
fpair = p.lambda * exp0 / r * Erep;
fpair1 = 2.0 * exp0 * frho1 * p.delta2inv;
fsum = fpair + fpair1;
// derivatives of the product of rij and ni, the result is a vector
dprodnorm1[0] = dnormdri[0][0][i]*delx + dnormdri[1][0][i]*dely + dnormdri[2][0][i]*delz;
dprodnorm1[1] = dnormdri[0][1][i]*delx + dnormdri[1][1][i]*dely + dnormdri[2][1][i]*delz;
dprodnorm1[2] = dnormdri[0][2][i]*delx + dnormdri[1][2][i]*dely + dnormdri[2][2][i]*delz;
fp1[0] = prodnorm1*normal[i][0]*fpair1;
fp1[1] = prodnorm1*normal[i][1]*fpair1;
fp1[2] = prodnorm1*normal[i][2]*fpair1;
fprod1[0] = prodnorm1*dprodnorm1[0]*fpair1;
fprod1[1] = prodnorm1*dprodnorm1[1]*fpair1;
fprod1[2] = prodnorm1*dprodnorm1[2]*fpair1;
dprodnorm1[0] =
dnormdri[0][0][i] * delx + dnormdri[1][0][i] * dely + dnormdri[2][0][i] * delz;
dprodnorm1[1] =
dnormdri[0][1][i] * delx + dnormdri[1][1][i] * dely + dnormdri[2][1][i] * delz;
dprodnorm1[2] =
dnormdri[0][2][i] * delx + dnormdri[1][2][i] * dely + dnormdri[2][2][i] * delz;
fp1[0] = prodnorm1 * normal[i][0] * fpair1;
fp1[1] = prodnorm1 * normal[i][1] * fpair1;
fp1[2] = prodnorm1 * normal[i][2] * fpair1;
fprod1[0] = prodnorm1 * dprodnorm1[0] * fpair1;
fprod1[1] = prodnorm1 * dprodnorm1[1] * fpair1;
fprod1[2] = prodnorm1 * dprodnorm1[2] * fpair1;
fkcx = (delx*fsum - fp1[0])*Tap - Vilp*dTap*delx/r;
fkcy = (dely*fsum - fp1[1])*Tap - Vilp*dTap*dely/r;
fkcz = (delz*fsum - fp1[2])*Tap - Vilp*dTap*delz/r;
fkcx = (delx * fsum - fp1[0]) * Tap - Vilp * dTap * delx / r;
fkcy = (dely * fsum - fp1[1]) * Tap - Vilp * dTap * dely / r;
fkcz = (delz * fsum - fp1[2]) * Tap - Vilp * dTap * delz / r;
f[i][0] += fkcx - fprod1[0]*Tap;
f[i][1] += fkcy - fprod1[1]*Tap;
f[i][2] += fkcz - fprod1[2]*Tap;
f[i][0] += fkcx - fprod1[0] * Tap;
f[i][1] += fkcy - fprod1[1] * Tap;
f[i][2] += fkcz - fprod1[2] * Tap;
f[j][0] -= fkcx;
f[j][1] -= fkcy;
f[j][2] -= fkcz;
@ -209,25 +214,31 @@ void PairSAIPMETAL::calc_FRep(int eflag, int /* vflag */)
k = ILP_neighs_i[kk];
if (k == i) continue;
// derivatives of the product of rij and ni respect to rk, k=0,1,2, where atom k is the neighbors of atom i
dprodnorm1[0] = dnormal[0][0][kk][i]*delx + dnormal[1][0][kk][i]*dely + dnormal[2][0][kk][i]*delz;
dprodnorm1[1] = dnormal[0][1][kk][i]*delx + dnormal[1][1][kk][i]*dely + dnormal[2][1][kk][i]*delz;
dprodnorm1[2] = dnormal[0][2][kk][i]*delx + dnormal[1][2][kk][i]*dely + dnormal[2][2][kk][i]*delz;
fk[0] = (-prodnorm1*dprodnorm1[0]*fpair1)*Tap;
fk[1] = (-prodnorm1*dprodnorm1[1]*fpair1)*Tap;
fk[2] = (-prodnorm1*dprodnorm1[2]*fpair1)*Tap;
dprodnorm1[0] = dnormal[0][0][kk][i] * delx + dnormal[1][0][kk][i] * dely +
dnormal[2][0][kk][i] * delz;
dprodnorm1[1] = dnormal[0][1][kk][i] * delx + dnormal[1][1][kk][i] * dely +
dnormal[2][1][kk][i] * delz;
dprodnorm1[2] = dnormal[0][2][kk][i] * delx + dnormal[1][2][kk][i] * dely +
dnormal[2][2][kk][i] * delz;
fk[0] = (-prodnorm1 * dprodnorm1[0] * fpair1) * Tap;
fk[1] = (-prodnorm1 * dprodnorm1[1] * fpair1) * Tap;
fk[2] = (-prodnorm1 * dprodnorm1[2] * fpair1) * Tap;
f[k][0] += fk[0];
f[k][1] += fk[1];
f[k][2] += fk[2];
delki[0] = x[k][0] - x[i][0];
delki[1] = x[k][1] - x[i][1];
delki[2] = x[k][2] - x[i][2];
if (evflag) ev_tally_xyz(k,j,nlocal,newton_pair,0.0,0.0,fk[0],fk[1],fk[2],delki[0],delki[1],delki[2]);
if (evflag)
ev_tally_xyz(k, j, nlocal, newton_pair, 0.0, 0.0, fk[0], fk[1], fk[2], delki[0],
delki[1], delki[2]);
}
if (eflag) pvector[1] += evdwl = Tap*Vilp;
if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair,evdwl,0.0,fkcx,fkcy,fkcz,delx,dely,delz);
} // end of speration atoms for bulk materials
} // end of rsq < cutoff
} // loop over jj
} // loop over ii
if (eflag) pvector[1] += evdwl = Tap * Vilp;
if (evflag)
ev_tally_xyz(i, j, nlocal, newton_pair, evdwl, 0.0, fkcx, fkcy, fkcz, delx, dely, delz);
} // end of speration atoms for bulk materials
} // end of rsq < cutoff
} // loop over jj
} // loop over ii
}