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

Propagated changes in PairTersoff to tersoff/mod and tersoff/mod/c

Browse Source
This commit is contained in:
Aidan Thompson
2020-12-29 12:54:29 -07:00
parent 2e623ff0a5
commit 3e6bd3ef32
8 changed files with 35 additions and 62 deletions
Download Patch File Download Diff File Expand all files Collapse all files

47
src/MANYBODY/pair_tersoff.cpp
View File

@ -43,7 +43,6 @@ using namespace MathSpecial;
PairTersoff::PairTersoff(LAMMPS *lmp) : Pair(lmp)
{
shift_flag = 0;
single_enable = 0;
restartinfo = 0;
one_coeff = 1;
@ -91,7 +90,8 @@ void PairTersoff::compute(int eflag, int vflag)
int i,j,k,ii,jj,kk,inum,jnum;
int itype,jtype,ktype,iparam_ij,iparam_ijk;
tagint itag,jtag;
double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair, fforce;
double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
double fforce;
double rsq,rsq1,rsq2;
double delr1[3],delr2[3],fi[3],fj[3],fk[3];
double r1_hat[3],r2_hat[3];
@ -143,7 +143,9 @@ void PairTersoff::compute(int eflag, int vflag)
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
if(shift_flag == 1) {
// shift rsq and store correction for force
if(shift_flag) {
double rsqtmp = rsq + shift*shift + 2*sqrt(rsq)*shift;
forceshiftfac = sqrt(rsqtmp/rsq);
rsq = rsqtmp;
@ -174,7 +176,9 @@ void PairTersoff::compute(int eflag, int vflag)
repulsive(&params[iparam_ij],rsq,fpair,eflag,evdwl);
if(shift_flag == 1) fpair *= forceshiftfac;
// correct force for shift in rsq
if(shift_flag) fpair *= forceshiftfac;
fxtmp += delx*fpair;
fytmp += dely*fpair;
@ -201,14 +205,13 @@ void PairTersoff::compute(int eflag, int vflag)
delr1[2] = x[j][2] - ztmp;
rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2];
if(shift_flag == 1)
if(shift_flag)
rsq1 += shift*shift + 2*sqrt(rsq1)*shift;
if (rsq1 >= params[iparam_ij].cutsq) continue;
double r1inv = 1.0/sqrt(vec3_dot(delr1, delr1));
vec3_scale(r1inv, delr1, r1_hat);
// vec3_norm(delr1, r1_hat);
// accumulate bondorder zeta for each i-j interaction via loop over k
@ -226,14 +229,13 @@ void PairTersoff::compute(int eflag, int vflag)
delr2[2] = x[k][2] - ztmp;
rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2];
if(shift_flag == 1)
if(shift_flag)
rsq2 += shift*shift + 2*sqrt(rsq2)*shift;
if (rsq2 >= params[iparam_ijk].cutsq) continue;
double r2inv = 1.0/sqrt(vec3_dot(delr2, delr2));
vec3_scale(r2inv, delr2, r2_hat);
// vec3_norm(delr2, r2_hat);
zeta_ij += zeta(&params[iparam_ijk],rsq1,rsq2,r1_hat,r2_hat);
}
@ -267,7 +269,7 @@ void PairTersoff::compute(int eflag, int vflag)
delr2[2] = x[k][2] - ztmp;
rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2];
if(shift_flag == 1)
if(shift_flag)
rsq2 += shift*shift + 2*sqrt(rsq2)*shift;
if (rsq2 >= params[iparam_ijk].cutsq) continue;
@ -321,16 +323,23 @@ void PairTersoff::allocate()
void PairTersoff::settings(int narg, char **arg)
{
if (narg == 0) {
shift_flag = 0;
}
else if (narg == 2) {
if (strcmp(arg[0],"shift") == 0) {
// default values
shift_flag = 0;
// process optional keywords
int iarg = 0;
while (iarg < narg) {
if (strcmp(arg[iarg],"shift") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal option for pair_tersoff");
shift = utils::numeric(FLERR,arg[iarg+1],false,lmp);
shift_flag = 1;
shift = utils::numeric(FLERR,arg[1],false,lmp);
iarg += 2;
} else error->all(FLERR,"Illegal option for pair_tersoff");
}
else error->all(FLERR,"Illegal pair_style command");
}
/* ----------------------------------------------------------------------
@ -620,11 +629,10 @@ void PairTersoff::repulsive(Param *param, double rsq, double &fforce,
double PairTersoff::zeta(Param *param, double rsqij, double rsqik,
double *rij_hat, double *rik_hat)
{
double rij,rik,rijs,riks,costheta,arg,ex_delr;
double rij,rik,costheta,arg,ex_delr;
rij = sqrt(rsqij);
rik = sqrt(rsqik);
costheta = vec3_dot(rij_hat,rik_hat);
if (param->powermint == 3) arg = cube(param->lam3 * (rij-rik));
@ -672,6 +680,8 @@ void PairTersoff::attractive(Param *param, double prefactor,
rijinv = 1.0/rij;
rikinv = 1.0/rik;
// correct 1/r for shift in rsq
if(shift_flag == 1) {
rijinv *= sqrt((rsqij + shift*shift + 2*sqrt(rsqij)*shift)/rsqij);
rikinv *= sqrt((rsqik + shift*shift + 2*sqrt(rsqik)*shift)/rsqik);
@ -830,5 +840,4 @@ void PairTersoff::costheta_d(double *rij_hat, double rijinv,
vec3_scale(rikinv,drk,drk);
vec3_add(drj,drk,dri);
vec3_scale(-1.0,dri,dri);
// printf("costheta_d = %g %g %g\n",rik_hat[0], cos_theta, dri[0]);
}

5
src/MANYBODY/pair_tersoff.h
View File

@ -134,11 +134,6 @@ class PairTersoff : public Pair {
z[2] = k*x[2]+y[2];
}
inline void vec3_norm(const double x[3], double * const y) const {
double f = 1.0/sqrt(vec3_dot(x, x));
vec3_scale(f, x, y);
}
};
}

20
src/MANYBODY/pair_tersoff_mod.cpp
View File

@ -39,15 +39,6 @@ using namespace MathSpecial;
PairTersoffMOD::PairTersoffMOD(LAMMPS *lmp) : PairTersoff(lmp) {}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairTersoffMOD::settings(int narg, char **/*arg*/)
{
if (narg != 0) error->all(FLERR,"Illegal pair_style command");
shift_flag = 0;
}
/* ---------------------------------------------------------------------- */
void PairTersoffMOD::read_file(char *file)
@ -215,14 +206,13 @@ void PairTersoffMOD::setup_params()
/* ---------------------------------------------------------------------- */
double PairTersoffMOD::zeta(Param *param, double rsqij, double rsqik,
double *delrij, double *delrik)
double *rij_hat, double *rik_hat)
{
double rij,rik,costheta,arg,ex_delr;
rij = sqrt(rsqij);
rik = sqrt(rsqik);
costheta = (delrij[0]*delrik[0] + delrij[1]*delrik[1] +
delrij[2]*delrik[2]) / (rij*rik);
costheta = vec3_dot(rij_hat,rik_hat);
if (param->powermint == 3) arg = cube(param->lam3 * (rij-rik));
else arg = param->lam3 * (rij-rik);
@ -287,8 +277,8 @@ double PairTersoffMOD::ters_bij_d(double zeta, Param *param)
/* ---------------------------------------------------------------------- */
void PairTersoffMOD::ters_zetaterm_d(double prefactor,
double *rij_hat, double rij,
double *rik_hat, double rik,
double *rij_hat, double rij, double rijinv,
double *rik_hat, double rik, double rikinv,
double *dri, double *drj, double *drk,
Param *param)
{
@ -311,7 +301,7 @@ void PairTersoffMOD::ters_zetaterm_d(double prefactor,
cos_theta = vec3_dot(rij_hat,rik_hat);
gijk = ters_gijk_mod(cos_theta,param);
gijk_d = ters_gijk_d_mod(cos_theta,param);
costheta_d(rij_hat,rij,rik_hat,rik,dcosdri,dcosdrj,dcosdrk);
costheta_d(rij_hat,rijinv,rik_hat,rikinv,dcosdri,dcosdrj,dcosdrk);
// compute the derivative wrt Ri
// dri = -dfc*gijk*ex_delr*rik_hat;

4
src/MANYBODY/pair_tersoff_mod.h
View File

@ -33,7 +33,6 @@ class PairTersoffMOD : public PairTersoff {
static const int NPARAMS_PER_LINE = 20;
protected:
void settings(int, char **);
virtual void read_file(char *);
virtual void setup_params();
double zeta(Param *, double, double, double *, double *);
@ -42,7 +41,8 @@ class PairTersoffMOD : public PairTersoff {
double ters_fc_d(double, Param *);
double ters_bij(double, Param *);
double ters_bij_d(double, Param *);
void ters_zetaterm_d(double, double *, double, double *, double,
void ters_zetaterm_d(double, double *, double, double,
double *, double, double,
double *, double *, double *, Param *);
// inlined functions for efficiency

9
src/MANYBODY/pair_tersoff_mod_c.cpp
View File

@ -30,15 +30,6 @@ using namespace LAMMPS_NS;
#define DELTA 4
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairTersoffMODC::settings(int narg, char **/*arg*/)
{
if (narg != 0) error->all(FLERR,"Illegal pair_style command");
shift_flag = 0;
}
/* ---------------------------------------------------------------------- */
void PairTersoffMODC::read_file(char *file)

1
src/MANYBODY/pair_tersoff_mod_c.h
View File

@ -28,7 +28,6 @@ class PairTersoffMODC : public PairTersoffMOD {
public:
PairTersoffMODC(class LAMMPS *lmp) : PairTersoffMOD(lmp) {};
~PairTersoffMODC() {}
void settings(int, char **);
static const int NPARAMS_PER_LINE = 21;

10
src/MANYBODY/pair_tersoff_zbl.cpp
View File

@ -57,16 +57,6 @@ PairTersoffZBL::PairTersoffZBL(LAMMPS *lmp) : PairTersoff(lmp)
} else error->all(FLERR,"Pair tersoff/zbl requires metal or real units");
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairTersoffZBL::settings(int narg, char **/*arg*/)
{
if (narg != 0) error->all(FLERR,"Illegal pair_style command");
shift_flag = 0;
}
/* ---------------------------------------------------------------------- */
void PairTersoffZBL::read_file(char *file)

1
src/MANYBODY/pair_tersoff_zbl.h
View File

@ -28,7 +28,6 @@ class PairTersoffZBL : public PairTersoff {
public:
PairTersoffZBL(class LAMMPS *);
~PairTersoffZBL() {}
void settings(int, char **);
static const int NPARAMS_PER_LINE = 21;