tx258/lammps-gran-kokkos
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Remove atomics units and distance function

Browse Source
This commit is contained in:
Navraj Lalli
2025-03-31 17:44:17 +01:00
parent 578b1cf936
commit 084ba674a5
3 changed files with 31 additions and 38 deletions
Download Patch File Download Diff File Expand all files Collapse all files

21
src/REAXFF/fix_qeq_rel_reaxff.cpp
View File

@ -35,8 +35,6 @@
using namespace LAMMPS_NS;
using namespace FixConst;
static constexpr double ANGSTROM_TO_BOHRRADIUS = 1.8897261259;
/* ---------------------------------------------------------------------- */
FixQEqRelReaxFF::FixQEqRelReaxFF(LAMMPS *lmp, int narg, char **arg) : FixQtpieReaxFF(lmp, narg, arg)
@ -52,16 +50,16 @@ void FixQEqRelReaxFF::calc_chi_eff()
const auto x = (const double *const *) atom->x;
const int *type = atom->type;
double dist, overlap, sum_n, sum_d, expa, expb, chia, phia, phib, p, m;
double dx, dy, dz, dist_sq, overlap, sum_n, sum_d, expa, expb, chia, phia, phib, p, m;
int i, j;
// check ghost atoms are stored up to the distance cutoff for overlap integrals
const double comm_cutoff = MAX(neighbor->cutneighmax, comm->cutghostuser);
if (comm_cutoff < dist_cutoff / ANGSTROM_TO_BOHRRADIUS) {
if(comm_cutoff*comm_cutoff < dist_cutoff_sq) {
error->all(FLERR, Error::NOLASTLINE,
"Comm cutoff {} is smaller than distance cutoff {} for overlap integrals in fix {}. "
"Increase accordingly using comm_modify cutoff",
comm_cutoff, dist_cutoff / ANGSTROM_TO_BOHRRADIUS, style);
comm_cutoff, sqrt(dist_cutoff_sq), style);
}
// efield energy is in real units of kcal/mol, factor needed for conversion to eV
@ -85,26 +83,29 @@ void FixQEqRelReaxFF::calc_chi_eff()
sum_d = 0.0;
for (j = 0; j < nt; j++) {
dist = distance(x[i], x[j]) * ANGSTROM_TO_BOHRRADIUS; // in atomic units
dx = x[i][0] - x[j][0];
dy = x[i][1] - x[j][1];
dz = x[i][2] - x[j][2];
dist_sq = (dx*dx + dy*dy + dz*dz);
if (dist < dist_cutoff) {
if (dist_sq < dist_cutoff_sq) {
expb = gauss_exp[type[j]];
// overlap integral of two normalised 1s Gaussian type orbitals
p = expa + expb;
m = expa * expb / p;
overlap = pow((4.0 * m / p), 0.75) * exp(-m * dist * dist);
overlap = pow((4.0 * m / p), 0.75) * exp(-m * dist_sq);
if (efield->varflag != FixEfield::ATOM) {
phib = -factor * (x[j][0] * efield->ex + x[j][1] * efield->ey + x[j][2] * efield->ez);
} else { // atom-style potential from FixEfield
phib = efield->efield[j][3];
}
sum_n += (chia + scale * (phia - phib)) * overlap;
sum_n += phib * overlap;
sum_d += overlap;
}
}
chi_eff[i] = sum_n / sum_d;
chi_eff[i] = chia + scale * (phia - sum_n / sum_d);
}
} else {
for (i = 0; i < nn; i++) { chi_eff[i] = chi[type[i]]; }

45
src/REAXFF/fix_qtpie_reaxff.cpp
View File

@ -51,7 +51,7 @@ using namespace FixConst;
static constexpr double CONV_TO_EV = 14.4;
static constexpr double QSUMSMALL = 0.00001;
static constexpr double ANGSTROM_TO_BOHRRADIUS = 1.8897261259;
static constexpr double ANGSTROM_TO_BOHRRADIUS_SQ = 3.571064831;
static const char cite_fix_qtpie_reax[] =
"fix qtpie/reaxff command: \n\n"
@ -244,7 +244,7 @@ void FixQtpieReaxFF::pertype_parameters(char *arg)
if (exp < 0)
throw TokenizerException("Fix qtpie/reaxff: Invalid orbital exponent in gauss file",
std::to_string(exp));
gauss_exp[itype] = exp;
gauss_exp[itype] = exp * ANGSTROM_TO_BOHRRADIUS_SQ;
}
fclose(fp);
} catch (std::exception &e) {
@ -254,11 +254,11 @@ void FixQtpieReaxFF::pertype_parameters(char *arg)
MPI_Bcast(gauss_exp,ntypes+1,MPI_DOUBLE,0,world);
// define a cutoff distance (in atomic units) beyond which overlap integrals are neglected
// in calc_chi_eff()
const double exp_min = find_min_exp(gauss_exp,ntypes+1);
const int olap_cut = 10; // overlap integrals are neglected if less than pow(10,-olap_cut)
dist_cutoff = sqrt(2*olap_cut/exp_min*log(10.0));
// calculate a cutoff distance to neglect overlap integrals in calc_chi_eff()
// when less than pow(10, -olap_cut)
const double exp_min = find_min_exp(gauss_exp, ntypes+1);
const int olap_cut = 10;
dist_cutoff_sq = 2 * olap_cut * log(10.0) / exp_min;
// read chi, eta and gamma
@ -1131,15 +1131,16 @@ void FixQtpieReaxFF::calc_chi_eff()
const auto x = (const double * const *)atom->x;
const int *type = atom->type;
double dist,overlap,sum_n,sum_d,expa,expb,chia,chib,phia,phib,p,m;
double dx,dy,dz,dist_sq,overlap,sum_n,sum_d,expa,expb,chia,chib,phia,phib,p,m;
int i,j;
// check ghost atoms are stored up to the distance cutoff for overlap integrals
const double comm_cutoff = MAX(neighbor->cutneighmax,comm->cutghostuser);
if(comm_cutoff < dist_cutoff/ANGSTROM_TO_BOHRRADIUS) {
error->all(FLERR,"comm cutoff = {} Angstrom is smaller than distance cutoff = {} Angstrom "
"for overlap integrals in {}. Increase comm cutoff with comm_modify",
comm_cutoff, dist_cutoff/ANGSTROM_TO_BOHRRADIUS, style);
if(comm_cutoff*comm_cutoff < dist_cutoff_sq) {
error->all(FLERR, Error::NOLASTLINE,
"Comm cutoff {} is smaller than distance cutoff {} for overlap integrals in fix {}. "
"Increase accordingly using comm_modify cutoff",
comm_cutoff, sqrt(dist_cutoff_sq), style);
}
// efield energy is in real units of kcal/mol, factor needed for conversion to eV
@ -1167,16 +1168,19 @@ void FixQtpieReaxFF::calc_chi_eff()
sum_d = 0.0;
for (j = 0; j < nt; j++) {
dist = distance(x[i],x[j])*ANGSTROM_TO_BOHRRADIUS; // in atomic units
dx = x[i][0] - x[j][0];
dy = x[i][1] - x[j][1];
dz = x[i][2] - x[j][2];
dist_sq = (dx*dx + dy*dy + dz*dz);
if (dist < dist_cutoff) {
if (dist_sq < dist_cutoff_sq) {
expb = gauss_exp[type[j]];
chib = chi[type[j]];
// overlap integral of two normalised 1s Gaussian type orbitals
p = expa + expb;
m = expa * expb / p;
overlap = pow((4.0*m/p),0.75) * exp(-m*dist*dist);
overlap = pow((4.0 * m / p), 0.75) * exp(-m * dist_sq);
if (efield) {
if (efield->varflag != FixEfield::ATOM) {
@ -1209,14 +1213,3 @@ double FixQtpieReaxFF::find_min_exp(const double *array, const int array_length)
}
return exp_min;
}
/* ---------------------------------------------------------------------- */
double FixQtpieReaxFF::distance(const double *posa, const double *posb)
{
double dx, dy, dz;
dx = posb[0] - posa[0];
dy = posb[1] - posa[1];
dz = posb[2] - posa[2];
return sqrt(dx*dx + dy*dy + dz*dz);
}

3
src/REAXFF/fix_qtpie_reaxff.h
View File

@ -90,7 +90,7 @@ class FixQtpieReaxFF : public Fix {
char *pertype_option; // argument to determine how per-type info is obtained
char *gauss_file; // input file for gaussian orbital exponents
double *gauss_exp; // array of gaussian orbital exponents for each atom type
double dist_cutoff; // separation distance beyond which to neglect overlap integrals
double dist_cutoff_sq; // separation distance squared beyond which overlap integrals are neglected
double scale; // scaling factor for electric polarization effects
void pertype_parameters(char *);
@ -131,7 +131,6 @@ class FixQtpieReaxFF : public Fix {
virtual void calc_chi_eff();
double find_min_exp(const double*, const int);
double distance(const double*, const double*);
int matvecs_s, matvecs_t; // Iteration count for each system
};