diff --git a/src/REAXFF/fix_qeq_rel_reaxff.cpp b/src/REAXFF/fix_qeq_rel_reaxff.cpp index fbdbcbf674..2cad5e0924 100644 --- a/src/REAXFF/fix_qeq_rel_reaxff.cpp +++ b/src/REAXFF/fix_qeq_rel_reaxff.cpp @@ -50,12 +50,12 @@ void FixQEqRelReaxFF::calc_chi_eff() const auto x = (const double *const *) atom->x; const int *type = atom->type; - double dx, dy, dz, dist_sq, overlap, sum_n, sum_d, expa, expb, chia, phia, phib, p, m; + double dx, dy, dz, dist_sq, overlap, sum_n, sum_d, chia, phia, phib; 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*comm_cutoff < dist_cutoff_sq) { + 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", @@ -71,7 +71,6 @@ void FixQEqRelReaxFF::calc_chi_eff() // compute chi_eff for each local atom for (i = 0; i < nn; i++) { - expa = gauss_exp[type[i]]; chia = chi[type[i]]; if (efield->varflag != FixEfield::ATOM) { phia = -factor * (x[i][0] * efield->ex + x[i][1] * efield->ey + x[i][2] * efield->ez); @@ -86,15 +85,12 @@ void FixQEqRelReaxFF::calc_chi_eff() 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); + dist_sq = (dx*dx + dy*dy + dz*dz); 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_sq); + overlap = prefactor[type[i]][type[j]] * exp(-expfactor[type[i]][type[j]] * dist_sq); if (efield->varflag != FixEfield::ATOM) { phib = -factor * (x[j][0] * efield->ex + x[j][1] * efield->ey + x[j][2] * efield->ez); diff --git a/src/REAXFF/fix_qtpie_reaxff.cpp b/src/REAXFF/fix_qtpie_reaxff.cpp index eb6ae9265f..5ef25eb973 100644 --- a/src/REAXFF/fix_qtpie_reaxff.cpp +++ b/src/REAXFF/fix_qtpie_reaxff.cpp @@ -111,6 +111,8 @@ FixQtpieReaxFF::FixQtpieReaxFF(LAMMPS *lmp, int narg, char **arg) : iarg++; } shld = nullptr; + prefactor = nullptr; + expfactor = nullptr; nn = nt = n_cap = 0; nmax = 0; @@ -172,8 +174,10 @@ FixQtpieReaxFF::~FixQtpieReaxFF() FixQtpieReaxFF::deallocate_matrix(); memory->destroy(shld); - memory->destroy(gauss_exp); + memory->destroy(prefactor); + memory->destroy(expfactor); + if (!reaxflag) { memory->destroy(chi); memory->destroy(eta); @@ -502,6 +506,7 @@ void FixQtpieReaxFF::init() init_shielding(); init_taper(); + init_olap(); if (utils::strmatch(update->integrate_style,"^respa")) nlevels_respa = (dynamic_cast(update->integrate))->nlevels; @@ -569,6 +574,31 @@ void FixQtpieReaxFF::init_taper() /* ---------------------------------------------------------------------- */ +void FixQtpieReaxFF::init_olap() +{ + int i,j; + int ntypes; + double expa,expb,expsum,expnorm; + + ntypes = atom->ntypes; + if (prefactor == nullptr) + memory->create(prefactor,ntypes+1,ntypes+1,"qtpie:overlap_prefactor"); + if (expfactor == nullptr) + memory->create(expfactor,ntypes+1,ntypes+1,"qtpie:overlap_expfactor"); + + for (i = 1; i <= ntypes; ++i) + for (j = 1; j <= ntypes; ++j) { + expa = gauss_exp[i]; + expb = gauss_exp[j]; + expsum = expa + expb; + expnorm = expa * expb / expsum; + prefactor[i][j] = pow((4.0 * expnorm / expsum), 0.75); + expfactor[i][j] = expnorm; + } +} + +/* ---------------------------------------------------------------------- */ + void FixQtpieReaxFF::setup_pre_force(int vflag) { if (reaxff) { @@ -1131,7 +1161,7 @@ void FixQtpieReaxFF::calc_chi_eff() const auto x = (const double * const *)atom->x; const int *type = atom->type; - double dx,dy,dz,dist_sq,overlap,sum_n,sum_d,expa,expb,chia,chib,phia,phib,p,m; + double dx,dy,dz,dist_sq,overlap,sum_n,sum_d,chia,chib,phia,phib; int i,j; // check ghost atoms are stored up to the distance cutoff for overlap integrals @@ -1154,7 +1184,6 @@ void FixQtpieReaxFF::calc_chi_eff() // compute chi_eff for each local atom for (i = 0; i < nn; i++) { - expa = gauss_exp[type[i]]; chia = chi[type[i]]; if (efield) { if (efield->varflag != FixEfield::ATOM) { @@ -1174,13 +1203,10 @@ void FixQtpieReaxFF::calc_chi_eff() dist_sq = (dx*dx + dy*dy + dz*dz); 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_sq); + overlap = prefactor[type[i]][type[j]] * exp(-expfactor[type[i]][type[j]] * dist_sq); if (efield) { if (efield->varflag != FixEfield::ATOM) { diff --git a/src/REAXFF/fix_qtpie_reaxff.h b/src/REAXFF/fix_qtpie_reaxff.h index ae29aaf375..32c3fe5508 100644 --- a/src/REAXFF/fix_qtpie_reaxff.h +++ b/src/REAXFF/fix_qtpie_reaxff.h @@ -90,12 +90,15 @@ 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 **prefactor; // factor used in computation of overlap integrals + double **expfactor; // factor used in exponential term of 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 *); void init_shielding(); void init_taper(); + void init_olap(); void allocate_storage(); void deallocate_storage(); void reallocate_storage();