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
2024-09-09 21:12:03 -04:00
parent 68e134773c
commit 29bd3fc11b
2 changed files with 103 additions and 79 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
src/.gitignore vendored
View File

@ -437,6 +437,8 @@
/pair_coul_cut_global.h
/pair_coul_streitz.cpp
/pair_coul_streitz.h
/pair_coul_ctip.cpp
/pair_coul_ctip.h
/pair_gauss.cpp
/pair_gauss.h
/pair_lj96_cut.cpp
@ -976,6 +978,8 @@
/fix_pour.h
/fix_qeq_comb.cpp
/fix_qeq_comb.h
/fix_qeq_ctip.cpp
/fix_qeq_ctip.h
/fix_qeq_fire.cpp
/fix_qeq_fire.h
/fix_qeq_reaxff.cpp

178
src/EXTRA-PAIR/pair_coul_ctip.cpp
View File

@ -20,6 +20,7 @@
#include "atom.h"
#include "comm.h"
#include "error.h"
#include "ewald_const.h"
#include "force.h"
#include "math_const.h"
#include "memory.h"
@ -32,21 +33,15 @@
using namespace LAMMPS_NS;
using namespace MathConst;
using namespace EwaldConst;
#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 DELTA 4
static constexpr int DELTA = 4;
/* ---------------------------------------------------------------------- */
PairCoulCTIP::PairCoulCTIP(LAMMPS *lmp) : Pair(lmp)
{
params= nullptr;
params = nullptr;
}
/* ---------------------------------------------------------------------- */
@ -87,8 +82,10 @@ void PairCoulCTIP::compute(int eflag, int vflag)
double cut_coulcu, cut_coul4, alphacu, erfcd_cut, t_cut, erfcc_cut;
int elt1, elt2;
double shield[nelements][nelements], shieldcu[nelements][nelements], reffc[nelements][nelements], reffcsq[nelements][nelements], reffc4[nelements][nelements], reffc7[nelements][nelements];
double s2d_shift[nelements][nelements], f_shift[nelements][nelements], e_shift[nelements][nelements], self_factor[nelements][nelements];
double shield[nelements][nelements], shieldcu[nelements][nelements], reffc[nelements][nelements],
reffcsq[nelements][nelements], reffc4[nelements][nelements], reffc7[nelements][nelements];
double s2d_shift[nelements][nelements], f_shift[nelements][nelements],
e_shift[nelements][nelements], self_factor[nelements][nelements];
ecoul = 0.0;
selfion = 0.0;
@ -123,9 +120,14 @@ void PairCoulCTIP::compute(int eflag, int vflag)
reffcsq[elt1][elt2] = reffc[elt1][elt2] * reffc[elt1][elt2];
reffc4[elt1][elt2] = reffcsq[elt1][elt2] * reffcsq[elt1][elt2];
reffc7[elt1][elt2] = reffc4[elt1][elt2] * reffcsq[elt1][elt2] * reffc[elt1][elt2];
s2d_shift[elt1][elt2] = 2.0 * erfcc_cut / cut_coulcu + 4.0 * alpha / MY_PIS * erfcd_cut / cut_coulsq + 4.0 * alphacu / MY_PIS * erfcd_cut - 2 / cut_coulcu + 4 * cut_coul4 / reffc7[elt1][elt2] - 2 * cut_coul / reffc4[elt1][elt2];
f_shift[elt1][elt2] = erfcc_cut / cut_coulsq + 2.0 * alpha / MY_PIS * erfcd_cut / cut_coul - 1 / cut_coulsq + cut_coulsq / reffc4[elt1][elt2];
e_shift[elt1][elt2] = 2.0 * erfcc_cut / cut_coul + 2.0 * alpha / MY_PIS * erfcd_cut - 2 / cut_coul + 1 / reffc[elt1][elt2] + cut_coulcu / reffc4[elt1][elt2] + s2d_shift[elt1][elt2] * cut_coulsq * 0.5;
s2d_shift[elt1][elt2] = 2.0 * erfcc_cut / cut_coulcu +
4.0 * alpha / MY_PIS * erfcd_cut / cut_coulsq + 4.0 * alphacu / MY_PIS * erfcd_cut -
2 / cut_coulcu + 4 * cut_coul4 / reffc7[elt1][elt2] - 2 * cut_coul / reffc4[elt1][elt2];
f_shift[elt1][elt2] = erfcc_cut / cut_coulsq + 2.0 * alpha / MY_PIS * erfcd_cut / cut_coul -
1 / cut_coulsq + cut_coulsq / reffc4[elt1][elt2];
e_shift[elt1][elt2] = 2.0 * erfcc_cut / cut_coul + 2.0 * alpha / MY_PIS * erfcd_cut -
2 / cut_coul + 1 / reffc[elt1][elt2] + cut_coulcu / reffc4[elt1][elt2] +
s2d_shift[elt1][elt2] * cut_coulsq * 0.5;
self_factor[elt1][elt2] = -(e_shift[elt1][elt2] * 0.5 + alpha / MY_PIS) * qqrd2e;
}
}
@ -143,9 +145,9 @@ void PairCoulCTIP::compute(int eflag, int vflag)
jnum = numneigh[i];
iparam_i = elem1param[itype];
selfion = self(&params[iparam_i],qtmp);
selfion = self(&params[iparam_i], qtmp);
if (evflag) ev_tally(i,i,nlocal,0,0.0,selfion,0.0,0.0,0.0,0.0);
if (evflag) ev_tally(i, i, nlocal, 0, 0.0, selfion, 0.0, 0.0, 0.0, 0.0);
if (eflag) {
double e_self = self_factor[iparam_i][iparam_i] * qtmp * qtmp;
@ -166,7 +168,7 @@ void PairCoulCTIP::compute(int eflag, int vflag)
jtype = map[type[j]];
jparam_j = elem1param[jtype];
r = sqrt(rsq);
reff = cbrt(rsq * r + 1/shieldcu[iparam_i][jparam_j]);
reff = cbrt(rsq * r + 1 / shieldcu[iparam_i][jparam_j]);
reffsq = reff * reff;
reff4 = reffsq * reffsq;
prefactor = qqrd2e * qtmp * q[j] / r;
@ -174,7 +176,10 @@ void PairCoulCTIP::compute(int eflag, int vflag)
t = 1.0 / (1.0 + EWALD_P * alpha * r);
erfcc = t * (A1 + t * (A2 + t * (A3 + t * (A4 + t * A5)))) * erfcd;
forcecoul = prefactor * (erfcc / r + 2.0 * alpha / MY_PIS * erfcd + rsq * r / reff4 - 1 / r - r * f_shift[iparam_i][jparam_j] + r * s2d_shift[iparam_i][jparam_j] * (r - cut_coul)) * r;
forcecoul = prefactor *
(erfcc / r + 2.0 * alpha / MY_PIS * erfcd + rsq * r / reff4 - 1 / r -
r * f_shift[iparam_i][jparam_j] + r * s2d_shift[iparam_i][jparam_j] * (r - cut_coul)) *
r;
if (factor_coul < 1.0) forcecoul -= (1.0 - factor_coul) * prefactor;
fpair = forcecoul / rsq;
@ -188,7 +193,10 @@ void PairCoulCTIP::compute(int eflag, int vflag)
}
if (eflag) {
ecoul = prefactor * (erfcc + r / reff - 1 - r * erfcc_cut / cut_coul - r / reffc[iparam_i][jparam_j] + r / cut_coul + r * f_shift[iparam_i][jparam_j] * (r - cut_coul) - r * s2d_shift[iparam_i][jparam_j] * (r - cut_coul) * (r - cut_coul)* 0.5);
ecoul = prefactor *
(erfcc + r / reff - 1 - r * erfcc_cut / cut_coul - r / reffc[iparam_i][jparam_j] +
r / cut_coul + r * f_shift[iparam_i][jparam_j] * (r - cut_coul) -
r * s2d_shift[iparam_i][jparam_j] * (r - cut_coul) * (r - cut_coul) * 0.5);
if (factor_coul < 1.0) ecoul -= (1.0 - factor_coul) * prefactor;
} else
ecoul = 0.0;
@ -205,17 +213,17 @@ void PairCoulCTIP::compute(int eflag, int vflag)
double PairCoulCTIP::self(Param *param, double qi)
{
double s1=param->chi, s2=param->eta, s3=param->qmin, s4=param->qmax, s5=param->omega;
double s1 = param->chi, s2 = param->eta, s3 = param->qmin, s4 = param->qmax, s5 = param->omega;
if ( qi < s3 ) {
return qi*((s1-2*s3*s5)+qi*(0.50*s2+s5))+s3*s3*s5;
} else if ( qi < s4 ) {
return qi*(s1+qi*(0.50*s2));
} else {
return qi*((s1-2*s4*s5)+qi*(0.50*s2+s5))+s4*s4*s5;
}
if (qi < s3) {
return qi * ((s1 - 2 * s3 * s5) + qi * (0.50 * s2 + s5)) + s3 * s3 * s5;
} else if (qi < s4) {
return qi * (s1 + qi * (0.50 * s2));
} else {
return qi * ((s1 - 2 * s4 * s5) + qi * (0.50 * s2 + s5)) + s4 * s4 * s5;
}
return 0.0;
return 0.0;
}
/* ----------------------------------------------------------------------
@ -233,16 +241,16 @@ void PairCoulCTIP::allocate()
memory->create(cutsq, n + 1, n + 1, "pair:cutsq");
memory->create(qeq_x,n+1,"pair:qeq_x");
memory->create(qeq_j,n+1,"pair:qeq_j");
memory->create(qeq_g,n+1,"pair:qeq_g");
memory->create(qeq_z,n+1,"pair:qeq_z");
memory->create(qeq_c,n+1,"pair:qeq_c");
memory->create(qeq_q1,n+1,"pair:qeq_q1");
memory->create(qeq_q2,n+1,"pair:qeq_q2");
memory->create(qeq_w,n+1,"pair:qeq_w");
memory->create(qeq_x, n + 1, "pair:qeq_x");
memory->create(qeq_j, n + 1, "pair:qeq_j");
memory->create(qeq_g, n + 1, "pair:qeq_g");
memory->create(qeq_z, n + 1, "pair:qeq_z");
memory->create(qeq_c, n + 1, "pair:qeq_c");
memory->create(qeq_q1, n + 1, "pair:qeq_q1");
memory->create(qeq_q2, n + 1, "pair:qeq_q2");
memory->create(qeq_w, n + 1, "pair:qeq_w");
map = new int[n+1];
map = new int[n + 1];
}
/* ----------------------------------------------------------------------
@ -265,13 +273,12 @@ void PairCoulCTIP::coeff(int narg, char **arg)
{
if (!allocated) allocate();
map_element2type(narg-3, arg+3);
map_element2type(narg - 3, arg + 3);
// read potential file and initialize potential parameters
read_file(arg[2]);
setup_params();
}
/* ----------------------------------------------------------------------
@ -285,7 +292,6 @@ void PairCoulCTIP::init_style()
neighbor->add_request(this);
cut_coulsq = cut_coul * cut_coul;
}
/* ----------------------------------------------------------------------
@ -309,7 +315,7 @@ void PairCoulCTIP::read_file(char *file)
// open file on proc 0
if (comm->me == 0) {
PotentialFileReader reader(lmp, file, "coul/ctip");
char * line;
char *line;
while ((line = reader.next_line(NPARAMS_PER_LINE))) {
try {
@ -327,16 +333,14 @@ void PairCoulCTIP::read_file(char *file)
if (nparams == maxparam) {
maxparam += DELTA;
params = (Param *) memory->srealloc(params,maxparam*sizeof(Param),
"pair:params");
params = (Param *) memory->srealloc(params, maxparam * sizeof(Param), "pair:params");
// make certain all addional allocated storage is initialized
// to avoid false positives when checking with valgrind
memset(params + nparams, 0, DELTA*sizeof(Param));
memset(params + nparams, 0, DELTA * sizeof(Param));
}
params[nparams].ielement = ielement;
params[nparams].chi = values.next_double();
params[nparams].eta = values.next_double();
@ -353,9 +357,10 @@ void PairCoulCTIP::read_file(char *file)
// parameter sanity check
if (params[nparams].eta < 0.0 || params[nparams].zeta < 0.0 ||
params[nparams].zcore < 0.0 || params[nparams].gamma < 0.0 || params[nparams].qmin > params[nparams].qmax || params[nparams].omega < 0.0 )
error->one(FLERR,"Illegal coul/ctip parameter");
if (params[nparams].eta < 0.0 || params[nparams].zeta < 0.0 || params[nparams].zcore < 0.0 ||
params[nparams].gamma < 0.0 || params[nparams].qmin > params[nparams].qmax ||
params[nparams].omega < 0.0)
error->one(FLERR, "Illegal coul/ctip parameter");
nparams++;
}
@ -365,35 +370,34 @@ void PairCoulCTIP::read_file(char *file)
MPI_Bcast(&maxparam, 1, MPI_INT, 0, world);
if (comm->me != 0) {
params = (Param *) memory->srealloc(params,maxparam*sizeof(Param), "pair:params");
params = (Param *) memory->srealloc(params, maxparam * sizeof(Param), "pair:params");
}
MPI_Bcast(params, maxparam*sizeof(Param), MPI_BYTE, 0, world);
MPI_Bcast(params, maxparam * sizeof(Param), MPI_BYTE, 0, world);
}
/* ---------------------------------------------------------------------- */
void PairCoulCTIP::setup_params()
{
int i,m,n;
int i, m, n;
// set elem1param
memory->destroy(elem1param);
memory->create(elem1param,nelements,"pair:elem1param");
memory->create(elem1param, nelements, "pair:elem1param");
for (i = 0; i < nelements; i++) {
n = -1;
for (m = 0; m < nparams; m++) {
if (i == params[m].ielement) {
if (n >= 0) error->all(FLERR,"Potential file has duplicate entry");
if (n >= 0) error->all(FLERR, "Potential file has duplicate entry");
n = m;
}
}
if (n < 0) error->all(FLERR,"Potential file is missing an entry");
if (n < 0) error->all(FLERR, "Potential file is missing an entry");
elem1param[i] = n;
}
}
/* ----------------------------------------------------------------------
@ -465,60 +469,76 @@ void *PairCoulCTIP::extract(const char *str, int &dim)
dim = 0;
return (void *) &cut_coul;
}
if (strcmp(str,"chi") == 0 && qeq_x) {
if (strcmp(str, "chi") == 0 && qeq_x) {
dim = 1;
for (int i = 1; i <= atom->ntypes; i++)
if (map[i] >= 0) qeq_x[i] = params[map[i]].chi;
else qeq_x[i] = 0.0;
if (map[i] >= 0)
qeq_x[i] = params[map[i]].chi;
else
qeq_x[i] = 0.0;
return (void *) qeq_x;
}
if (strcmp(str,"eta") == 0 && qeq_j) {
if (strcmp(str, "eta") == 0 && qeq_j) {
dim = 1;
for (int i = 1; i <= atom->ntypes; i++)
if (map[i] >= 0) qeq_j[i] = params[map[i]].eta;
else qeq_j[i] = 0.0;
if (map[i] >= 0)
qeq_j[i] = params[map[i]].eta;
else
qeq_j[i] = 0.0;
return (void *) qeq_j;
}
if (strcmp(str,"gamma") == 0 && qeq_g) {
if (strcmp(str, "gamma") == 0 && qeq_g) {
dim = 1;
for (int i = 1; i <= atom->ntypes; i++)
if (map[i] >= 0) qeq_g[i] = params[map[i]].gamma;
else qeq_g[i] = 0.0;
if (map[i] >= 0)
qeq_g[i] = params[map[i]].gamma;
else
qeq_g[i] = 0.0;
return (void *) qeq_g;
}
if (strcmp(str,"zeta") == 0 && qeq_z) {
if (strcmp(str, "zeta") == 0 && qeq_z) {
dim = 1;
for (int i = 1; i <= atom->ntypes; i++)
if (map[i] >= 0) qeq_z[i] = params[map[i]].zeta;
else qeq_z[i] = 0.0;
if (map[i] >= 0)
qeq_z[i] = params[map[i]].zeta;
else
qeq_z[i] = 0.0;
return (void *) qeq_z;
}
if (strcmp(str,"zcore") == 0 && qeq_c) {
if (strcmp(str, "zcore") == 0 && qeq_c) {
dim = 1;
for (int i = 1; i <= atom->ntypes; i++)
if (map[i] >= 0) qeq_c[i] = params[map[i]].zcore;
else qeq_c[i] = 0.0;
if (map[i] >= 0)
qeq_c[i] = params[map[i]].zcore;
else
qeq_c[i] = 0.0;
return (void *) qeq_c;
}
if (strcmp(str,"qmin") == 0 && qeq_q1) {
if (strcmp(str, "qmin") == 0 && qeq_q1) {
dim = 1;
for (int i = 1; i <= atom->ntypes; i++)
if (map[i] >= 0) qeq_q1[i] = params[map[i]].qmin;
else qeq_q1[i] = 0.0;
if (map[i] >= 0)
qeq_q1[i] = params[map[i]].qmin;
else
qeq_q1[i] = 0.0;
return (void *) qeq_q1;
}
if (strcmp(str,"qmax") == 0 && qeq_q2) {
if (strcmp(str, "qmax") == 0 && qeq_q2) {
dim = 1;
for (int i = 1; i <= atom->ntypes; i++)
if (map[i] >= 0) qeq_q2[i] = params[map[i]].qmax;
else qeq_q2[i] = 0.0;
if (map[i] >= 0)
qeq_q2[i] = params[map[i]].qmax;
else
qeq_q2[i] = 0.0;
return (void *) qeq_q2;
}
if (strcmp(str,"omega") == 0 && qeq_w) {
if (strcmp(str, "omega") == 0 && qeq_w) {
dim = 1;
for (int i = 1; i <= atom->ntypes; i++)
if (map[i] >= 0) qeq_w[i] = params[map[i]].omega;
else qeq_w[i] = 0.0;
if (map[i] >= 0)
qeq_w[i] = params[map[i]].omega;
else
qeq_w[i] = 0.0;
return (void *) qeq_w;
}
return nullptr;