From d9151d745abb3dd2b0e119805624aa577cf87e4d Mon Sep 17 00:00:00 2001 From: Navraj Lalli Date: Thu, 1 Aug 2024 17:11:06 +0100 Subject: [PATCH 01/40] Create fix qtpie/reaxff --- src/REAXFF/fix_qtpie_reaxff.cpp | 1160 +++++++++++++++++++++++++++++++ src/REAXFF/fix_qtpie_reaxff.h | 149 ++++ src/REAXFF/pair_reaxff.cpp | 6 +- src/fix_efield.h | 1 + 4 files changed, 1314 insertions(+), 2 deletions(-) create mode 100644 src/REAXFF/fix_qtpie_reaxff.cpp create mode 100644 src/REAXFF/fix_qtpie_reaxff.h diff --git a/src/REAXFF/fix_qtpie_reaxff.cpp b/src/REAXFF/fix_qtpie_reaxff.cpp new file mode 100644 index 0000000000..695bdb4316 --- /dev/null +++ b/src/REAXFF/fix_qtpie_reaxff.cpp @@ -0,0 +1,1160 @@ +// clang-format off +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Navraj S Lalli & Efstratios Kritikos (Imperial College London) +------------------------------------------------------------------------- */ + +#include "fix_qtpie_reaxff.h" + +#include "atom.h" +#include "citeme.h" +#include "comm.h" +#include "domain.h" +#include "error.h" +#include "fix_efield.h" +#include "force.h" +#include "group.h" +#include "memory.h" +#include "modify.h" +#include "neigh_list.h" +#include "neighbor.h" +#include "pair.h" +#include "region.h" +#include "respa.h" +#include "text_file_reader.h" +#include "update.h" + +#include "pair_reaxff.h" +#include "reaxff_api.h" + +#include +#include +#include + +using namespace LAMMPS_NS; +using namespace FixConst; + +static constexpr double EV_TO_KCAL_PER_MOL = 14.4; +static constexpr double SMALL = 1.0e-14; +static constexpr double QSUMSMALL = 0.00001; + +static const char cite_fix_qtpie_reaxff[] = + "fix qtpie/reaxff command: doi:https://doi.org/10.1016/j.cplett.2007.02.065\n\n" + "@article{chen2007qtpie,\n" + "title={QTPIE: Charge transfer with polarization current equalization. A fluctuating charge model with correct asymptotics},\n" + "author={Chen, Jiahao and Martinez, Todd J},\n" + "journal={Chemical physics letters},\n" + "volume={438},\n" + "number={4-6},\n" + "pages={315--320},\n" + "year={2007},\n" + "publisher={Elsevier}\n" + "}\n\n"; + +/* ---------------------------------------------------------------------- */ + +FixQtpieReaxFF::FixQtpieReaxFF(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg), matvecs(0), pertype_option(nullptr) +{ + scalar_flag = 1; + extscalar = 0; + imax = 200; + maxwarn = 1; + + if ((narg < 9) || (narg > 13)) error->all(FLERR,"Illegal fix {} command", style); + + nevery = utils::inumeric(FLERR,arg[3],false,lmp); + if (nevery <= 0) error->all(FLERR,"Illegal fix {} command", style); + + swa = utils::numeric(FLERR,arg[4],false,lmp); + swb = utils::numeric(FLERR,arg[5],false,lmp); + tolerance = utils::numeric(FLERR,arg[6],false,lmp); + pertype_option = utils::strdup(arg[7]); + + // dual CG support only available for OPENMP variant + // check for compatibility is in Fix::post_constructor() + + dual_enabled = 0; + + int iarg = 9; + while (iarg < narg) { + if (strcmp(arg[iarg],"dual") == 0) dual_enabled = 1; + else if (strcmp(arg[iarg],"nowarn") == 0) maxwarn = 0; + else if (strcmp(arg[iarg],"maxiter") == 0) { + if (iarg+1 > narg-1) + error->all(FLERR,"Illegal fix {} command", style); + imax = utils::numeric(FLERR,arg[iarg+1],false,lmp); + iarg++; + } else error->all(FLERR,"Illegal fix {} command", style); + iarg++; + } + shld = nullptr; + + nn = n_cap = 0; + nmax = 0; + m_fill = m_cap = 0; + pack_flag = 0; + s = nullptr; + t = nullptr; + nprev = 4; + + Hdia_inv = nullptr; + b_s = nullptr; + chi_field = nullptr; + b_t = nullptr; + b_prc = nullptr; + b_prm = nullptr; + + // CG + + p = nullptr; + q = nullptr; + r = nullptr; + d = nullptr; + + // H matrix + + H.firstnbr = nullptr; + H.numnbrs = nullptr; + H.jlist = nullptr; + H.val = nullptr; + + // dual CG support + // Update comm sizes for this fix + + if (dual_enabled) comm_forward = comm_reverse = 2; + else comm_forward = comm_reverse = 1; + + // perform initial allocation of atom-based arrays + // register with Atom class + + reaxff = dynamic_cast(force->pair_match("^reaxff",0)); + + s_hist = t_hist = nullptr; + atom->add_callback(Atom::GROW); +} + +/* ---------------------------------------------------------------------- */ + +FixQtpieReaxFF::~FixQtpieReaxFF() +{ + if (copymode) return; + + delete[] pertype_option; + + // unregister callbacks to this fix from Atom class + + atom->delete_callback(id,Atom::GROW); + + memory->destroy(s_hist); + memory->destroy(t_hist); + + FixQtpieReaxFF::deallocate_storage(); + FixQtpieReaxFF::deallocate_matrix(); + + memory->destroy(shld); + + if (!reaxflag) { + memory->destroy(chi); + memory->destroy(eta); + memory->destroy(gamma); + } +} + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::post_constructor() +{ + if (lmp->citeme) lmp->citeme->add(cite_fix_qtpie_reaxff); + + grow_arrays(atom->nmax); + for (int i = 0; i < atom->nmax; i++) + for (int j = 0; j < nprev; ++j) + s_hist[i][j] = t_hist[i][j] = 0; + + pertype_parameters(pertype_option); + if (dual_enabled) + error->all(FLERR,"Dual keyword only supported with fix qeq/reaxff/omp"); +} + +/* ---------------------------------------------------------------------- */ + +int FixQtpieReaxFF::setmask() +{ + int mask = 0; + mask |= PRE_FORCE; + mask |= PRE_FORCE_RESPA; + mask |= MIN_PRE_FORCE; + return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::pertype_parameters(char *arg) +{ + const int nlocal = atom->nlocal; + const int *mask = atom->mask; + const int *type = atom->type; + + if (utils::strmatch(arg,"^reaxff")) { + reaxflag = 1; + Pair *pair = force->pair_match("^reaxff",0); + if (!pair) error->all(FLERR,"No reaxff pair style for fix qeq/reaxff"); + + int tmp, tmp_all; + chi = (double *) pair->extract("chi",tmp); + eta = (double *) pair->extract("eta",tmp); + gamma = (double *) pair->extract("gamma",tmp); + if ((chi == nullptr) || (eta == nullptr) || (gamma == nullptr)) + error->all(FLERR, "Fix qeq/reaxff could not extract all Qtpie parameters from pair reaxff"); + tmp = tmp_all = 0; + for (int i = 0; i < nlocal; ++i) { + if (mask[i] & groupbit) { + if ((chi[type[i]] == 0.0) && (eta[type[i]] == 0.0) && (gamma[type[i]] == 0.0)) + tmp = type[i]; + } + } + MPI_Allreduce(&tmp, &tmp_all, 1, MPI_INT, MPI_MAX, world); + if (tmp_all) + error->all(FLERR, "No Qtpie parameters for atom type {} provided by pair reaxff", tmp_all); + return; + } else if (utils::strmatch(arg,"^reax/c")) { + error->all(FLERR, "Fix qeq/reaxff keyword 'reax/c' is obsolete; please use 'reaxff'"); + } else if (platform::file_is_readable(arg)) { + ; // arg is readable file. will read below + } else { + error->all(FLERR, "Unknown fix qeq/reaxff keyword {}", arg); + } + + reaxflag = 0; + + const int ntypes = atom->ntypes; + memory->create(chi,ntypes+1,"qeq/reaxff:chi"); + memory->create(eta,ntypes+1,"qeq/reaxff:eta"); + memory->create(gamma,ntypes+1,"qeq/reaxff:gamma"); + + if (comm->me == 0) { + chi[0] = eta[0] = gamma[0] = 0.0; + try { + TextFileReader reader(arg,"qeq/reaxff parameter"); + reader.ignore_comments = false; + for (int i = 1; i <= ntypes; i++) { + const char *line = reader.next_line(); + if (!line) + throw TokenizerException("Fix qeq/reaxff: Invalid param file format",""); + ValueTokenizer values(line); + + if (values.count() != 4) + throw TokenizerException("Fix qeq/reaxff: Incorrect format of param file",""); + + int itype = values.next_int(); + if ((itype < 1) || (itype > ntypes)) + throw TokenizerException("Fix qeq/reaxff: invalid atom type in param file", + std::to_string(itype)); + + chi[itype] = values.next_double(); + eta[itype] = values.next_double(); + gamma[itype] = values.next_double(); + } + } catch (std::exception &e) { + error->one(FLERR,e.what()); + } + } + + MPI_Bcast(chi,ntypes+1,MPI_DOUBLE,0,world); + MPI_Bcast(eta,ntypes+1,MPI_DOUBLE,0,world); + MPI_Bcast(gamma,ntypes+1,MPI_DOUBLE,0,world); +} + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::allocate_storage() +{ + nmax = atom->nmax; + + memory->create(s,nmax,"qeq:s"); + memory->create(t,nmax,"qeq:t"); + + memory->create(Hdia_inv,nmax,"qeq:Hdia_inv"); + memory->create(b_s,nmax,"qeq:b_s"); + memory->create(chi_field,nmax,"qeq:chi_field"); + memory->create(b_t,nmax,"qeq:b_t"); + memory->create(b_prc,nmax,"qeq:b_prc"); + memory->create(b_prm,nmax,"qeq:b_prm"); + + // dual CG support + int size = nmax; + if (dual_enabled) size*= 2; + + memory->create(p,size,"qeq:p"); + memory->create(q,size,"qeq:q"); + memory->create(r,size,"qeq:r"); + memory->create(d,size,"qeq:d"); +} + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::deallocate_storage() +{ + memory->destroy(s); + memory->destroy(t); + + memory->destroy(Hdia_inv); + memory->destroy(b_s); + memory->destroy(b_t); + memory->destroy(b_prc); + memory->destroy(b_prm); + memory->destroy(chi_field); + + memory->destroy(p); + memory->destroy(q); + memory->destroy(r); + memory->destroy(d); +} + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::reallocate_storage() +{ + deallocate_storage(); + allocate_storage(); + init_storage(); +} + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::allocate_matrix() +{ + int i,ii,m; + + int mincap; + double safezone; + + if (reaxflag) { + mincap = reaxff->api->system->mincap; + safezone = reaxff->api->system->safezone; + } else { + mincap = REAX_MIN_CAP; + safezone = REAX_SAFE_ZONE; + } + + n_cap = MAX((int)(atom->nlocal * safezone), mincap); + + // determine the total space for the H matrix + + m = 0; + for (ii = 0; ii < nn; ii++) { + i = ilist[ii]; + m += numneigh[i]; + } + m_cap = MAX((int)(m * safezone), mincap * REAX_MIN_NBRS); + + H.n = n_cap; + H.m = m_cap; + memory->create(H.firstnbr,n_cap,"qeq:H.firstnbr"); + memory->create(H.numnbrs,n_cap,"qeq:H.numnbrs"); + memory->create(H.jlist,m_cap,"qeq:H.jlist"); + memory->create(H.val,m_cap,"qeq:H.val"); +} + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::deallocate_matrix() +{ + memory->destroy(H.firstnbr); + memory->destroy(H.numnbrs); + memory->destroy(H.jlist); + memory->destroy(H.val); +} + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::reallocate_matrix() +{ + deallocate_matrix(); + allocate_matrix(); +} + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::init() +{ + if (!atom->q_flag) + error->all(FLERR,"Fix {} requires atom attribute q", style); + + if (group->count(igroup) == 0) + error->all(FLERR,"Fix {} group has no atoms", style); + + // compute net charge and print warning if too large + + double qsum_local = 0.0, qsum = 0.0; + for (int i = 0; i < atom->nlocal; i++) { + if (atom->mask[i] & groupbit) + qsum_local += atom->q[i]; + } + MPI_Allreduce(&qsum_local,&qsum,1,MPI_DOUBLE,MPI_SUM,world); + + if ((comm->me == 0) && (fabs(qsum) > QSUMSMALL)) + error->warning(FLERR,"Fix {} group is not charge neutral, net charge = {:.8}", style, qsum); + + // get pointer to fix efield if present. there may be at most one instance of fix efield in use. + + efield = nullptr; + auto fixes = modify->get_fix_by_style("^efield"); + if (fixes.size() == 1) efield = dynamic_cast(fixes.front()); + else if (fixes.size() > 1) + error->all(FLERR, "There may be only one fix efield instance used with fix {}", style); + + // ensure that fix efield is properly initialized before accessing its data and check some settings + if (efield) { + efield->init(); + if (strcmp(update->unit_style,"real") != 0) + error->all(FLERR,"Must use unit_style real with fix {} and external fields", style); + + if (efield->varflag == FixEfield::ATOM && efield->pstyle != FixEfield::ATOM) + error->all(FLERR,"Atom-style external electric field requires atom-style " + "potential variable when used with fix {}", style); + if (((efield->xstyle != FixEfield::CONSTANT) && domain->xperiodic) || + ((efield->ystyle != FixEfield::CONSTANT) && domain->yperiodic) || + ((efield->zstyle != FixEfield::CONSTANT) && domain->zperiodic)) + error->all(FLERR,"Must not have electric field component in direction of periodic " + "boundary when using charge equilibration with ReaxFF."); + if (((fabs(efield->ex) > SMALL) && domain->xperiodic) || + ((fabs(efield->ey) > SMALL) && domain->yperiodic) || + ((fabs(efield->ez) > SMALL) && domain->zperiodic)) + error->all(FLERR,"Must not have electric field component in direction of periodic " + "boundary when using charge equilibration with ReaxFF."); + } + + // we need a half neighbor list w/ Newton off + // built whenever re-neighboring occurs + + neighbor->add_request(this, NeighConst::REQ_NEWTON_OFF); + + init_shielding(); + init_taper(); + + if (utils::strmatch(update->integrate_style,"^respa")) + nlevels_respa = (dynamic_cast(update->integrate))->nlevels; +} + +/* ---------------------------------------------------------------------- */ + +double FixQtpieReaxFF::compute_scalar() +{ + return matvecs/2.0; +} + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::init_list(int /*id*/, NeighList *ptr) +{ + list = ptr; +} + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::init_shielding() +{ + int i,j; + int ntypes; + + ntypes = atom->ntypes; + if (shld == nullptr) + memory->create(shld,ntypes+1,ntypes+1,"qeq:shielding"); + + for (i = 1; i <= ntypes; ++i) + for (j = 1; j <= ntypes; ++j) + shld[i][j] = pow(gamma[i] * gamma[j], -1.5); +} + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::init_taper() +{ + double d7, swa2, swa3, swb2, swb3; + + if (fabs(swa) > 0.01 && comm->me == 0) + error->warning(FLERR,"Fix qeq/reaxff has non-zero lower Taper radius cutoff"); + if (swb < 0) + error->all(FLERR, "Fix qeq/reaxff has negative upper Taper radius cutoff"); + else if (swb < 5 && comm->me == 0) + error->warning(FLERR,"Fix qeq/reaxff has very low Taper radius cutoff"); + + d7 = pow(swb - swa, 7); + swa2 = SQR(swa); + swa3 = CUBE(swa); + swb2 = SQR(swb); + swb3 = CUBE(swb); + + Tap[7] = 20.0 / d7; + Tap[6] = -70.0 * (swa + swb) / d7; + Tap[5] = 84.0 * (swa2 + 3.0*swa*swb + swb2) / d7; + Tap[4] = -35.0 * (swa3 + 9.0*swa2*swb + 9.0*swa*swb2 + swb3) / d7; + Tap[3] = 140.0 * (swa3*swb + 3.0*swa2*swb2 + swa*swb3) / d7; + Tap[2] =-210.0 * (swa3*swb2 + swa2*swb3) / d7; + Tap[1] = 140.0 * swa3 * swb3 / d7; + Tap[0] = (-35.0*swa3*swb2*swb2 + 21.0*swa2*swb3*swb2 - + 7.0*swa*swb3*swb3 + swb3*swb3*swb) / d7; +} + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::setup_pre_force(int vflag) +{ + if (reaxff) { + nn = reaxff->list->inum; + ilist = reaxff->list->ilist; + numneigh = reaxff->list->numneigh; + firstneigh = reaxff->list->firstneigh; + } else { + nn = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + } + + deallocate_storage(); + allocate_storage(); + + init_storage(); + + deallocate_matrix(); + allocate_matrix(); + + pre_force(vflag); +} + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::setup_pre_force_respa(int vflag, int ilevel) +{ + if (ilevel < nlevels_respa-1) return; + setup_pre_force(vflag); +} + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::min_setup_pre_force(int vflag) +{ + setup_pre_force(vflag); +} + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::init_storage() +{ + if (efield) get_chi_field(); + + for (int ii = 0; ii < nn; ii++) { + int i = ilist[ii]; + if (atom->mask[i] & groupbit) { + Hdia_inv[i] = 1. / eta[atom->type[i]]; + b_s[i] = -chi[atom->type[i]]; + if (efield) b_s[i] -= chi_field[i]; + b_t[i] = -1.0; + b_prc[i] = 0; + b_prm[i] = 0; + s[i] = t[i] = 0; + } + } +} + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::pre_force(int /*vflag*/) +{ + if (update->ntimestep % nevery) return; + + int n = atom->nlocal; + + if (reaxff) { + nn = reaxff->list->inum; + ilist = reaxff->list->ilist; + numneigh = reaxff->list->numneigh; + firstneigh = reaxff->list->firstneigh; + } else { + nn = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + } + + // grow arrays if necessary + // need to be atom->nmax in length + + if (atom->nmax > nmax) reallocate_storage(); + if (n > n_cap*DANGER_ZONE || m_fill > m_cap*DANGER_ZONE) + reallocate_matrix(); + + if (efield) get_chi_field(); + + init_matvec(); + + matvecs_s = CG(b_s, s); // CG on s - parallel + matvecs_t = CG(b_t, t); // CG on t - parallel + matvecs = matvecs_s + matvecs_t; + + calculate_Q(); +} + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::pre_force_respa(int vflag, int ilevel, int /*iloop*/) +{ + if (ilevel == nlevels_respa-1) pre_force(vflag); +} + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::min_pre_force(int vflag) +{ + pre_force(vflag); +} + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::init_matvec() +{ + /* fill-in H matrix */ + compute_H(); + + int ii, i; + + for (ii = 0; ii < nn; ++ii) { + i = ilist[ii]; + if (atom->mask[i] & groupbit) { + + /* init pre-conditioner for H and init solution vectors */ + Hdia_inv[i] = 1. / eta[atom->type[i]]; + b_s[i] = -chi[atom->type[i]]; + if (efield) b_s[i] -= chi_field[i]; + b_t[i] = -1.0; + + /* quadratic extrapolation for s & t from previous solutions */ + t[i] = t_hist[i][2] + 3 * (t_hist[i][0] - t_hist[i][1]); + + /* cubic extrapolation for s & t from previous solutions */ + s[i] = 4*(s_hist[i][0]+s_hist[i][2])-(6*s_hist[i][1]+s_hist[i][3]); + } + } + + pack_flag = 2; + comm->forward_comm(this); //Dist_vector(s); + pack_flag = 3; + comm->forward_comm(this); //Dist_vector(t); +} + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::compute_H() +{ + int jnum; + int i, j, ii, jj, flag; + double dx, dy, dz, r_sqr; + constexpr double EPSILON = 0.0001; + + int *type = atom->type; + tagint *tag = atom->tag; + double **x = atom->x; + int *mask = atom->mask; + + // fill in the H matrix + m_fill = 0; + r_sqr = 0; + for (ii = 0; ii < nn; ii++) { + i = ilist[ii]; + if (mask[i] & groupbit) { + jlist = firstneigh[i]; + jnum = numneigh[i]; + H.firstnbr[i] = m_fill; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + + dx = x[j][0] - x[i][0]; + dy = x[j][1] - x[i][1]; + dz = x[j][2] - x[i][2]; + r_sqr = SQR(dx) + SQR(dy) + SQR(dz); + + flag = 0; + if (r_sqr <= SQR(swb)) { + if (j < atom->nlocal) flag = 1; + else if (tag[i] < tag[j]) flag = 1; + else if (tag[i] == tag[j]) { + if (dz > EPSILON) flag = 1; + else if (fabs(dz) < EPSILON) { + if (dy > EPSILON) flag = 1; + else if (fabs(dy) < EPSILON && dx > EPSILON) + flag = 1; + } + } + } + + if (flag) { + H.jlist[m_fill] = j; + H.val[m_fill] = calculate_H(sqrt(r_sqr), shld[type[i]][type[j]]); + m_fill++; + } + } + H.numnbrs[i] = m_fill - H.firstnbr[i]; + } + } + + if (m_fill >= H.m) + error->all(FLERR,"Fix qeq/reaxff H matrix size has been exceeded: m_fill={} H.m={}\n", + m_fill, H.m); +} + +/* ---------------------------------------------------------------------- */ + +double FixQtpieReaxFF::calculate_H(double r, double gamma) +{ + double Taper, denom; + + Taper = Tap[7] * r + Tap[6]; + Taper = Taper * r + Tap[5]; + Taper = Taper * r + Tap[4]; + Taper = Taper * r + Tap[3]; + Taper = Taper * r + Tap[2]; + Taper = Taper * r + Tap[1]; + Taper = Taper * r + Tap[0]; + + denom = r * r * r + gamma; + denom = pow(denom,1.0/3.0); + + return Taper * EV_TO_KCAL_PER_MOL / denom; +} + +/* ---------------------------------------------------------------------- */ + +int FixQtpieReaxFF::CG(double *b, double *x) +{ + int i, j; + double tmp, alpha, beta, b_norm; + double sig_old, sig_new; + + int jj; + + pack_flag = 1; + sparse_matvec(&H, x, q); + comm->reverse_comm(this); //Coll_Vector(q); + + vector_sum(r , 1., b, -1., q, nn); + + for (jj = 0; jj < nn; ++jj) { + j = ilist[jj]; + if (atom->mask[j] & groupbit) + d[j] = r[j] * Hdia_inv[j]; //pre-condition + } + + b_norm = parallel_norm(b, nn); + sig_new = parallel_dot(r, d, nn); + + for (i = 1; i < imax && sqrt(sig_new) / b_norm > tolerance; ++i) { + comm->forward_comm(this); //Dist_vector(d); + sparse_matvec(&H, d, q); + comm->reverse_comm(this); //Coll_vector(q); + + tmp = parallel_dot(d, q, nn); + alpha = sig_new / tmp; + + vector_add(x, alpha, d, nn); + vector_add(r, -alpha, q, nn); + + // pre-conditioning + for (jj = 0; jj < nn; ++jj) { + j = ilist[jj]; + if (atom->mask[j] & groupbit) + p[j] = r[j] * Hdia_inv[j]; + } + + sig_old = sig_new; + sig_new = parallel_dot(r, p, nn); + + beta = sig_new / sig_old; + vector_sum(d, 1., p, beta, d, nn); + } + + if ((i >= imax) && maxwarn && (comm->me == 0)) + error->warning(FLERR, "Fix qeq/reaxff CG convergence failed after {} iterations at step {}", + i,update->ntimestep); + return i; +} + + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::sparse_matvec(sparse_matrix *A, double *x, double *b) +{ + int i, j, itr_j; + int ii; + + for (ii = 0; ii < nn; ++ii) { + i = ilist[ii]; + if (atom->mask[i] & groupbit) + b[i] = eta[atom->type[i]] * x[i]; + } + + int nall = atom->nlocal + atom->nghost; + for (i = atom->nlocal; i < nall; ++i) + b[i] = 0; + + for (ii = 0; ii < nn; ++ii) { + i = ilist[ii]; + if (atom->mask[i] & groupbit) { + for (itr_j=A->firstnbr[i]; itr_jfirstnbr[i]+A->numnbrs[i]; itr_j++) { + j = A->jlist[itr_j]; + b[i] += A->val[itr_j] * x[j]; + b[j] += A->val[itr_j] * x[i]; + } + } + } + +} + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::calculate_Q() +{ + int i, k; + double u, s_sum, t_sum; + double *q = atom->q; + + int ii; + + s_sum = parallel_vector_acc(s, nn); + t_sum = parallel_vector_acc(t, nn); + u = s_sum / t_sum; + + for (ii = 0; ii < nn; ++ii) { + i = ilist[ii]; + if (atom->mask[i] & groupbit) { + q[i] = s[i] - u * t[i]; + + /* backup s & t */ + for (k = nprev-1; k > 0; --k) { + s_hist[i][k] = s_hist[i][k-1]; + t_hist[i][k] = t_hist[i][k-1]; + } + s_hist[i][0] = s[i]; + t_hist[i][0] = t[i]; + } + } + + pack_flag = 4; + comm->forward_comm(this); //Dist_vector(atom->q); +} + +/* ---------------------------------------------------------------------- */ + +int FixQtpieReaxFF::pack_forward_comm(int n, int *list, double *buf, + int /*pbc_flag*/, int * /*pbc*/) +{ + int m; + + if (pack_flag == 1) + for (m = 0; m < n; m++) buf[m] = d[list[m]]; + else if (pack_flag == 2) + for (m = 0; m < n; m++) buf[m] = s[list[m]]; + else if (pack_flag == 3) + for (m = 0; m < n; m++) buf[m] = t[list[m]]; + else if (pack_flag == 4) + for (m = 0; m < n; m++) buf[m] = atom->q[list[m]]; + else if (pack_flag == 5) { + m = 0; + for (int i = 0; i < n; i++) { + int j = 2 * list[i]; + buf[m++] = d[j]; + buf[m++] = d[j+1]; + } + return m; + } + return n; +} + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::unpack_forward_comm(int n, int first, double *buf) +{ + int i, m; + + if (pack_flag == 1) + for (m = 0, i = first; m < n; m++, i++) d[i] = buf[m]; + else if (pack_flag == 2) + for (m = 0, i = first; m < n; m++, i++) s[i] = buf[m]; + else if (pack_flag == 3) + for (m = 0, i = first; m < n; m++, i++) t[i] = buf[m]; + else if (pack_flag == 4) + for (m = 0, i = first; m < n; m++, i++) atom->q[i] = buf[m]; + else if (pack_flag == 5) { + int last = first + n; + m = 0; + for (i = first; i < last; i++) { + int j = 2 * i; + d[j] = buf[m++]; + d[j+1] = buf[m++]; + } + } +} + +/* ---------------------------------------------------------------------- */ + +int FixQtpieReaxFF::pack_reverse_comm(int n, int first, double *buf) +{ + int i, m; + if (pack_flag == 5) { + m = 0; + int last = first + n; + for (i = first; i < last; i++) { + int indxI = 2 * i; + buf[m++] = q[indxI]; + buf[m++] = q[indxI+1]; + } + return m; + } else { + for (m = 0, i = first; m < n; m++, i++) buf[m] = q[i]; + return n; + } +} + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::unpack_reverse_comm(int n, int *list, double *buf) +{ + if (pack_flag == 5) { + int m = 0; + for (int i = 0; i < n; i++) { + int indxI = 2 * list[i]; + q[indxI] += buf[m++]; + q[indxI+1] += buf[m++]; + } + } else { + for (int m = 0; m < n; m++) q[list[m]] += buf[m]; + } +} + +/* ---------------------------------------------------------------------- + memory usage of local atom-based arrays +------------------------------------------------------------------------- */ + +double FixQtpieReaxFF::memory_usage() +{ + double bytes; + + bytes = (double)atom->nmax*nprev*2 * sizeof(double); // s_hist & t_hist + bytes += (double)atom->nmax*11 * sizeof(double); // storage + bytes += (double)n_cap*2 * sizeof(int); // matrix... + bytes += (double)m_cap * sizeof(int); + bytes += (double)m_cap * sizeof(double); + + if (dual_enabled) + bytes += (double)atom->nmax*4 * sizeof(double); // double size for q, d, r, and p + + return bytes; +} + +/* ---------------------------------------------------------------------- + allocate fictitious charge arrays +------------------------------------------------------------------------- */ + +void FixQtpieReaxFF::grow_arrays(int nmax) +{ + memory->grow(s_hist,nmax,nprev,"qeq:s_hist"); + memory->grow(t_hist,nmax,nprev,"qeq:t_hist"); +} + +/* ---------------------------------------------------------------------- + copy values within fictitious charge arrays +------------------------------------------------------------------------- */ + +void FixQtpieReaxFF::copy_arrays(int i, int j, int /*delflag*/) +{ + for (int m = 0; m < nprev; m++) { + s_hist[j][m] = s_hist[i][m]; + t_hist[j][m] = t_hist[i][m]; + } +} + +/* ---------------------------------------------------------------------- + pack values in local atom-based array for exchange with another proc +------------------------------------------------------------------------- */ + +int FixQtpieReaxFF::pack_exchange(int i, double *buf) +{ + for (int m = 0; m < nprev; m++) buf[m] = s_hist[i][m]; + for (int m = 0; m < nprev; m++) buf[nprev+m] = t_hist[i][m]; + return nprev*2; +} + +/* ---------------------------------------------------------------------- + unpack values in local atom-based array from exchange with another proc +------------------------------------------------------------------------- */ + +int FixQtpieReaxFF::unpack_exchange(int nlocal, double *buf) +{ + for (int m = 0; m < nprev; m++) s_hist[nlocal][m] = buf[m]; + for (int m = 0; m < nprev; m++) t_hist[nlocal][m] = buf[nprev+m]; + return nprev*2; +} + +/* ---------------------------------------------------------------------- */ + +double FixQtpieReaxFF::parallel_norm(double *v, int n) +{ + int i; + double my_sum, norm_sqr; + + int ii; + + my_sum = 0.0; + norm_sqr = 0.0; + for (ii = 0; ii < n; ++ii) { + i = ilist[ii]; + if (atom->mask[i] & groupbit) + my_sum += SQR(v[i]); + } + + MPI_Allreduce(&my_sum, &norm_sqr, 1, MPI_DOUBLE, MPI_SUM, world); + + return sqrt(norm_sqr); +} + +/* ---------------------------------------------------------------------- */ + +double FixQtpieReaxFF::parallel_dot(double *v1, double *v2, int n) +{ + int i; + double my_dot, res; + + int ii; + + my_dot = 0.0; + res = 0.0; + for (ii = 0; ii < n; ++ii) { + i = ilist[ii]; + if (atom->mask[i] & groupbit) + my_dot += v1[i] * v2[i]; + } + + MPI_Allreduce(&my_dot, &res, 1, MPI_DOUBLE, MPI_SUM, world); + + return res; +} + +/* ---------------------------------------------------------------------- */ + +double FixQtpieReaxFF::parallel_vector_acc(double *v, int n) +{ + int i; + double my_acc, res; + + int ii; + + my_acc = 0.0; + res = 0.0; + for (ii = 0; ii < n; ++ii) { + i = ilist[ii]; + if (atom->mask[i] & groupbit) + my_acc += v[i]; + } + + MPI_Allreduce(&my_acc, &res, 1, MPI_DOUBLE, MPI_SUM, world); + + return res; +} + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::vector_sum(double* dest, double c, double* v, + double d, double* y, int k) +{ + int kk; + + for (--k; k>=0; --k) { + kk = ilist[k]; + if (atom->mask[kk] & groupbit) + dest[kk] = c * v[kk] + d * y[kk]; + } +} + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::vector_add(double* dest, double c, double* v, int k) +{ + int kk; + + for (--k; k>=0; --k) { + kk = ilist[k]; + if (atom->mask[kk] & groupbit) + dest[kk] += c * v[kk]; + } +} + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::get_chi_field() +{ + memset(&chi_field[0],0,atom->nmax*sizeof(double)); + if (!efield) return; + + const auto x = (const double * const *)atom->x; + const int *mask = atom->mask; + const imageint *image = atom->image; + const int nlocal = atom->nlocal; + + + // update electric field region if necessary + + Region *region = efield->region; + if (region) region->prematch(); + + // efield energy is in real units of kcal/mol/angstrom, need to convert to eV + + const double qe2f = force->qe2f; + const double factor = -1.0/qe2f; + + + if (efield->varflag != FixEfield::CONSTANT) + efield->update_efield_variables(); + + // atom selection is for the group of fix efield + + double unwrap[3]; + const double ex = efield->ex; + const double ey = efield->ey; + const double ez = efield->ez; + const int efgroupbit = efield->groupbit; + + // charge interactions + // force = qE, potential energy = F dot x in unwrapped coords + if (efield->varflag != FixEfield::ATOM) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & efgroupbit) { + if (region && !region->match(x[i][0],x[i][1],x[i][2])) continue; + domain->unmap(x[i],image[i],unwrap); + chi_field[i] = factor*(ex*unwrap[0] + ey*unwrap[1] + ez*unwrap[2]); + } + } + } else { // must use atom-style potential from FixEfield + for (int i = 0; i < nlocal; i++) { + if (mask[i] & efgroupbit) { + if (region && !region->match(x[i][0],x[i][1],x[i][2])) continue; + chi_field[i] = -efield->efield[i][3]; + } + } + } +} diff --git a/src/REAXFF/fix_qtpie_reaxff.h b/src/REAXFF/fix_qtpie_reaxff.h new file mode 100644 index 0000000000..fcdedddf26 --- /dev/null +++ b/src/REAXFF/fix_qtpie_reaxff.h @@ -0,0 +1,149 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef FIX_CLASS +// clang-format off +FixStyle(qtpie/reax,FixQtpieReaxFF); +FixStyle(qtpie/reaxff,FixQtpieReaxFF); +// clang-format on +#else + +#ifndef LMP_FIX_QTPIE_REAXFF_H +#define LMP_FIX_QTPIE_REAXFF_H + +#include "fix.h" + +namespace LAMMPS_NS { + +class FixQtpieReaxFF : public Fix { + public: + FixQtpieReaxFF(class LAMMPS *, int, char **); + ~FixQtpieReaxFF() override; + int setmask() override; + void post_constructor() override; + void init() override; + void init_list(int, class NeighList *) override; + virtual void init_storage(); + void setup_pre_force(int) override; + void pre_force(int) override; + + void setup_pre_force_respa(int, int) override; + void pre_force_respa(int, int, int) override; + + void min_setup_pre_force(int); + void min_pre_force(int) override; + + double compute_scalar() override; + + protected: + int nevery, reaxflag; + int matvecs; + int nn, m_fill; + int n_cap, nmax, m_cap; + int pack_flag; + int nlevels_respa; + class NeighList *list; + class PairReaxFF *reaxff; + class FixEfield *efield; + int *ilist, *jlist, *numneigh, **firstneigh; + + double swa, swb; // lower/upper Taper cutoff radius + double Tap[8]; // Taper function + double tolerance; // tolerance for the norm of the rel residual in CG + + double *chi, *eta, *gamma; // qeq parameters + double **shld; + + // fictitious charges + + double *s, *t; + double **s_hist, **t_hist; + int nprev; + + typedef struct { + int n, m; + int *firstnbr; + int *numnbrs; + int *jlist; + double *val; + } sparse_matrix; + + sparse_matrix H; + double *Hdia_inv; + double *b_s, *b_t; + double *b_prc, *b_prm; + double *chi_field; + + //CG storage + double *p, *q, *r, *d; + int imax, maxwarn; + + char *pertype_option; // argument to determine how per-type info is obtained + + // Params from Kritikos - could rename or move to protected later + char *gauss_file; // input file for gaussian exponents for each type of REAXFF file + double cutghost; // ghost atoms cutoff (used for check) + int nn_prev; // number of local atoms; needed for memory reallocation of chi_eff (when multiprocessing) + double *gauss_exp; // array of gaussian exponents + double *chi_eff; // array of effective electronegativities + double *chi_eff_init; // array of effective electronegativities for FixQEqReax::init_storage() + + // void calculate_chi_eff(LAMMPS_NS::Atom *atom, reax_system *system, double *chi, + // int ni, int nj, double *lchi_eff); + virtual void pertype_parameters(char *); + void init_shielding(); + void init_taper(); + virtual void allocate_storage(); + virtual void deallocate_storage(); + void reallocate_storage(); + virtual void allocate_matrix(); + virtual void deallocate_matrix(); + void reallocate_matrix(); + + virtual void init_matvec(); + void init_H(); + virtual void compute_H(); + double calculate_H(double, double); + virtual void calculate_Q(); + + virtual int CG(double *, double *); + virtual void sparse_matvec(sparse_matrix *, double *, double *); + + int pack_forward_comm(int, int *, double *, int, int *) override; + void unpack_forward_comm(int, int, double *) override; + int pack_reverse_comm(int, int, double *) override; + void unpack_reverse_comm(int, int *, double *) override; + double memory_usage() override; + void grow_arrays(int) override; + void copy_arrays(int, int, int) override; + int pack_exchange(int, double *) override; + int unpack_exchange(int, double *) override; + + virtual double parallel_norm(double *, int); + virtual double parallel_dot(double *, double *, int); + virtual double parallel_vector_acc(double *, int); + + virtual void vector_sum(double *, double, double *, double, double *, int); + virtual void vector_add(double *, double, double *, int); + + virtual void get_chi_field(); + + // dual CG support + int dual_enabled; // 0: Original, separate s & t optimization; 1: dual optimization + int matvecs_s, matvecs_t; // Iteration count for each system +}; + +} // namespace LAMMPS_NS + +#endif +#endif diff --git a/src/REAXFF/pair_reaxff.cpp b/src/REAXFF/pair_reaxff.cpp index b9f4f6c838..06ad172a38 100644 --- a/src/REAXFF/pair_reaxff.cpp +++ b/src/REAXFF/pair_reaxff.cpp @@ -339,11 +339,13 @@ void PairReaxFF::init_style() auto acks2_fixes = modify->get_fix_by_style("^acks2/reax"); int have_qeq = modify->get_fix_by_style("^qeq/reax").size() - + modify->get_fix_by_style("^qeq/shielded").size() + acks2_fixes.size(); + + modify->get_fix_by_style("^qeq/shielded").size() + acks2_fixes.size() + + modify->get_fix_by_style("^qtpie/reax").size(); if (qeqflag && (have_qeq != 1)) error->all(FLERR,"Pair style reaxff requires use of exactly one of the " - "fix qeq/reaxff or fix qeq/shielded or fix acks2/reaxff commands"); + "fix qeq/reaxff or fix qeq/shielded or fix acks2/reaxff or " + "fix qtpie/reaxff commands"); api->system->acks2_flag = acks2_fixes.size(); if (api->system->acks2_flag) diff --git a/src/fix_efield.h b/src/fix_efield.h index 72fd204898..108395cc2c 100644 --- a/src/fix_efield.h +++ b/src/fix_efield.h @@ -26,6 +26,7 @@ namespace LAMMPS_NS { class FixEfield : public Fix { friend class FixQEqReaxFF; + friend class FixQtpieReaxFF; public: FixEfield(class LAMMPS *, int, char **); From be43a2bdeb2e614aaaf816cfbf94c7ce0368a0ed Mon Sep 17 00:00:00 2001 From: Navraj Lalli Date: Fri, 2 Aug 2024 18:11:20 +0100 Subject: [PATCH 02/40] Allow for reading of Gaussian exponents from file --- src/REAXFF/fix_qtpie_reaxff.cpp | 53 ++++++++++++++++++++++++++++----- 1 file changed, 45 insertions(+), 8 deletions(-) diff --git a/src/REAXFF/fix_qtpie_reaxff.cpp b/src/REAXFF/fix_qtpie_reaxff.cpp index 695bdb4316..77b4708f7a 100644 --- a/src/REAXFF/fix_qtpie_reaxff.cpp +++ b/src/REAXFF/fix_qtpie_reaxff.cpp @@ -42,6 +42,7 @@ #include #include #include +#include using namespace LAMMPS_NS; using namespace FixConst; @@ -66,7 +67,7 @@ static const char cite_fix_qtpie_reaxff[] = /* ---------------------------------------------------------------------- */ FixQtpieReaxFF::FixQtpieReaxFF(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg), matvecs(0), pertype_option(nullptr) + Fix(lmp, narg, arg), matvecs(0), pertype_option(nullptr), gauss_file(nullptr) { scalar_flag = 1; extscalar = 0; @@ -82,6 +83,7 @@ FixQtpieReaxFF::FixQtpieReaxFF(LAMMPS *lmp, int narg, char **arg) : swb = utils::numeric(FLERR,arg[5],false,lmp); tolerance = utils::numeric(FLERR,arg[6],false,lmp); pertype_option = utils::strdup(arg[7]); + gauss_file = utils::strdup(arg[8]); // dual CG support only available for OPENMP variant // check for compatibility is in Fix::post_constructor() @@ -153,6 +155,7 @@ FixQtpieReaxFF::~FixQtpieReaxFF() if (copymode) return; delete[] pertype_option; + delete[] gauss_file; // unregister callbacks to this fix from Atom class @@ -166,6 +169,7 @@ FixQtpieReaxFF::~FixQtpieReaxFF() memory->destroy(shld); + memory->destroy(gauss_exp); if (!reaxflag) { memory->destroy(chi); memory->destroy(eta); @@ -207,18 +211,52 @@ void FixQtpieReaxFF::pertype_parameters(char *arg) const int nlocal = atom->nlocal; const int *mask = atom->mask; const int *type = atom->type; + const int ntypes = atom->ntypes; + + // read gaussian exponents + memory->create(gauss_exp,ntypes+1,"qtpie/reaxff:gauss_exp"); + if (comm->me == 0) { + gauss_exp[0] = 0.0; + try { + TextFileReader reader(gauss_file,"qtpie/reaxff gaussian exponents"); + reader.ignore_comments = false; + for (int i = 1; i <= ntypes; i++) { + const char *line = reader.next_line(); + std::cout << line; + if (!line) + throw TokenizerException("Fix qtpie/reaxff: Invalid param file format",""); + ValueTokenizer values(line); + + if (values.count() != 2) + throw TokenizerException("Fix qtpie/reaxff: Incorrect format of param file",""); + + int itype = values.next_int(); + if ((itype < 1) || (itype > ntypes)) + throw TokenizerException("Fix qtpie/reaxff: Invalid atom type in param file", + std::to_string(itype)); + + gauss_exp[itype] = values.next_double(); + } + } catch (std::exception &e) { + error->one(FLERR,e.what()); + } + } + + MPI_Bcast(gauss_exp,ntypes+1,MPI_DOUBLE,0,world); + + // read chi, eta and gamma if (utils::strmatch(arg,"^reaxff")) { reaxflag = 1; Pair *pair = force->pair_match("^reaxff",0); - if (!pair) error->all(FLERR,"No reaxff pair style for fix qeq/reaxff"); + if (!pair) error->all(FLERR,"No reaxff pair style for fix qtpie/reaxff"); int tmp, tmp_all; chi = (double *) pair->extract("chi",tmp); eta = (double *) pair->extract("eta",tmp); gamma = (double *) pair->extract("gamma",tmp); if ((chi == nullptr) || (eta == nullptr) || (gamma == nullptr)) - error->all(FLERR, "Fix qeq/reaxff could not extract all Qtpie parameters from pair reaxff"); + error->all(FLERR, "Fix qtpie/reaxff could not extract qtpie parameters from pair reaxff"); tmp = tmp_all = 0; for (int i = 0; i < nlocal; ++i) { if (mask[i] & groupbit) { @@ -228,19 +266,18 @@ void FixQtpieReaxFF::pertype_parameters(char *arg) } MPI_Allreduce(&tmp, &tmp_all, 1, MPI_INT, MPI_MAX, world); if (tmp_all) - error->all(FLERR, "No Qtpie parameters for atom type {} provided by pair reaxff", tmp_all); + error->all(FLERR, "No qtpie parameters for atom type {} provided by pair reaxff", tmp_all); return; } else if (utils::strmatch(arg,"^reax/c")) { - error->all(FLERR, "Fix qeq/reaxff keyword 'reax/c' is obsolete; please use 'reaxff'"); + error->all(FLERR, "Fix qtpie/reaxff keyword 'reax/c' is obsolete; please use 'reaxff'"); } else if (platform::file_is_readable(arg)) { ; // arg is readable file. will read below } else { - error->all(FLERR, "Unknown fix qeq/reaxff keyword {}", arg); + error->all(FLERR, "Unknown fix qtpie/reaxff keyword {}", arg); } reaxflag = 0; - const int ntypes = atom->ntypes; memory->create(chi,ntypes+1,"qeq/reaxff:chi"); memory->create(eta,ntypes+1,"qeq/reaxff:eta"); memory->create(gamma,ntypes+1,"qeq/reaxff:gamma"); @@ -261,7 +298,7 @@ void FixQtpieReaxFF::pertype_parameters(char *arg) int itype = values.next_int(); if ((itype < 1) || (itype > ntypes)) - throw TokenizerException("Fix qeq/reaxff: invalid atom type in param file", + throw TokenizerException("Fix qeq/reaxff: Invalid atom type in param file", std::to_string(itype)); chi[itype] = values.next_double(); From bfb1c64b647219e5927b65c932e96aeb321b5401 Mon Sep 17 00:00:00 2001 From: Navraj Lalli Date: Tue, 6 Aug 2024 12:43:19 +0100 Subject: [PATCH 03/40] Add functionality for calculating chi_eff --- src/REAXFF/fix_qtpie_reaxff.cpp | 142 ++++++++++++++++++++++++++++++-- src/REAXFF/fix_qtpie_reaxff.h | 14 ++-- 2 files changed, 142 insertions(+), 14 deletions(-) diff --git a/src/REAXFF/fix_qtpie_reaxff.cpp b/src/REAXFF/fix_qtpie_reaxff.cpp index 77b4708f7a..9f3714c32d 100644 --- a/src/REAXFF/fix_qtpie_reaxff.cpp +++ b/src/REAXFF/fix_qtpie_reaxff.cpp @@ -104,7 +104,7 @@ FixQtpieReaxFF::FixQtpieReaxFF(LAMMPS *lmp, int narg, char **arg) : } shld = nullptr; - nn = n_cap = 0; + nn = ng = n_cap = 0; nmax = 0; m_fill = m_cap = 0; pack_flag = 0; @@ -115,6 +115,7 @@ FixQtpieReaxFF::FixQtpieReaxFF(LAMMPS *lmp, int narg, char **arg) : Hdia_inv = nullptr; b_s = nullptr; chi_field = nullptr; + chi_eff = nullptr; b_t = nullptr; b_prc = nullptr; b_prm = nullptr; @@ -327,6 +328,7 @@ void FixQtpieReaxFF::allocate_storage() memory->create(Hdia_inv,nmax,"qeq:Hdia_inv"); memory->create(b_s,nmax,"qeq:b_s"); memory->create(chi_field,nmax,"qeq:chi_field"); + memory->create(chi_eff,nmax,"qtpie:chi_eff"); memory->create(b_t,nmax,"qeq:b_t"); memory->create(b_prc,nmax,"qeq:b_prc"); memory->create(b_prm,nmax,"qeq:b_prm"); @@ -354,6 +356,7 @@ void FixQtpieReaxFF::deallocate_storage() memory->destroy(b_prc); memory->destroy(b_prm); memory->destroy(chi_field); + memory->destroy(chi_eff); memory->destroy(p); memory->destroy(q); @@ -553,11 +556,13 @@ void FixQtpieReaxFF::setup_pre_force(int vflag) { if (reaxff) { nn = reaxff->list->inum; + ng = reaxff->list->inum + reaxff->list->gnum; ilist = reaxff->list->ilist; numneigh = reaxff->list->numneigh; firstneigh = reaxff->list->firstneigh; } else { nn = list->inum; + ng = list->inum + list->gnum; ilist = list->ilist; numneigh = list->numneigh; firstneigh = list->firstneigh; @@ -593,14 +598,15 @@ void FixQtpieReaxFF::min_setup_pre_force(int vflag) void FixQtpieReaxFF::init_storage() { - if (efield) get_chi_field(); + // if (efield) get_chi_field(); + calc_chi_eff(); for (int ii = 0; ii < nn; ii++) { int i = ilist[ii]; if (atom->mask[i] & groupbit) { Hdia_inv[i] = 1. / eta[atom->type[i]]; - b_s[i] = -chi[atom->type[i]]; - if (efield) b_s[i] -= chi_field[i]; + b_s[i] = -chi_eff[i]; + // if (efield) b_s[i] -= chi_field[i]; b_t[i] = -1.0; b_prc[i] = 0; b_prm[i] = 0; @@ -619,11 +625,13 @@ void FixQtpieReaxFF::pre_force(int /*vflag*/) if (reaxff) { nn = reaxff->list->inum; + ng = reaxff->list->inum + reaxff->list->gnum; ilist = reaxff->list->ilist; numneigh = reaxff->list->numneigh; firstneigh = reaxff->list->firstneigh; } else { nn = list->inum; + ng = list->inum + list->gnum; ilist = list->ilist; numneigh = list->numneigh; firstneigh = list->firstneigh; @@ -636,7 +644,8 @@ void FixQtpieReaxFF::pre_force(int /*vflag*/) if (n > n_cap*DANGER_ZONE || m_fill > m_cap*DANGER_ZONE) reallocate_matrix(); - if (efield) get_chi_field(); + // if (efield) get_chi_field(); + calc_chi_eff(); init_matvec(); @@ -676,8 +685,8 @@ void FixQtpieReaxFF::init_matvec() /* init pre-conditioner for H and init solution vectors */ Hdia_inv[i] = 1. / eta[atom->type[i]]; - b_s[i] = -chi[atom->type[i]]; - if (efield) b_s[i] -= chi_field[i]; + b_s[i] = -chi_eff[i]; + // if (efield) b_s[i] -= chi_field[i]; b_t[i] = -1.0; /* quadratic extrapolation for s & t from previous solutions */ @@ -1195,3 +1204,122 @@ void FixQtpieReaxFF::get_chi_field() } } } + +void FixQtpieReaxFF::calc_chi_eff() +{ + memset(&chi_eff[0],0,atom->nmax*sizeof(double)); + const int KSCREEN = 10; + const double ZERONAME = 1.0e-50; + const double ANG_TO_BOHRRAD = 1.8897259886; // 1 Ang = 1.8897259886 Bohr radius + + double R,a_min,OvIntMaxR,voltage,overlap,nominator,denominator; + double ea,eb,chia,chib,p,m; + double phia,phib; + int i,j; + const int ntypes = atom->ntypes; + const int *type = atom->type; + // const auto x = (const double * const *)atom->x; + double **x = atom->x; + + // Use integral pre-screening for overlap calculations + a_min = find_min(gauss_exp,ntypes+1); + OvIntMaxR = sqrt(pow(a_min,-1.)*log(pow(M_PI/(2.*a_min),3.)*pow(10.,2.*KSCREEN))); + + ghost_cutoff = MAX(neighbor->cutneighmax,comm->cutghostuser); + if(ghost_cutoff < OvIntMaxR/ANG_TO_BOHRRAD) { + error->all(FLERR,"ghost cutoff error"); + // char errmsg[256]; + // snprintf(errmsg, 256,"qtpie/reaxff: limit distance for overlap integral: %f " + // "Angstrom > ghost cutoff: %f Angstrom. Increase the ghost atom cutoff " + // "with comm_modify.",OvIntMaxR/ANG_TO_BOHRRAD,ghost_cutoff); + // // system->error_ptr->all(FLERR,errmsg); + // error->all(FLERR,errmsg); + } + + for (i = 0; i < nn; i++) { + + // type_i = type[i]; + ea = gauss_exp[type[i]]; + chia = chi[type[i]]; + + nominator = denominator = 0.0; + + for (j = 0; j < ng; j++) { + + R = distance(x[i],x[j])*ANG_TO_BOHRRAD; + overlap = voltage = 0.0; + + if (R < OvIntMaxR) + { + // type_j = type[j]; + eb = gauss_exp[type[j]]; + chib = chi[type[j]]; + + // The expressions below are in atomic units + // Implementation from Chen Jiahao, Theory and applications of fluctuating-charge models, 2009 (with normalization constants added) + p = ea + eb; + m = ea * eb / p; + overlap = pow((4. * m / p), 0.75) * exp(-m * R * R); + + // Implementation from T. Halgaker et al., Molecular electronic-structure theory, 2000 +// p = ea + eb; +// m = ea * eb / p; +// Overlap = pow((M_PI / p), 1.5) * exp(-m * R * R); + + if (efield) { + phib = efield_potential(x[j]); + voltage = chia - chib + phib; + } else { + voltage = chia - chib; + } + nominator += voltage * overlap; + denominator += overlap; + } + } + if (denominator != 0.0 && nominator != 0.0) + chi_eff[i] = nominator / denominator; + else + chi_eff[i] = ZERONAME; + + if (efield) { + phia = efield_potential(x[i]); + chi_eff[i] -= phia; + } + } +} + +/* ---------------------------------------------------------------------- */ + +double FixQtpieReaxFF::find_min(double *array, int array_length) +{ + // since types start from 1, gaussian exponents start from 1 + double smallest = array[1]; + for (int i = 1; i < array_length; i++) + { + if (array[i] < smallest) + smallest = array[i]; + } + return smallest; +} + +/* ---------------------------------------------------------------------- */ + +double FixQtpieReaxFF::distance(double *loc1, double *loc2) +{ + double distx, disty, distz; + distx = loc2[0] - loc1[0]; + disty = loc2[1] - loc1[1]; + distz = loc2[2] - loc1[2]; + return sqrt(distx*distx + disty*disty + distz*distz); +} + +/* ---------------------------------------------------------------------- */ + +double FixQtpieReaxFF::efield_potential(double *x) +{ + double x_efcomp, y_efcomp, z_efcomp; + x_efcomp = x[0] * efield->ex; + y_efcomp = x[1] * efield->ey; + z_efcomp = x[2] * efield->ez; + return x_efcomp + y_efcomp + z_efcomp; +} diff --git a/src/REAXFF/fix_qtpie_reaxff.h b/src/REAXFF/fix_qtpie_reaxff.h index fcdedddf26..216fc56468 100644 --- a/src/REAXFF/fix_qtpie_reaxff.h +++ b/src/REAXFF/fix_qtpie_reaxff.h @@ -48,7 +48,7 @@ class FixQtpieReaxFF : public Fix { protected: int nevery, reaxflag; int matvecs; - int nn, m_fill; + int nn, ng, m_fill; int n_cap, nmax, m_cap; int pack_flag; int nlevels_respa; @@ -90,16 +90,11 @@ class FixQtpieReaxFF : public Fix { char *pertype_option; // argument to determine how per-type info is obtained - // Params from Kritikos - could rename or move to protected later char *gauss_file; // input file for gaussian exponents for each type of REAXFF file - double cutghost; // ghost atoms cutoff (used for check) - int nn_prev; // number of local atoms; needed for memory reallocation of chi_eff (when multiprocessing) + double ghost_cutoff; // ghost atoms cutoff double *gauss_exp; // array of gaussian exponents double *chi_eff; // array of effective electronegativities - double *chi_eff_init; // array of effective electronegativities for FixQEqReax::init_storage() - // void calculate_chi_eff(LAMMPS_NS::Atom *atom, reax_system *system, double *chi, - // int ni, int nj, double *lchi_eff); virtual void pertype_parameters(char *); void init_shielding(); void init_taper(); @@ -110,6 +105,11 @@ class FixQtpieReaxFF : public Fix { virtual void deallocate_matrix(); void reallocate_matrix(); + void calc_chi_eff(); + double find_min(double*, int); + double distance(double*, double*); + double efield_potential(double*); + virtual void init_matvec(); void init_H(); virtual void compute_H(); From 8c8882927c0a2db1c35d9b22d5b74ee942886172 Mon Sep 17 00:00:00 2001 From: Navraj Lalli Date: Thu, 15 Aug 2024 16:20:45 +0100 Subject: [PATCH 04/40] Rename variables in calc_chi_eff() --- src/REAXFF/fix_qtpie_reaxff.cpp | 43 +++++++++++++++++++++++---------- 1 file changed, 30 insertions(+), 13 deletions(-) diff --git a/src/REAXFF/fix_qtpie_reaxff.cpp b/src/REAXFF/fix_qtpie_reaxff.cpp index 9f3714c32d..f3151a85ba 100644 --- a/src/REAXFF/fix_qtpie_reaxff.cpp +++ b/src/REAXFF/fix_qtpie_reaxff.cpp @@ -13,7 +13,7 @@ ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- - Contributing author: Navraj S Lalli & Efstratios Kritikos (Imperial College London) + Contributing author: Navraj S Lalli & Efstratios M Kritikos (Imperial College London) ------------------------------------------------------------------------- */ #include "fix_qtpie_reaxff.h" @@ -1209,10 +1209,10 @@ void FixQtpieReaxFF::calc_chi_eff() { memset(&chi_eff[0],0,atom->nmax*sizeof(double)); const int KSCREEN = 10; - const double ZERONAME = 1.0e-50; + // const double ZERONAME = 1.0e-50; const double ANG_TO_BOHRRAD = 1.8897259886; // 1 Ang = 1.8897259886 Bohr radius - double R,a_min,OvIntMaxR,voltage,overlap,nominator,denominator; + double R,a_min,OvIntMaxR,overlap,sum_n,sum_d; double ea,eb,chia,chib,p,m; double phia,phib; int i,j; @@ -1224,8 +1224,15 @@ void FixQtpieReaxFF::calc_chi_eff() // Use integral pre-screening for overlap calculations a_min = find_min(gauss_exp,ntypes+1); OvIntMaxR = sqrt(pow(a_min,-1.)*log(pow(M_PI/(2.*a_min),3.)*pow(10.,2.*KSCREEN))); + // OvIntMaxR = 0.5; + // if (comm->me == 0) { + // std::cout << "OvIntMaxR is " << OvIntMaxR/ANG_TO_BOHRRAD; + // } ghost_cutoff = MAX(neighbor->cutneighmax,comm->cutghostuser); + // if (comm->me == 0) { + // std::cout << "ghost_cutoff is " << ghost_cutoff; + // } if(ghost_cutoff < OvIntMaxR/ANG_TO_BOHRRAD) { error->all(FLERR,"ghost cutoff error"); // char errmsg[256]; @@ -1242,12 +1249,14 @@ void FixQtpieReaxFF::calc_chi_eff() ea = gauss_exp[type[i]]; chia = chi[type[i]]; - nominator = denominator = 0.0; + // nominator = denominator = 0.0; + sum_n = 0.0; + sum_d = 0.0; for (j = 0; j < ng; j++) { - R = distance(x[i],x[j])*ANG_TO_BOHRRAD; - overlap = voltage = 0.0; + R = distance(x[i],x[j])*ANG_TO_BOHRRAD; // Distance between atoms as a multiple of Bohr radius + // overlap = voltage = 0.0; if (R < OvIntMaxR) { @@ -1268,18 +1277,26 @@ void FixQtpieReaxFF::calc_chi_eff() if (efield) { phib = efield_potential(x[j]); - voltage = chia - chib + phib; + sum_n += (chia - chib + phib) * overlap; + // voltage = chia - chib + phib; } else { - voltage = chia - chib; + sum_n += (chia - chib) * overlap; + // voltage = chia - chib; } - nominator += voltage * overlap; - denominator += overlap; + sum_d += overlap; + // nominator += voltage * overlap; + // denominator += overlap; } } - if (denominator != 0.0 && nominator != 0.0) - chi_eff[i] = nominator / denominator; + + if (fabs(sum_n) < SMALL && fabs(sum_d) < SMALL) + chi_eff[i] = 0.0; // SMALL; else - chi_eff[i] = ZERONAME; + chi_eff[i] = sum_n / sum_d; + // if (denominator != 0.0 && nominator != 0.0) + // chi_eff[i] = nominator / denominator; + // else + // chi_eff[i] = ZERONAME; if (efield) { phia = efield_potential(x[i]); From 6d47e417415d3dc4cb6f28d909e81be8972e8a9c Mon Sep 17 00:00:00 2001 From: Navraj Lalli Date: Tue, 20 Aug 2024 16:00:51 +0100 Subject: [PATCH 05/40] Update calculation of chi_eff chi_eff can now be calculated when atom or equal style variables are used for the electric field, under the restriction that the electric field is applied to all atoms. --- src/REAXFF/fix_qtpie_reaxff.cpp | 156 ++++++++++++++------------------ src/REAXFF/fix_qtpie_reaxff.h | 4 +- 2 files changed, 71 insertions(+), 89 deletions(-) diff --git a/src/REAXFF/fix_qtpie_reaxff.cpp b/src/REAXFF/fix_qtpie_reaxff.cpp index f3151a85ba..6c8e15d965 100644 --- a/src/REAXFF/fix_qtpie_reaxff.cpp +++ b/src/REAXFF/fix_qtpie_reaxff.cpp @@ -223,7 +223,7 @@ void FixQtpieReaxFF::pertype_parameters(char *arg) reader.ignore_comments = false; for (int i = 1; i <= ntypes; i++) { const char *line = reader.next_line(); - std::cout << line; + std::cout << "Orbital exponent " << line; if (!line) throw TokenizerException("Fix qtpie/reaxff: Invalid param file format",""); ValueTokenizer values(line); @@ -438,7 +438,6 @@ void FixQtpieReaxFF::init() error->all(FLERR,"Fix {} group has no atoms", style); // compute net charge and print warning if too large - double qsum_local = 0.0, qsum = 0.0; for (int i = 0; i < atom->nlocal; i++) { if (atom->mask[i] & groupbit) @@ -450,7 +449,6 @@ void FixQtpieReaxFF::init() error->warning(FLERR,"Fix {} group is not charge neutral, net charge = {:.8}", style, qsum); // get pointer to fix efield if present. there may be at most one instance of fix efield in use. - efield = nullptr; auto fixes = modify->get_fix_by_style("^efield"); if (fixes.size() == 1) efield = dynamic_cast(fixes.front()); @@ -463,19 +461,27 @@ void FixQtpieReaxFF::init() if (strcmp(update->unit_style,"real") != 0) error->all(FLERR,"Must use unit_style real with fix {} and external fields", style); + if (efield->groupbit != 1){ // if efield is not applied to all atoms + error->all(FLERR,"Must use group id all for fix efield when using fix {}", style); + } + + if (efield->region){ // if efield is not applied to all atoms + error->all(FLERR,"Keyword region not supported for fix efield when using fix {}", style); + } + if (efield->varflag == FixEfield::ATOM && efield->pstyle != FixEfield::ATOM) error->all(FLERR,"Atom-style external electric field requires atom-style " "potential variable when used with fix {}", style); - if (((efield->xstyle != FixEfield::CONSTANT) && domain->xperiodic) || - ((efield->ystyle != FixEfield::CONSTANT) && domain->yperiodic) || - ((efield->zstyle != FixEfield::CONSTANT) && domain->zperiodic)) - error->all(FLERR,"Must not have electric field component in direction of periodic " - "boundary when using charge equilibration with ReaxFF."); - if (((fabs(efield->ex) > SMALL) && domain->xperiodic) || - ((fabs(efield->ey) > SMALL) && domain->yperiodic) || - ((fabs(efield->ez) > SMALL) && domain->zperiodic)) - error->all(FLERR,"Must not have electric field component in direction of periodic " - "boundary when using charge equilibration with ReaxFF."); + // if (((efield->xstyle != FixEfield::CONSTANT) && domain->xperiodic) || + // ((efield->ystyle != FixEfield::CONSTANT) && domain->yperiodic) || + // ((efield->zstyle != FixEfield::CONSTANT) && domain->zperiodic)) + // error->all(FLERR,"Must not have electric field component in direction of periodic " + // "boundary when using charge equilibration with ReaxFF."); + // if (((fabs(efield->ex) > SMALL) && domain->xperiodic) || + // ((fabs(efield->ey) > SMALL) && domain->yperiodic) || + // ((fabs(efield->ez) > SMALL) && domain->zperiodic)) + // error->all(FLERR,"Must not have electric field component in direction of periodic " + // "boundary when using charge equilibration with ReaxFF."); } // we need a half neighbor list w/ Newton off @@ -1208,59 +1214,55 @@ void FixQtpieReaxFF::get_chi_field() void FixQtpieReaxFF::calc_chi_eff() { memset(&chi_eff[0],0,atom->nmax*sizeof(double)); - const int KSCREEN = 10; - // const double ZERONAME = 1.0e-50; - const double ANG_TO_BOHRRAD = 1.8897259886; // 1 Ang = 1.8897259886 Bohr radius - double R,a_min,OvIntMaxR,overlap,sum_n,sum_d; - double ea,eb,chia,chib,p,m; - double phia,phib; - int i,j; + const int KSCREEN = 10; + const double ANG_TO_BOHRRAD = 1.8897259886; // 1 Ang = 1.8897259886 Bohr radius + const auto x = (const double * const *)atom->x; const int ntypes = atom->ntypes; const int *type = atom->type; - // const auto x = (const double * const *)atom->x; - double **x = atom->x; - // Use integral pre-screening for overlap calculations - a_min = find_min(gauss_exp,ntypes+1); - OvIntMaxR = sqrt(pow(a_min,-1.)*log(pow(M_PI/(2.*a_min),3.)*pow(10.,2.*KSCREEN))); - // OvIntMaxR = 0.5; - // if (comm->me == 0) { - // std::cout << "OvIntMaxR is " << OvIntMaxR/ANG_TO_BOHRRAD; - // } + double dist,overlap,sum_n,sum_d,ea,eb,chia,chib,phia,phib,p,m; + int i,j; - ghost_cutoff = MAX(neighbor->cutneighmax,comm->cutghostuser); - // if (comm->me == 0) { - // std::cout << "ghost_cutoff is " << ghost_cutoff; - // } - if(ghost_cutoff < OvIntMaxR/ANG_TO_BOHRRAD) { - error->all(FLERR,"ghost cutoff error"); - // char errmsg[256]; - // snprintf(errmsg, 256,"qtpie/reaxff: limit distance for overlap integral: %f " - // "Angstrom > ghost cutoff: %f Angstrom. Increase the ghost atom cutoff " - // "with comm_modify.",OvIntMaxR/ANG_TO_BOHRRAD,ghost_cutoff); - // // system->error_ptr->all(FLERR,errmsg); - // error->all(FLERR,errmsg); + // efield energy is in real units of kcal/mol, factor needed for conversion to eV + const double qe2f = force->qe2f; + const double factor = 1.0/qe2f; + + if (efield) { + if (efield->varflag != FixEfield::CONSTANT) + efield->update_efield_variables(); } - for (i = 0; i < nn; i++) { + // use integral pre-screening for overlap calculations + const double emin = find_min(gauss_exp,ntypes+1); + const double dist_cutoff = sqrt(pow(emin,-1.)*log(pow(M_PI/(2.*emin),3.)*pow(10.,2.*KSCREEN))); - // type_i = type[i]; + const double comm_cutoff = MAX(neighbor->cutneighmax,comm->cutghostuser); + if(comm_cutoff < dist_cutoff/ANG_TO_BOHRRAD) { + error->all(FLERR,"comm cutoff = {} Angstrom is smaller than distance cutoff = {} Angstrom " + "for overlap integral in {}. Increase comm cutoff with comm_modify", + comm_cutoff, dist_cutoff/ANG_TO_BOHRRAD, style); + } + + // compute chi_eff for each local atom + for (i = 0; i < nn; i++) { ea = gauss_exp[type[i]]; chia = chi[type[i]]; + if (efield) { + if (efield->varflag != FixEfield::ATOM) { + phia = factor*(x[i][0]*efield->ex + x[i][1]*efield->ey + x[i][2]*efield->ez); + } else { // atom-style potential from FixEfield + phia = efield->efield[i][3]; + } + } - // nominator = denominator = 0.0; sum_n = 0.0; sum_d = 0.0; for (j = 0; j < ng; j++) { + dist = distance(x[i],x[j])*ANG_TO_BOHRRAD; // distance between atoms as a multiple of Bohr radius - R = distance(x[i],x[j])*ANG_TO_BOHRRAD; // Distance between atoms as a multiple of Bohr radius - // overlap = voltage = 0.0; - - if (R < OvIntMaxR) - { - // type_j = type[j]; + if (dist < dist_cutoff) { eb = gauss_exp[type[j]]; chib = chi[type[j]]; @@ -1268,7 +1270,7 @@ void FixQtpieReaxFF::calc_chi_eff() // Implementation from Chen Jiahao, Theory and applications of fluctuating-charge models, 2009 (with normalization constants added) p = ea + eb; m = ea * eb / p; - overlap = pow((4. * m / p), 0.75) * exp(-m * R * R); + overlap = pow((4.0*m/p),0.75) * exp(-m*dist*dist); // Implementation from T. Halgaker et al., Molecular electronic-structure theory, 2000 // p = ea + eb; @@ -1276,32 +1278,25 @@ void FixQtpieReaxFF::calc_chi_eff() // Overlap = pow((M_PI / p), 1.5) * exp(-m * R * R); if (efield) { - phib = efield_potential(x[j]); - sum_n += (chia - chib + phib) * overlap; - // voltage = chia - chib + phib; + 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 - chib + phib - phia) * overlap; } else { sum_n += (chia - chib) * overlap; - // voltage = chia - chib; } sum_d += overlap; - // nominator += voltage * overlap; - // denominator += overlap; } } - if (fabs(sum_n) < SMALL && fabs(sum_d) < SMALL) - chi_eff[i] = 0.0; // SMALL; - else - chi_eff[i] = sum_n / sum_d; - // if (denominator != 0.0 && nominator != 0.0) - // chi_eff[i] = nominator / denominator; - // else - // chi_eff[i] = ZERONAME; + chi_eff[i] = sum_n / sum_d; - if (efield) { - phia = efield_potential(x[i]); - chi_eff[i] -= phia; - } + if (fabs(sum_n) < SMALL && fabs(sum_d) < SMALL) + error->all(FLERR,"Unexpected value: fabs(sum_d) is {}", fabs(sum_d)); + if (fabs(sum_d) < 1.0) + error->all(FLERR,"Unexpected value: fabs(sum_d) is {}", fabs(sum_d)); } } @@ -1321,22 +1316,11 @@ double FixQtpieReaxFF::find_min(double *array, int array_length) /* ---------------------------------------------------------------------- */ -double FixQtpieReaxFF::distance(double *loc1, double *loc2) +double FixQtpieReaxFF::distance(const double *posa, const double *posb) { - double distx, disty, distz; - distx = loc2[0] - loc1[0]; - disty = loc2[1] - loc1[1]; - distz = loc2[2] - loc1[2]; - return sqrt(distx*distx + disty*disty + distz*distz); -} - -/* ---------------------------------------------------------------------- */ - -double FixQtpieReaxFF::efield_potential(double *x) -{ - double x_efcomp, y_efcomp, z_efcomp; - x_efcomp = x[0] * efield->ex; - y_efcomp = x[1] * efield->ey; - z_efcomp = x[2] * efield->ez; - return x_efcomp + y_efcomp + z_efcomp; + 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); } diff --git a/src/REAXFF/fix_qtpie_reaxff.h b/src/REAXFF/fix_qtpie_reaxff.h index 216fc56468..7b6feacabb 100644 --- a/src/REAXFF/fix_qtpie_reaxff.h +++ b/src/REAXFF/fix_qtpie_reaxff.h @@ -91,7 +91,6 @@ class FixQtpieReaxFF : public Fix { char *pertype_option; // argument to determine how per-type info is obtained char *gauss_file; // input file for gaussian exponents for each type of REAXFF file - double ghost_cutoff; // ghost atoms cutoff double *gauss_exp; // array of gaussian exponents double *chi_eff; // array of effective electronegativities @@ -107,8 +106,7 @@ class FixQtpieReaxFF : public Fix { void calc_chi_eff(); double find_min(double*, int); - double distance(double*, double*); - double efield_potential(double*); + double distance(const double*, const double*); virtual void init_matvec(); void init_H(); From 27e911cd10fb976a16d39b80980de3755cadb080 Mon Sep 17 00:00:00 2001 From: Navraj Lalli Date: Tue, 20 Aug 2024 17:36:56 +0100 Subject: [PATCH 06/40] Remove chi_field --- src/REAXFF/fix_qtpie_reaxff.cpp | 60 --------------------------------- src/REAXFF/fix_qtpie_reaxff.h | 15 ++++----- 2 files changed, 6 insertions(+), 69 deletions(-) diff --git a/src/REAXFF/fix_qtpie_reaxff.cpp b/src/REAXFF/fix_qtpie_reaxff.cpp index 6c8e15d965..82b0f8b595 100644 --- a/src/REAXFF/fix_qtpie_reaxff.cpp +++ b/src/REAXFF/fix_qtpie_reaxff.cpp @@ -114,7 +114,6 @@ FixQtpieReaxFF::FixQtpieReaxFF(LAMMPS *lmp, int narg, char **arg) : Hdia_inv = nullptr; b_s = nullptr; - chi_field = nullptr; chi_eff = nullptr; b_t = nullptr; b_prc = nullptr; @@ -327,7 +326,6 @@ void FixQtpieReaxFF::allocate_storage() memory->create(Hdia_inv,nmax,"qeq:Hdia_inv"); memory->create(b_s,nmax,"qeq:b_s"); - memory->create(chi_field,nmax,"qeq:chi_field"); memory->create(chi_eff,nmax,"qtpie:chi_eff"); memory->create(b_t,nmax,"qeq:b_t"); memory->create(b_prc,nmax,"qeq:b_prc"); @@ -355,7 +353,6 @@ void FixQtpieReaxFF::deallocate_storage() memory->destroy(b_t); memory->destroy(b_prc); memory->destroy(b_prm); - memory->destroy(chi_field); memory->destroy(chi_eff); memory->destroy(p); @@ -604,7 +601,6 @@ void FixQtpieReaxFF::min_setup_pre_force(int vflag) void FixQtpieReaxFF::init_storage() { - // if (efield) get_chi_field(); calc_chi_eff(); for (int ii = 0; ii < nn; ii++) { @@ -612,7 +608,6 @@ void FixQtpieReaxFF::init_storage() if (atom->mask[i] & groupbit) { Hdia_inv[i] = 1. / eta[atom->type[i]]; b_s[i] = -chi_eff[i]; - // if (efield) b_s[i] -= chi_field[i]; b_t[i] = -1.0; b_prc[i] = 0; b_prm[i] = 0; @@ -650,7 +645,6 @@ void FixQtpieReaxFF::pre_force(int /*vflag*/) if (n > n_cap*DANGER_ZONE || m_fill > m_cap*DANGER_ZONE) reallocate_matrix(); - // if (efield) get_chi_field(); calc_chi_eff(); init_matvec(); @@ -692,7 +686,6 @@ void FixQtpieReaxFF::init_matvec() /* init pre-conditioner for H and init solution vectors */ Hdia_inv[i] = 1. / eta[atom->type[i]]; b_s[i] = -chi_eff[i]; - // if (efield) b_s[i] -= chi_field[i]; b_t[i] = -1.0; /* quadratic extrapolation for s & t from previous solutions */ @@ -1158,59 +1151,6 @@ void FixQtpieReaxFF::vector_add(double* dest, double c, double* v, int k) /* ---------------------------------------------------------------------- */ -void FixQtpieReaxFF::get_chi_field() -{ - memset(&chi_field[0],0,atom->nmax*sizeof(double)); - if (!efield) return; - - const auto x = (const double * const *)atom->x; - const int *mask = atom->mask; - const imageint *image = atom->image; - const int nlocal = atom->nlocal; - - - // update electric field region if necessary - - Region *region = efield->region; - if (region) region->prematch(); - - // efield energy is in real units of kcal/mol/angstrom, need to convert to eV - - const double qe2f = force->qe2f; - const double factor = -1.0/qe2f; - - - if (efield->varflag != FixEfield::CONSTANT) - efield->update_efield_variables(); - - // atom selection is for the group of fix efield - - double unwrap[3]; - const double ex = efield->ex; - const double ey = efield->ey; - const double ez = efield->ez; - const int efgroupbit = efield->groupbit; - - // charge interactions - // force = qE, potential energy = F dot x in unwrapped coords - if (efield->varflag != FixEfield::ATOM) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & efgroupbit) { - if (region && !region->match(x[i][0],x[i][1],x[i][2])) continue; - domain->unmap(x[i],image[i],unwrap); - chi_field[i] = factor*(ex*unwrap[0] + ey*unwrap[1] + ez*unwrap[2]); - } - } - } else { // must use atom-style potential from FixEfield - for (int i = 0; i < nlocal; i++) { - if (mask[i] & efgroupbit) { - if (region && !region->match(x[i][0],x[i][1],x[i][2])) continue; - chi_field[i] = -efield->efield[i][3]; - } - } - } -} - void FixQtpieReaxFF::calc_chi_eff() { memset(&chi_eff[0],0,atom->nmax*sizeof(double)); diff --git a/src/REAXFF/fix_qtpie_reaxff.h b/src/REAXFF/fix_qtpie_reaxff.h index 7b6feacabb..07afdf9439 100644 --- a/src/REAXFF/fix_qtpie_reaxff.h +++ b/src/REAXFF/fix_qtpie_reaxff.h @@ -82,17 +82,15 @@ class FixQtpieReaxFF : public Fix { double *Hdia_inv; double *b_s, *b_t; double *b_prc, *b_prm; - double *chi_field; + double *chi_eff; // array of effective electronegativities //CG storage double *p, *q, *r, *d; int imax, maxwarn; char *pertype_option; // argument to determine how per-type info is obtained - - char *gauss_file; // input file for gaussian exponents for each type of REAXFF file + char *gauss_file; // input file for gaussian exponents double *gauss_exp; // array of gaussian exponents - double *chi_eff; // array of effective electronegativities virtual void pertype_parameters(char *); void init_shielding(); @@ -104,10 +102,6 @@ class FixQtpieReaxFF : public Fix { virtual void deallocate_matrix(); void reallocate_matrix(); - void calc_chi_eff(); - double find_min(double*, int); - double distance(const double*, const double*); - virtual void init_matvec(); void init_H(); virtual void compute_H(); @@ -134,7 +128,10 @@ class FixQtpieReaxFF : public Fix { virtual void vector_sum(double *, double, double *, double, double *, int); virtual void vector_add(double *, double, double *, int); - virtual void get_chi_field(); + void calc_chi_eff(); + double find_min(double*, int); + double distance(const double*, const double*); + // dual CG support int dual_enabled; // 0: Original, separate s & t optimization; 1: dual optimization From 5021c8c971d8f60acd53109a221bae0e863106b2 Mon Sep 17 00:00:00 2001 From: Navraj Lalli Date: Tue, 20 Aug 2024 17:39:44 +0100 Subject: [PATCH 07/40] Replace qeq with qtpie --- src/REAXFF/fix_qtpie_reaxff.cpp | 60 ++++++++++++++++----------------- 1 file changed, 30 insertions(+), 30 deletions(-) diff --git a/src/REAXFF/fix_qtpie_reaxff.cpp b/src/REAXFF/fix_qtpie_reaxff.cpp index 82b0f8b595..aecb4352fa 100644 --- a/src/REAXFF/fix_qtpie_reaxff.cpp +++ b/src/REAXFF/fix_qtpie_reaxff.cpp @@ -278,27 +278,27 @@ void FixQtpieReaxFF::pertype_parameters(char *arg) reaxflag = 0; - memory->create(chi,ntypes+1,"qeq/reaxff:chi"); - memory->create(eta,ntypes+1,"qeq/reaxff:eta"); - memory->create(gamma,ntypes+1,"qeq/reaxff:gamma"); + memory->create(chi,ntypes+1,"qtpie/reaxff:chi"); + memory->create(eta,ntypes+1,"qtpie/reaxff:eta"); + memory->create(gamma,ntypes+1,"qtpie/reaxff:gamma"); if (comm->me == 0) { chi[0] = eta[0] = gamma[0] = 0.0; try { - TextFileReader reader(arg,"qeq/reaxff parameter"); + TextFileReader reader(arg,"qtpie/reaxff parameter"); reader.ignore_comments = false; for (int i = 1; i <= ntypes; i++) { const char *line = reader.next_line(); if (!line) - throw TokenizerException("Fix qeq/reaxff: Invalid param file format",""); + throw TokenizerException("Fix qtpie/reaxff: Invalid param file format",""); ValueTokenizer values(line); if (values.count() != 4) - throw TokenizerException("Fix qeq/reaxff: Incorrect format of param file",""); + throw TokenizerException("Fix qtpie/reaxff: Incorrect format of param file",""); int itype = values.next_int(); if ((itype < 1) || (itype > ntypes)) - throw TokenizerException("Fix qeq/reaxff: Invalid atom type in param file", + throw TokenizerException("Fix qtpie/reaxff: Invalid atom type in param file", std::to_string(itype)); chi[itype] = values.next_double(); @@ -321,24 +321,24 @@ void FixQtpieReaxFF::allocate_storage() { nmax = atom->nmax; - memory->create(s,nmax,"qeq:s"); - memory->create(t,nmax,"qeq:t"); + memory->create(s,nmax,"qtpie:s"); + memory->create(t,nmax,"qtpie:t"); - memory->create(Hdia_inv,nmax,"qeq:Hdia_inv"); - memory->create(b_s,nmax,"qeq:b_s"); + memory->create(Hdia_inv,nmax,"qtpie:Hdia_inv"); + memory->create(b_s,nmax,"qtpie:b_s"); memory->create(chi_eff,nmax,"qtpie:chi_eff"); - memory->create(b_t,nmax,"qeq:b_t"); - memory->create(b_prc,nmax,"qeq:b_prc"); - memory->create(b_prm,nmax,"qeq:b_prm"); + memory->create(b_t,nmax,"qtpie:b_t"); + memory->create(b_prc,nmax,"qtpie:b_prc"); + memory->create(b_prm,nmax,"qtpie:b_prm"); // dual CG support int size = nmax; if (dual_enabled) size*= 2; - memory->create(p,size,"qeq:p"); - memory->create(q,size,"qeq:q"); - memory->create(r,size,"qeq:r"); - memory->create(d,size,"qeq:d"); + memory->create(p,size,"qtpie:p"); + memory->create(q,size,"qtpie:q"); + memory->create(r,size,"qtpie:r"); + memory->create(d,size,"qtpie:d"); } /* ---------------------------------------------------------------------- */ @@ -400,10 +400,10 @@ void FixQtpieReaxFF::allocate_matrix() H.n = n_cap; H.m = m_cap; - memory->create(H.firstnbr,n_cap,"qeq:H.firstnbr"); - memory->create(H.numnbrs,n_cap,"qeq:H.numnbrs"); - memory->create(H.jlist,m_cap,"qeq:H.jlist"); - memory->create(H.val,m_cap,"qeq:H.val"); + memory->create(H.firstnbr,n_cap,"qtpie:H.firstnbr"); + memory->create(H.numnbrs,n_cap,"qtpie:H.numnbrs"); + memory->create(H.jlist,m_cap,"qtpie:H.jlist"); + memory->create(H.val,m_cap,"qtpie:H.val"); } /* ---------------------------------------------------------------------- */ @@ -516,7 +516,7 @@ void FixQtpieReaxFF::init_shielding() ntypes = atom->ntypes; if (shld == nullptr) - memory->create(shld,ntypes+1,ntypes+1,"qeq:shielding"); + memory->create(shld,ntypes+1,ntypes+1,"qtpie:shielding"); for (i = 1; i <= ntypes; ++i) for (j = 1; j <= ntypes; ++j) @@ -530,11 +530,11 @@ void FixQtpieReaxFF::init_taper() double d7, swa2, swa3, swb2, swb3; if (fabs(swa) > 0.01 && comm->me == 0) - error->warning(FLERR,"Fix qeq/reaxff has non-zero lower Taper radius cutoff"); + error->warning(FLERR,"Fix qtpie/reaxff has non-zero lower Taper radius cutoff"); if (swb < 0) - error->all(FLERR, "Fix qeq/reaxff has negative upper Taper radius cutoff"); + error->all(FLERR, "Fix qtpie/reaxff has negative upper Taper radius cutoff"); else if (swb < 5 && comm->me == 0) - error->warning(FLERR,"Fix qeq/reaxff has very low Taper radius cutoff"); + error->warning(FLERR,"Fix qtpie/reaxff has very low Taper radius cutoff"); d7 = pow(swb - swa, 7); swa2 = SQR(swa); @@ -760,7 +760,7 @@ void FixQtpieReaxFF::compute_H() } if (m_fill >= H.m) - error->all(FLERR,"Fix qeq/reaxff H matrix size has been exceeded: m_fill={} H.m={}\n", + error->all(FLERR,"Fix qtpie/reaxff H matrix size has been exceeded: m_fill={} H.m={}\n", m_fill, H.m); } @@ -835,7 +835,7 @@ int FixQtpieReaxFF::CG(double *b, double *x) } if ((i >= imax) && maxwarn && (comm->me == 0)) - error->warning(FLERR, "Fix qeq/reaxff CG convergence failed after {} iterations at step {}", + error->warning(FLERR, "Fix qtpie/reaxff CG convergence failed after {} iterations at step {}", i,update->ntimestep); return i; } @@ -1018,8 +1018,8 @@ double FixQtpieReaxFF::memory_usage() void FixQtpieReaxFF::grow_arrays(int nmax) { - memory->grow(s_hist,nmax,nprev,"qeq:s_hist"); - memory->grow(t_hist,nmax,nprev,"qeq:t_hist"); + memory->grow(s_hist,nmax,nprev,"qtpie:s_hist"); + memory->grow(t_hist,nmax,nprev,"qtpie:t_hist"); } /* ---------------------------------------------------------------------- From 149d9b310d0b144065c9525fa5444f3bad804090 Mon Sep 17 00:00:00 2001 From: Navraj Lalli Date: Tue, 20 Aug 2024 18:03:50 +0100 Subject: [PATCH 08/40] Remove dual as a possible keyword argument --- src/REAXFF/fix_qtpie_reaxff.cpp | 23 +++-------------------- src/REAXFF/fix_qtpie_reaxff.h | 3 --- 2 files changed, 3 insertions(+), 23 deletions(-) diff --git a/src/REAXFF/fix_qtpie_reaxff.cpp b/src/REAXFF/fix_qtpie_reaxff.cpp index aecb4352fa..4066103152 100644 --- a/src/REAXFF/fix_qtpie_reaxff.cpp +++ b/src/REAXFF/fix_qtpie_reaxff.cpp @@ -74,7 +74,7 @@ FixQtpieReaxFF::FixQtpieReaxFF(LAMMPS *lmp, int narg, char **arg) : imax = 200; maxwarn = 1; - if ((narg < 9) || (narg > 13)) error->all(FLERR,"Illegal fix {} command", style); + if ((narg < 9) || (narg > 12)) error->all(FLERR,"Illegal fix {} command", style); nevery = utils::inumeric(FLERR,arg[3],false,lmp); if (nevery <= 0) error->all(FLERR,"Illegal fix {} command", style); @@ -85,15 +85,9 @@ FixQtpieReaxFF::FixQtpieReaxFF(LAMMPS *lmp, int narg, char **arg) : pertype_option = utils::strdup(arg[7]); gauss_file = utils::strdup(arg[8]); - // dual CG support only available for OPENMP variant - // check for compatibility is in Fix::post_constructor() - - dual_enabled = 0; - int iarg = 9; while (iarg < narg) { - if (strcmp(arg[iarg],"dual") == 0) dual_enabled = 1; - else if (strcmp(arg[iarg],"nowarn") == 0) maxwarn = 0; + if (strcmp(arg[iarg],"nowarn") == 0) maxwarn = 0; else if (strcmp(arg[iarg],"maxiter") == 0) { if (iarg+1 > narg-1) error->all(FLERR,"Illegal fix {} command", style); @@ -133,11 +127,7 @@ FixQtpieReaxFF::FixQtpieReaxFF(LAMMPS *lmp, int narg, char **arg) : H.jlist = nullptr; H.val = nullptr; - // dual CG support - // Update comm sizes for this fix - - if (dual_enabled) comm_forward = comm_reverse = 2; - else comm_forward = comm_reverse = 1; + comm_forward = comm_reverse = 1; // perform initial allocation of atom-based arrays // register with Atom class @@ -189,8 +179,6 @@ void FixQtpieReaxFF::post_constructor() s_hist[i][j] = t_hist[i][j] = 0; pertype_parameters(pertype_option); - if (dual_enabled) - error->all(FLERR,"Dual keyword only supported with fix qeq/reaxff/omp"); } /* ---------------------------------------------------------------------- */ @@ -331,9 +319,7 @@ void FixQtpieReaxFF::allocate_storage() memory->create(b_prc,nmax,"qtpie:b_prc"); memory->create(b_prm,nmax,"qtpie:b_prm"); - // dual CG support int size = nmax; - if (dual_enabled) size*= 2; memory->create(p,size,"qtpie:p"); memory->create(q,size,"qtpie:q"); @@ -1006,9 +992,6 @@ double FixQtpieReaxFF::memory_usage() bytes += (double)m_cap * sizeof(int); bytes += (double)m_cap * sizeof(double); - if (dual_enabled) - bytes += (double)atom->nmax*4 * sizeof(double); // double size for q, d, r, and p - return bytes; } diff --git a/src/REAXFF/fix_qtpie_reaxff.h b/src/REAXFF/fix_qtpie_reaxff.h index 07afdf9439..827083c1aa 100644 --- a/src/REAXFF/fix_qtpie_reaxff.h +++ b/src/REAXFF/fix_qtpie_reaxff.h @@ -132,9 +132,6 @@ class FixQtpieReaxFF : public Fix { double find_min(double*, int); double distance(const double*, const double*); - - // dual CG support - int dual_enabled; // 0: Original, separate s & t optimization; 1: dual optimization int matvecs_s, matvecs_t; // Iteration count for each system }; From eb6e5b438abc9784a3b654d76aed364402e784fd Mon Sep 17 00:00:00 2001 From: Navraj Lalli Date: Wed, 21 Aug 2024 10:18:51 +0100 Subject: [PATCH 09/40] Remove virtual keyword --- src/REAXFF/fix_qtpie_reaxff.h | 36 +++++++++++++++++------------------ 1 file changed, 18 insertions(+), 18 deletions(-) diff --git a/src/REAXFF/fix_qtpie_reaxff.h b/src/REAXFF/fix_qtpie_reaxff.h index 827083c1aa..1609202693 100644 --- a/src/REAXFF/fix_qtpie_reaxff.h +++ b/src/REAXFF/fix_qtpie_reaxff.h @@ -33,7 +33,7 @@ class FixQtpieReaxFF : public Fix { void post_constructor() override; void init() override; void init_list(int, class NeighList *) override; - virtual void init_storage(); + void init_storage(); void setup_pre_force(int) override; void pre_force(int) override; @@ -89,27 +89,27 @@ class FixQtpieReaxFF : public Fix { int imax, maxwarn; char *pertype_option; // argument to determine how per-type info is obtained - char *gauss_file; // input file for gaussian exponents - double *gauss_exp; // array of gaussian exponents + char *gauss_file; // input file for gaussian exponents + double *gauss_exp; // array of gaussian exponents - virtual void pertype_parameters(char *); + void pertype_parameters(char *); void init_shielding(); void init_taper(); - virtual void allocate_storage(); - virtual void deallocate_storage(); + void allocate_storage(); + void deallocate_storage(); void reallocate_storage(); - virtual void allocate_matrix(); - virtual void deallocate_matrix(); + void allocate_matrix(); + void deallocate_matrix(); void reallocate_matrix(); - virtual void init_matvec(); + void init_matvec(); void init_H(); - virtual void compute_H(); + void compute_H(); double calculate_H(double, double); - virtual void calculate_Q(); + void calculate_Q(); - virtual int CG(double *, double *); - virtual void sparse_matvec(sparse_matrix *, double *, double *); + int CG(double *, double *); + void sparse_matvec(sparse_matrix *, double *, double *); int pack_forward_comm(int, int *, double *, int, int *) override; void unpack_forward_comm(int, int, double *) override; @@ -121,12 +121,12 @@ class FixQtpieReaxFF : public Fix { int pack_exchange(int, double *) override; int unpack_exchange(int, double *) override; - virtual double parallel_norm(double *, int); - virtual double parallel_dot(double *, double *, int); - virtual double parallel_vector_acc(double *, int); + double parallel_norm(double *, int); + double parallel_dot(double *, double *, int); + double parallel_vector_acc(double *, int); - virtual void vector_sum(double *, double, double *, double, double *, int); - virtual void vector_add(double *, double, double *, int); + void vector_sum(double *, double, double *, double, double *, int); + void vector_add(double *, double, double *, int); void calc_chi_eff(); double find_min(double*, int); From 6dd45ccfdb38dd4ae374a3f59a85e352bc3610e5 Mon Sep 17 00:00:00 2001 From: Navraj Lalli Date: Wed, 21 Aug 2024 10:57:00 +0100 Subject: [PATCH 10/40] Add fix_qtpie_reaxff --- src/.gitignore | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/.gitignore b/src/.gitignore index c26eaaba30..33595ed937 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -982,6 +982,8 @@ /fix_qeq_reaxff.h /fix_qmmm.cpp /fix_qmmm.h +/fix_qtpie_reaxff.cpp +/fix_qtpie_reaxff.h /fix_reaxff.cpp /fix_reaxff.h /fix_reaxff_bonds.cpp From c2e4816717fb132854628fd2a4cee082cb51f154 Mon Sep 17 00:00:00 2001 From: Navraj Lalli Date: Wed, 21 Aug 2024 11:32:18 +0100 Subject: [PATCH 11/40] Update contributing authors --- src/REAXFF/fix_qtpie_reaxff.cpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/REAXFF/fix_qtpie_reaxff.cpp b/src/REAXFF/fix_qtpie_reaxff.cpp index 4066103152..8acf9af50c 100644 --- a/src/REAXFF/fix_qtpie_reaxff.cpp +++ b/src/REAXFF/fix_qtpie_reaxff.cpp @@ -13,7 +13,9 @@ ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- - Contributing author: Navraj S Lalli & Efstratios M Kritikos (Imperial College London) + Contributing authors: + Navraj S Lalli (Imperial College London) + Efstratios M Kritikos (California Institute of Technology) ------------------------------------------------------------------------- */ #include "fix_qtpie_reaxff.h" @@ -1158,7 +1160,7 @@ void FixQtpieReaxFF::calc_chi_eff() // use integral pre-screening for overlap calculations const double emin = find_min(gauss_exp,ntypes+1); - const double dist_cutoff = sqrt(pow(emin,-1.)*log(pow(M_PI/(2.*emin),3.)*pow(10.,2.*KSCREEN))); + const double dist_cutoff = sqrt(pow(emin,-1.0)*log(pow(M_PI/(2.0*emin),3.0)*pow(10.0,2.0*KSCREEN))); const double comm_cutoff = MAX(neighbor->cutneighmax,comm->cutghostuser); if(comm_cutoff < dist_cutoff/ANG_TO_BOHRRAD) { From 62b14aa702208ee2e4eda276c2834e869bf4dc19 Mon Sep 17 00:00:00 2001 From: Navraj Lalli Date: Wed, 21 Aug 2024 12:27:14 +0100 Subject: [PATCH 12/40] Remove unused include --- src/REAXFF/fix_qtpie_reaxff.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/REAXFF/fix_qtpie_reaxff.cpp b/src/REAXFF/fix_qtpie_reaxff.cpp index 8acf9af50c..a1ec2e3c10 100644 --- a/src/REAXFF/fix_qtpie_reaxff.cpp +++ b/src/REAXFF/fix_qtpie_reaxff.cpp @@ -28,7 +28,6 @@ #include "fix_efield.h" #include "force.h" #include "group.h" -#include "memory.h" #include "modify.h" #include "neigh_list.h" #include "neighbor.h" From 79cc70c9dad7c43fa8ade7508bf92ff1cac332ee Mon Sep 17 00:00:00 2001 From: Navraj Lalli Date: Wed, 21 Aug 2024 15:37:26 +0100 Subject: [PATCH 13/40] Rename variable for sum of local and ghost atoms --- src/REAXFF/fix_qtpie_reaxff.cpp | 12 ++++++------ src/REAXFF/fix_qtpie_reaxff.h | 2 +- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/REAXFF/fix_qtpie_reaxff.cpp b/src/REAXFF/fix_qtpie_reaxff.cpp index a1ec2e3c10..ea0ca4b782 100644 --- a/src/REAXFF/fix_qtpie_reaxff.cpp +++ b/src/REAXFF/fix_qtpie_reaxff.cpp @@ -99,7 +99,7 @@ FixQtpieReaxFF::FixQtpieReaxFF(LAMMPS *lmp, int narg, char **arg) : } shld = nullptr; - nn = ng = n_cap = 0; + nn = nt = n_cap = 0; nmax = 0; m_fill = m_cap = 0; pack_flag = 0; @@ -546,13 +546,13 @@ void FixQtpieReaxFF::setup_pre_force(int vflag) { if (reaxff) { nn = reaxff->list->inum; - ng = reaxff->list->inum + reaxff->list->gnum; + nt = reaxff->list->inum + reaxff->list->gnum; ilist = reaxff->list->ilist; numneigh = reaxff->list->numneigh; firstneigh = reaxff->list->firstneigh; } else { nn = list->inum; - ng = list->inum + list->gnum; + nt = list->inum + list->gnum; ilist = list->ilist; numneigh = list->numneigh; firstneigh = list->firstneigh; @@ -613,13 +613,13 @@ void FixQtpieReaxFF::pre_force(int /*vflag*/) if (reaxff) { nn = reaxff->list->inum; - ng = reaxff->list->inum + reaxff->list->gnum; + nt = reaxff->list->inum + reaxff->list->gnum; ilist = reaxff->list->ilist; numneigh = reaxff->list->numneigh; firstneigh = reaxff->list->firstneigh; } else { nn = list->inum; - ng = list->inum + list->gnum; + nt = list->inum + list->gnum; ilist = list->ilist; numneigh = list->numneigh; firstneigh = list->firstneigh; @@ -1183,7 +1183,7 @@ void FixQtpieReaxFF::calc_chi_eff() sum_n = 0.0; sum_d = 0.0; - for (j = 0; j < ng; j++) { + for (j = 0; j < nt; j++) { dist = distance(x[i],x[j])*ANG_TO_BOHRRAD; // distance between atoms as a multiple of Bohr radius if (dist < dist_cutoff) { diff --git a/src/REAXFF/fix_qtpie_reaxff.h b/src/REAXFF/fix_qtpie_reaxff.h index 1609202693..523e39ecb5 100644 --- a/src/REAXFF/fix_qtpie_reaxff.h +++ b/src/REAXFF/fix_qtpie_reaxff.h @@ -48,7 +48,7 @@ class FixQtpieReaxFF : public Fix { protected: int nevery, reaxflag; int matvecs; - int nn, ng, m_fill; + int nn, nt, m_fill; int n_cap, nmax, m_cap; int pack_flag; int nlevels_respa; From f3e5e4b4c1d85409553317bf25d5823f40a1bb23 Mon Sep 17 00:00:00 2001 From: Navraj Lalli Date: Wed, 21 Aug 2024 17:31:06 +0100 Subject: [PATCH 14/40] Rename misleading variable name --- src/REAXFF/fix_qtpie_reaxff.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/REAXFF/fix_qtpie_reaxff.cpp b/src/REAXFF/fix_qtpie_reaxff.cpp index ea0ca4b782..cf9046d88d 100644 --- a/src/REAXFF/fix_qtpie_reaxff.cpp +++ b/src/REAXFF/fix_qtpie_reaxff.cpp @@ -48,7 +48,7 @@ using namespace LAMMPS_NS; using namespace FixConst; -static constexpr double EV_TO_KCAL_PER_MOL = 14.4; +static constexpr double CONV_TO_EV = 14.4; static constexpr double SMALL = 1.0e-14; static constexpr double QSUMSMALL = 0.00001; @@ -768,7 +768,7 @@ double FixQtpieReaxFF::calculate_H(double r, double gamma) denom = r * r * r + gamma; denom = pow(denom,1.0/3.0); - return Taper * EV_TO_KCAL_PER_MOL / denom; + return Taper * CONV_TO_EV / denom; } /* ---------------------------------------------------------------------- */ From dff91accb0f45b5605d196801b6b6de24aaf8309 Mon Sep 17 00:00:00 2001 From: Navraj Lalli Date: Fri, 23 Aug 2024 12:59:54 +0100 Subject: [PATCH 15/40] Correct calculation of cut off distance --- src/REAXFF/fix_qtpie_reaxff.cpp | 64 ++++++++++++++++----------------- src/REAXFF/fix_qtpie_reaxff.h | 9 ++--- 2 files changed, 35 insertions(+), 38 deletions(-) diff --git a/src/REAXFF/fix_qtpie_reaxff.cpp b/src/REAXFF/fix_qtpie_reaxff.cpp index cf9046d88d..de4ca30101 100644 --- a/src/REAXFF/fix_qtpie_reaxff.cpp +++ b/src/REAXFF/fix_qtpie_reaxff.cpp @@ -51,18 +51,19 @@ using namespace FixConst; static constexpr double CONV_TO_EV = 14.4; static constexpr double SMALL = 1.0e-14; static constexpr double QSUMSMALL = 0.00001; +static constexpr double ANGSTROM_TO_BOHRRADIUS = 1.8897261259; static const char cite_fix_qtpie_reaxff[] = - "fix qtpie/reaxff command: doi:https://doi.org/10.1016/j.cplett.2007.02.065\n\n" - "@article{chen2007qtpie,\n" - "title={QTPIE: Charge transfer with polarization current equalization. A fluctuating charge model with correct asymptotics},\n" - "author={Chen, Jiahao and Martinez, Todd J},\n" - "journal={Chemical physics letters},\n" - "volume={438},\n" - "number={4-6},\n" - "pages={315--320},\n" - "year={2007},\n" - "publisher={Elsevier}\n" + "fix qtpie/reaxff command: doi\n\n" + "@article{,\n" + "title={},\n" + "author={},\n" + "journal={},\n" + "volume={},\n" + "number={},\n" + "pages={},\n" + "year={},\n" + "publisher={}\n" "}\n\n"; /* ---------------------------------------------------------------------- */ @@ -233,6 +234,12 @@ 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 emin = find_min(gauss_exp,ntypes+1); + const int olap_cut = 10; // overlap integrals are neglected if less than pow(10,-olap_cut) + dist_cutoff = sqrt(1/emin*log(pow(10.0,2.0*olap_cut))); + // read chi, eta and gamma if (utils::strmatch(arg,"^reaxff")) { @@ -1139,8 +1146,6 @@ void FixQtpieReaxFF::calc_chi_eff() { memset(&chi_eff[0],0,atom->nmax*sizeof(double)); - const int KSCREEN = 10; - const double ANG_TO_BOHRRAD = 1.8897259886; // 1 Ang = 1.8897259886 Bohr radius const auto x = (const double * const *)atom->x; const int ntypes = atom->ntypes; const int *type = atom->type; @@ -1148,6 +1153,14 @@ void FixQtpieReaxFF::calc_chi_eff() double dist,overlap,sum_n,sum_d,ea,eb,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); + } + // efield energy is in real units of kcal/mol, factor needed for conversion to eV const double qe2f = force->qe2f; const double factor = 1.0/qe2f; @@ -1157,17 +1170,6 @@ void FixQtpieReaxFF::calc_chi_eff() efield->update_efield_variables(); } - // use integral pre-screening for overlap calculations - const double emin = find_min(gauss_exp,ntypes+1); - const double dist_cutoff = sqrt(pow(emin,-1.0)*log(pow(M_PI/(2.0*emin),3.0)*pow(10.0,2.0*KSCREEN))); - - const double comm_cutoff = MAX(neighbor->cutneighmax,comm->cutghostuser); - if(comm_cutoff < dist_cutoff/ANG_TO_BOHRRAD) { - error->all(FLERR,"comm cutoff = {} Angstrom is smaller than distance cutoff = {} Angstrom " - "for overlap integral in {}. Increase comm cutoff with comm_modify", - comm_cutoff, dist_cutoff/ANG_TO_BOHRRAD, style); - } - // compute chi_eff for each local atom for (i = 0; i < nn; i++) { ea = gauss_exp[type[i]]; @@ -1184,23 +1186,17 @@ void FixQtpieReaxFF::calc_chi_eff() sum_d = 0.0; for (j = 0; j < nt; j++) { - dist = distance(x[i],x[j])*ANG_TO_BOHRRAD; // distance between atoms as a multiple of Bohr radius + dist = distance(x[i],x[j])*ANGSTROM_TO_BOHRRADIUS; // in atomic units if (dist < dist_cutoff) { eb = gauss_exp[type[j]]; chib = chi[type[j]]; - // The expressions below are in atomic units - // Implementation from Chen Jiahao, Theory and applications of fluctuating-charge models, 2009 (with normalization constants added) + // overlap integral of two normalised 1s Gaussian type orbitals p = ea + eb; m = ea * eb / p; overlap = pow((4.0*m/p),0.75) * exp(-m*dist*dist); - // Implementation from T. Halgaker et al., Molecular electronic-structure theory, 2000 -// p = ea + eb; -// m = ea * eb / p; -// Overlap = pow((M_PI / p), 1.5) * exp(-m * R * R); - if (efield) { if (efield->varflag != FixEfield::ATOM) { phib = factor*(x[j][0]*efield->ex + x[j][1]*efield->ey + x[j][2]*efield->ez); @@ -1226,11 +1222,11 @@ void FixQtpieReaxFF::calc_chi_eff() /* ---------------------------------------------------------------------- */ -double FixQtpieReaxFF::find_min(double *array, int array_length) +double FixQtpieReaxFF::find_min(const double *array, const int array_length) { - // since types start from 1, gaussian exponents start from 1 + // index of first gaussian orbital exponent is 1 double smallest = array[1]; - for (int i = 1; i < array_length; i++) + for (int i = 2; i < array_length; i++) { if (array[i] < smallest) smallest = array[i]; diff --git a/src/REAXFF/fix_qtpie_reaxff.h b/src/REAXFF/fix_qtpie_reaxff.h index 523e39ecb5..8b98025a81 100644 --- a/src/REAXFF/fix_qtpie_reaxff.h +++ b/src/REAXFF/fix_qtpie_reaxff.h @@ -82,15 +82,16 @@ class FixQtpieReaxFF : public Fix { double *Hdia_inv; double *b_s, *b_t; double *b_prc, *b_prm; - double *chi_eff; // array of effective electronegativities + double *chi_eff; // array of effective electronegativities //CG storage double *p, *q, *r, *d; int imax, maxwarn; char *pertype_option; // argument to determine how per-type info is obtained - char *gauss_file; // input file for gaussian exponents - double *gauss_exp; // array of gaussian exponents + 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 void pertype_parameters(char *); void init_shielding(); @@ -129,7 +130,7 @@ class FixQtpieReaxFF : public Fix { void vector_add(double *, double, double *, int); void calc_chi_eff(); - double find_min(double*, int); + double find_min(const double*, const int); double distance(const double*, const double*); int matvecs_s, matvecs_t; // Iteration count for each system From ba2217a4b4413a7c4742ece76010e3d91375b0eb Mon Sep 17 00:00:00 2001 From: Navraj Lalli Date: Fri, 23 Aug 2024 15:20:35 +0100 Subject: [PATCH 16/40] Improve exceptions in reading of gauss file --- src/REAXFF/fix_qtpie_reaxff.cpp | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/src/REAXFF/fix_qtpie_reaxff.cpp b/src/REAXFF/fix_qtpie_reaxff.cpp index de4ca30101..ac9415c306 100644 --- a/src/REAXFF/fix_qtpie_reaxff.cpp +++ b/src/REAXFF/fix_qtpie_reaxff.cpp @@ -209,23 +209,27 @@ void FixQtpieReaxFF::pertype_parameters(char *arg) gauss_exp[0] = 0.0; try { TextFileReader reader(gauss_file,"qtpie/reaxff gaussian exponents"); - reader.ignore_comments = false; + reader.ignore_comments = true; for (int i = 1; i <= ntypes; i++) { const char *line = reader.next_line(); - std::cout << "Orbital exponent " << line; if (!line) - throw TokenizerException("Fix qtpie/reaxff: Invalid param file format",""); + throw TokenizerException("Fix qtpie/reaxff: Incorrect number of atom types in gauss file",""); ValueTokenizer values(line); if (values.count() != 2) - throw TokenizerException("Fix qtpie/reaxff: Incorrect format of param file",""); + throw TokenizerException("Fix qtpie/reaxff: Incorrect number of values per line " + "in gauss file",std::to_string(values.count())); int itype = values.next_int(); if ((itype < 1) || (itype > ntypes)) - throw TokenizerException("Fix qtpie/reaxff: Invalid atom type in param file", + throw TokenizerException("Fix qtpie/reaxff: Invalid atom type in gauss file", std::to_string(itype)); - gauss_exp[itype] = values.next_double(); + double expo = values.next_double(); + if (expo < 0) + throw TokenizerException("Fix qtpie/reaxff: Invalid orbital exponent in gauss file", + std::to_string(expo)); + gauss_exp[itype] = expo; } } catch (std::exception &e) { error->one(FLERR,e.what()); From 5fab9e665f3243977f216f072b4c285bbcc248c1 Mon Sep 17 00:00:00 2001 From: Navraj Lalli Date: Mon, 2 Sep 2024 16:59:07 +0100 Subject: [PATCH 17/40] Update with changes made to fix_qeq_reaxff.cpp --- src/REAXFF/fix_qtpie_reaxff.cpp | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/src/REAXFF/fix_qtpie_reaxff.cpp b/src/REAXFF/fix_qtpie_reaxff.cpp index ac9415c306..904fb19029 100644 --- a/src/REAXFF/fix_qtpie_reaxff.cpp +++ b/src/REAXFF/fix_qtpie_reaxff.cpp @@ -372,7 +372,8 @@ void FixQtpieReaxFF::reallocate_storage() void FixQtpieReaxFF::allocate_matrix() { - int i,ii,m; + int i,ii; + bigint m; int mincap; double safezone; @@ -394,14 +395,17 @@ void FixQtpieReaxFF::allocate_matrix() i = ilist[ii]; m += numneigh[i]; } - m_cap = MAX((int)(m * safezone), mincap * REAX_MIN_NBRS); + bigint m_cap_big = (bigint)MAX(m * safezone, mincap * REAX_MIN_NBRS); + if (m_cap_big > MAXSMALLINT) + error->one(FLERR,"Too many neighbors in fix qeq/reaxff"); + m_cap = m_cap_big; H.n = n_cap; H.m = m_cap; - memory->create(H.firstnbr,n_cap,"qtpie:H.firstnbr"); - memory->create(H.numnbrs,n_cap,"qtpie:H.numnbrs"); - memory->create(H.jlist,m_cap,"qtpie:H.jlist"); - memory->create(H.val,m_cap,"qtpie:H.val"); + memory->create(H.firstnbr,n_cap,"qeq:H.firstnbr"); + memory->create(H.numnbrs,n_cap,"qeq:H.numnbrs"); + memory->create(H.jlist,m_cap,"qeq:H.jlist"); + memory->create(H.val,m_cap,"qeq:H.val"); } /* ---------------------------------------------------------------------- */ From 49dcb679f642b6d40f4b1a3353ca5cf9b405611d Mon Sep 17 00:00:00 2001 From: Navraj Lalli Date: Tue, 3 Sep 2024 12:22:10 +0100 Subject: [PATCH 18/40] Change names of orbital exponents --- src/REAXFF/fix_qtpie_reaxff.cpp | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/src/REAXFF/fix_qtpie_reaxff.cpp b/src/REAXFF/fix_qtpie_reaxff.cpp index 904fb19029..ec573f6a57 100644 --- a/src/REAXFF/fix_qtpie_reaxff.cpp +++ b/src/REAXFF/fix_qtpie_reaxff.cpp @@ -203,7 +203,7 @@ void FixQtpieReaxFF::pertype_parameters(char *arg) const int *type = atom->type; const int ntypes = atom->ntypes; - // read gaussian exponents + // read gaussian orbital exponents memory->create(gauss_exp,ntypes+1,"qtpie/reaxff:gauss_exp"); if (comm->me == 0) { gauss_exp[0] = 0.0; @@ -225,11 +225,11 @@ void FixQtpieReaxFF::pertype_parameters(char *arg) throw TokenizerException("Fix qtpie/reaxff: Invalid atom type in gauss file", std::to_string(itype)); - double expo = values.next_double(); - if (expo < 0) + double exp = values.next_double(); + if (exp < 0) throw TokenizerException("Fix qtpie/reaxff: Invalid orbital exponent in gauss file", - std::to_string(expo)); - gauss_exp[itype] = expo; + std::to_string(exp)); + gauss_exp[itype] = exp; } } catch (std::exception &e) { error->one(FLERR,e.what()); @@ -240,9 +240,9 @@ void FixQtpieReaxFF::pertype_parameters(char *arg) // define a cutoff distance (in atomic units) beyond which overlap integrals are neglected // in calc_chi_eff() - const double emin = find_min(gauss_exp,ntypes+1); + const double expmin = find_min(gauss_exp,ntypes+1); const int olap_cut = 10; // overlap integrals are neglected if less than pow(10,-olap_cut) - dist_cutoff = sqrt(1/emin*log(pow(10.0,2.0*olap_cut))); + dist_cutoff = sqrt(2*olap_cut/expmin*log(10.0)); // read chi, eta and gamma @@ -1158,7 +1158,7 @@ void FixQtpieReaxFF::calc_chi_eff() const int ntypes = atom->ntypes; const int *type = atom->type; - double dist,overlap,sum_n,sum_d,ea,eb,chia,chib,phia,phib,p,m; + double dist,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 @@ -1180,7 +1180,7 @@ void FixQtpieReaxFF::calc_chi_eff() // compute chi_eff for each local atom for (i = 0; i < nn; i++) { - ea = gauss_exp[type[i]]; + expa = gauss_exp[type[i]]; chia = chi[type[i]]; if (efield) { if (efield->varflag != FixEfield::ATOM) { @@ -1197,12 +1197,12 @@ void FixQtpieReaxFF::calc_chi_eff() dist = distance(x[i],x[j])*ANGSTROM_TO_BOHRRADIUS; // in atomic units if (dist < dist_cutoff) { - eb = gauss_exp[type[j]]; + expb = gauss_exp[type[j]]; chib = chi[type[j]]; // overlap integral of two normalised 1s Gaussian type orbitals - p = ea + eb; - m = ea * eb / p; + p = expa + expb; + m = expa * expb / p; overlap = pow((4.0*m/p),0.75) * exp(-m*dist*dist); if (efield) { From 3a5e764730afc6a90dc6850109d2643586c73558 Mon Sep 17 00:00:00 2001 From: Navraj Lalli Date: Tue, 3 Sep 2024 17:50:14 +0100 Subject: [PATCH 19/40] Fix whitespace --- src/REAXFF/fix_qtpie_reaxff.cpp | 38 ++++++++++++++++----------------- 1 file changed, 19 insertions(+), 19 deletions(-) diff --git a/src/REAXFF/fix_qtpie_reaxff.cpp b/src/REAXFF/fix_qtpie_reaxff.cpp index ec573f6a57..5fd60e51c5 100644 --- a/src/REAXFF/fix_qtpie_reaxff.cpp +++ b/src/REAXFF/fix_qtpie_reaxff.cpp @@ -13,7 +13,7 @@ ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- - Contributing authors: + Contributing authors: Navraj S Lalli (Imperial College London) Efstratios M Kritikos (California Institute of Technology) ------------------------------------------------------------------------- */ @@ -218,15 +218,15 @@ void FixQtpieReaxFF::pertype_parameters(char *arg) if (values.count() != 2) throw TokenizerException("Fix qtpie/reaxff: Incorrect number of values per line " - "in gauss file",std::to_string(values.count())); + "in gauss file",std::to_string(values.count())); int itype = values.next_int(); if ((itype < 1) || (itype > ntypes)) throw TokenizerException("Fix qtpie/reaxff: Invalid atom type in gauss file", std::to_string(itype)); - double exp = values.next_double(); - if (exp < 0) + double exp = values.next_double(); + if (exp < 0) throw TokenizerException("Fix qtpie/reaxff: Invalid orbital exponent in gauss file", std::to_string(exp)); gauss_exp[itype] = exp; @@ -1165,14 +1165,14 @@ void FixQtpieReaxFF::calc_chi_eff() 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); + "for overlap integrals in {}. Increase comm cutoff with comm_modify", + comm_cutoff, dist_cutoff/ANGSTROM_TO_BOHRRADIUS, style); } // efield energy is in real units of kcal/mol, factor needed for conversion to eV const double qe2f = force->qe2f; const double factor = 1.0/qe2f; - + if (efield) { if (efield->varflag != FixEfield::CONSTANT) efield->update_efield_variables(); @@ -1184,9 +1184,9 @@ void FixQtpieReaxFF::calc_chi_eff() chia = chi[type[i]]; if (efield) { if (efield->varflag != FixEfield::ATOM) { - phia = factor*(x[i][0]*efield->ex + x[i][1]*efield->ey + x[i][2]*efield->ez); + phia = factor*(x[i][0]*efield->ex + x[i][1]*efield->ey + x[i][2]*efield->ez); } else { // atom-style potential from FixEfield - phia = efield->efield[i][3]; + phia = efield->efield[i][3]; } } @@ -1200,22 +1200,22 @@ void FixQtpieReaxFF::calc_chi_eff() expb = gauss_exp[type[j]]; chib = chi[type[j]]; - // overlap integral of two normalised 1s Gaussian type orbitals + // 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); if (efield) { - 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 - chib + phib - phia) * overlap; + 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 - chib + phib - phia) * overlap; } else { - sum_n += (chia - chib) * overlap; + sum_n += (chia - chib) * overlap; } - sum_d += overlap; + sum_d += overlap; } } @@ -1236,7 +1236,7 @@ double FixQtpieReaxFF::find_min(const double *array, const int array_length) double smallest = array[1]; for (int i = 2; i < array_length; i++) { - if (array[i] < smallest) + if (array[i] < smallest) smallest = array[i]; } return smallest; From 5e8ecf9cb4ebc63dcb5d8092b9c7cc940030bf80 Mon Sep 17 00:00:00 2001 From: Navraj Lalli Date: Tue, 10 Sep 2024 14:45:16 +0100 Subject: [PATCH 20/40] Rename variables and function for min exponent --- src/REAXFF/fix_qtpie_reaxff.cpp | 14 +++++++------- src/REAXFF/fix_qtpie_reaxff.h | 2 +- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/REAXFF/fix_qtpie_reaxff.cpp b/src/REAXFF/fix_qtpie_reaxff.cpp index 5fd60e51c5..06b5ff1660 100644 --- a/src/REAXFF/fix_qtpie_reaxff.cpp +++ b/src/REAXFF/fix_qtpie_reaxff.cpp @@ -240,9 +240,9 @@ void FixQtpieReaxFF::pertype_parameters(char *arg) // define a cutoff distance (in atomic units) beyond which overlap integrals are neglected // in calc_chi_eff() - const double expmin = find_min(gauss_exp,ntypes+1); + 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/expmin*log(10.0)); + dist_cutoff = sqrt(2*olap_cut/exp_min*log(10.0)); // read chi, eta and gamma @@ -1230,16 +1230,16 @@ void FixQtpieReaxFF::calc_chi_eff() /* ---------------------------------------------------------------------- */ -double FixQtpieReaxFF::find_min(const double *array, const int array_length) +double FixQtpieReaxFF::find_min_exp(const double *array, const int array_length) { // index of first gaussian orbital exponent is 1 - double smallest = array[1]; + double exp_min = array[1]; for (int i = 2; i < array_length; i++) { - if (array[i] < smallest) - smallest = array[i]; + if (array[i] < exp_min) + exp_min = array[i]; } - return smallest; + return exp_min; } /* ---------------------------------------------------------------------- */ diff --git a/src/REAXFF/fix_qtpie_reaxff.h b/src/REAXFF/fix_qtpie_reaxff.h index 8b98025a81..2a89e6f746 100644 --- a/src/REAXFF/fix_qtpie_reaxff.h +++ b/src/REAXFF/fix_qtpie_reaxff.h @@ -130,7 +130,7 @@ class FixQtpieReaxFF : public Fix { void vector_add(double *, double, double *, int); void calc_chi_eff(); - double find_min(const double*, const int); + double find_min_exp(const double*, const int); double distance(const double*, const double*); int matvecs_s, matvecs_t; // Iteration count for each system From 25f33e8721a801206db1bcc85b25ec6c1fbb0544 Mon Sep 17 00:00:00 2001 From: Navraj Lalli Date: Tue, 10 Sep 2024 16:20:51 +0100 Subject: [PATCH 21/40] Add water examples featuring fix qtpie/reaxff --- examples/reaxff/water/gauss_exp.txt | 5 ++++ examples/reaxff/water/in.water.qtpie | 29 +++++++++++++++++++++ examples/reaxff/water/in.water.qtpie.field | 30 ++++++++++++++++++++++ 3 files changed, 64 insertions(+) create mode 100644 examples/reaxff/water/gauss_exp.txt create mode 100644 examples/reaxff/water/in.water.qtpie create mode 100644 examples/reaxff/water/in.water.qtpie.field diff --git a/examples/reaxff/water/gauss_exp.txt b/examples/reaxff/water/gauss_exp.txt new file mode 100644 index 0000000000..4210471e9f --- /dev/null +++ b/examples/reaxff/water/gauss_exp.txt @@ -0,0 +1,5 @@ +# Gaussian orbital exponents (required for fix qtpie/reaxff) taken from Table 2.2 +# of Chen, J. (2009). Theory and applications of fluctuating-charge models. +# The units of the exponents are 1 / (Bohr radius)^2 . +1 0.2240 # O +2 0.5434 # H diff --git a/examples/reaxff/water/in.water.qtpie b/examples/reaxff/water/in.water.qtpie new file mode 100644 index 0000000000..a8f8759444 --- /dev/null +++ b/examples/reaxff/water/in.water.qtpie @@ -0,0 +1,29 @@ +# QTPIE Water + +boundary p p p +units real +atom_style charge + +read_data data.water + +variable x index 1 +variable y index 1 +variable z index 1 + +replicate $x $y $z + +pair_style reaxff NULL safezone 3.0 mincap 150 +pair_coeff * * qeq_ff.water O H +neighbor 0.5 bin +neigh_modify every 1 delay 0 check yes + +velocity all create 300.0 4928459 rot yes dist gaussian + +fix 1 all qtpie/reaxff 1 0.0 10.0 1.0e-6 reaxff gauss_exp.txt +fix 2 all nvt temp 300 300 50.0 + +timestep 0.5 +thermo 10 +thermo_style custom step temp press density vol + +run 20 diff --git a/examples/reaxff/water/in.water.qtpie.field b/examples/reaxff/water/in.water.qtpie.field new file mode 100644 index 0000000000..e5ac77484f --- /dev/null +++ b/examples/reaxff/water/in.water.qtpie.field @@ -0,0 +1,30 @@ +# QTPIE Water + +boundary p p p +units real +atom_style charge + +read_data data.water + +variable x index 1 +variable y index 1 +variable z index 1 + +replicate $x $y $z + +pair_style reaxff NULL safezone 3.0 mincap 150 +pair_coeff * * qeq_ff.water O H +neighbor 0.5 bin +neigh_modify every 1 delay 0 check yes + +velocity all create 300.0 4928459 rot yes dist gaussian + +fix 1 all qtpie/reaxff 1 0.0 10.0 1.0e-6 reaxff gauss_exp.txt +fix 2 all nvt temp 300 300 50.0 +fix 3 all efield 0.0 0.0 0.05 + +timestep 0.5 +thermo 10 +thermo_style custom step temp press density vol + +run 20 From 6f2c4aaf0b80ef3db943933ddc29291ba3c4f1b5 Mon Sep 17 00:00:00 2001 From: Navraj Lalli Date: Tue, 10 Sep 2024 16:40:19 +0100 Subject: [PATCH 22/40] Remove unused code --- src/REAXFF/fix_qtpie_reaxff.cpp | 26 -------------------------- 1 file changed, 26 deletions(-) diff --git a/src/REAXFF/fix_qtpie_reaxff.cpp b/src/REAXFF/fix_qtpie_reaxff.cpp index 06b5ff1660..5174cfc112 100644 --- a/src/REAXFF/fix_qtpie_reaxff.cpp +++ b/src/REAXFF/fix_qtpie_reaxff.cpp @@ -21,7 +21,6 @@ #include "fix_qtpie_reaxff.h" #include "atom.h" -#include "citeme.h" #include "comm.h" #include "domain.h" #include "error.h" @@ -53,19 +52,6 @@ static constexpr double SMALL = 1.0e-14; static constexpr double QSUMSMALL = 0.00001; static constexpr double ANGSTROM_TO_BOHRRADIUS = 1.8897261259; -static const char cite_fix_qtpie_reaxff[] = - "fix qtpie/reaxff command: doi\n\n" - "@article{,\n" - "title={},\n" - "author={},\n" - "journal={},\n" - "volume={},\n" - "number={},\n" - "pages={},\n" - "year={},\n" - "publisher={}\n" - "}\n\n"; - /* ---------------------------------------------------------------------- */ FixQtpieReaxFF::FixQtpieReaxFF(LAMMPS *lmp, int narg, char **arg) : @@ -173,8 +159,6 @@ FixQtpieReaxFF::~FixQtpieReaxFF() void FixQtpieReaxFF::post_constructor() { - if (lmp->citeme) lmp->citeme->add(cite_fix_qtpie_reaxff); - grow_arrays(atom->nmax); for (int i = 0; i < atom->nmax; i++) for (int j = 0; j < nprev; ++j) @@ -471,16 +455,6 @@ void FixQtpieReaxFF::init() if (efield->varflag == FixEfield::ATOM && efield->pstyle != FixEfield::ATOM) error->all(FLERR,"Atom-style external electric field requires atom-style " "potential variable when used with fix {}", style); - // if (((efield->xstyle != FixEfield::CONSTANT) && domain->xperiodic) || - // ((efield->ystyle != FixEfield::CONSTANT) && domain->yperiodic) || - // ((efield->zstyle != FixEfield::CONSTANT) && domain->zperiodic)) - // error->all(FLERR,"Must not have electric field component in direction of periodic " - // "boundary when using charge equilibration with ReaxFF."); - // if (((fabs(efield->ex) > SMALL) && domain->xperiodic) || - // ((fabs(efield->ey) > SMALL) && domain->yperiodic) || - // ((fabs(efield->ez) > SMALL) && domain->zperiodic)) - // error->all(FLERR,"Must not have electric field component in direction of periodic " - // "boundary when using charge equilibration with ReaxFF."); } // we need a half neighbor list w/ Newton off From d67d23738667519acfaf6a28ddce5c8ed2dc19b5 Mon Sep 17 00:00:00 2001 From: Navraj Lalli Date: Thu, 12 Sep 2024 19:41:12 +0100 Subject: [PATCH 23/40] Update author contributions --- src/REAXFF/fix_qtpie_reaxff.cpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/REAXFF/fix_qtpie_reaxff.cpp b/src/REAXFF/fix_qtpie_reaxff.cpp index 5174cfc112..548bce8cfb 100644 --- a/src/REAXFF/fix_qtpie_reaxff.cpp +++ b/src/REAXFF/fix_qtpie_reaxff.cpp @@ -14,8 +14,10 @@ /* ---------------------------------------------------------------------- Contributing authors: - Navraj S Lalli (Imperial College London) - Efstratios M Kritikos (California Institute of Technology) + Efstratios M Kritikos, California Institute of Technology + (Implemented original version in LAMMMPS Aug 2019) + Navraj S Lalli, Imperial College London + (Reimplemented QTPIE as a new fix in LAMMPS Aug 2024 and extended functionality) ------------------------------------------------------------------------- */ #include "fix_qtpie_reaxff.h" From 62f82a7fe12cd7ae99bdf9e7eb40a4359a5a4af1 Mon Sep 17 00:00:00 2001 From: Navraj Lalli Date: Fri, 13 Sep 2024 15:46:27 +0100 Subject: [PATCH 24/40] Remove additional fix name --- src/REAXFF/fix_qtpie_reaxff.h | 1 - 1 file changed, 1 deletion(-) diff --git a/src/REAXFF/fix_qtpie_reaxff.h b/src/REAXFF/fix_qtpie_reaxff.h index 2a89e6f746..2d82ad197c 100644 --- a/src/REAXFF/fix_qtpie_reaxff.h +++ b/src/REAXFF/fix_qtpie_reaxff.h @@ -13,7 +13,6 @@ #ifdef FIX_CLASS // clang-format off -FixStyle(qtpie/reax,FixQtpieReaxFF); FixStyle(qtpie/reaxff,FixQtpieReaxFF); // clang-format on #else From d56f43b4e616e0b06d64abb38fc9f5d2e9856b6b Mon Sep 17 00:00:00 2001 From: Navraj Lalli Date: Fri, 13 Sep 2024 15:50:44 +0100 Subject: [PATCH 25/40] Remove unnecessary tests --- src/REAXFF/fix_qtpie_reaxff.cpp | 5 ----- 1 file changed, 5 deletions(-) diff --git a/src/REAXFF/fix_qtpie_reaxff.cpp b/src/REAXFF/fix_qtpie_reaxff.cpp index 548bce8cfb..9701704972 100644 --- a/src/REAXFF/fix_qtpie_reaxff.cpp +++ b/src/REAXFF/fix_qtpie_reaxff.cpp @@ -1196,11 +1196,6 @@ void FixQtpieReaxFF::calc_chi_eff() } chi_eff[i] = sum_n / sum_d; - - if (fabs(sum_n) < SMALL && fabs(sum_d) < SMALL) - error->all(FLERR,"Unexpected value: fabs(sum_d) is {}", fabs(sum_d)); - if (fabs(sum_d) < 1.0) - error->all(FLERR,"Unexpected value: fabs(sum_d) is {}", fabs(sum_d)); } } From 8ec010f8ca695b2c7d8a650c982e2037ca8f849d Mon Sep 17 00:00:00 2001 From: Navraj Lalli Date: Fri, 13 Sep 2024 15:54:12 +0100 Subject: [PATCH 26/40] Remove unused header file --- src/REAXFF/fix_qtpie_reaxff.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/REAXFF/fix_qtpie_reaxff.cpp b/src/REAXFF/fix_qtpie_reaxff.cpp index 9701704972..0b2ab30294 100644 --- a/src/REAXFF/fix_qtpie_reaxff.cpp +++ b/src/REAXFF/fix_qtpie_reaxff.cpp @@ -44,7 +44,6 @@ #include #include #include -#include using namespace LAMMPS_NS; using namespace FixConst; From bd07f1e8e04eac8d93e890bd0c1fc2471120f793 Mon Sep 17 00:00:00 2001 From: Navraj Lalli Date: Fri, 13 Sep 2024 15:56:16 +0100 Subject: [PATCH 27/40] Change qeq to qtpie --- src/REAXFF/fix_qtpie_reaxff.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/REAXFF/fix_qtpie_reaxff.cpp b/src/REAXFF/fix_qtpie_reaxff.cpp index 0b2ab30294..c26c77c882 100644 --- a/src/REAXFF/fix_qtpie_reaxff.cpp +++ b/src/REAXFF/fix_qtpie_reaxff.cpp @@ -382,15 +382,15 @@ void FixQtpieReaxFF::allocate_matrix() } bigint m_cap_big = (bigint)MAX(m * safezone, mincap * REAX_MIN_NBRS); if (m_cap_big > MAXSMALLINT) - error->one(FLERR,"Too many neighbors in fix qeq/reaxff"); + error->one(FLERR,"Too many neighbors in fix {}",style); m_cap = m_cap_big; H.n = n_cap; H.m = m_cap; - memory->create(H.firstnbr,n_cap,"qeq:H.firstnbr"); - memory->create(H.numnbrs,n_cap,"qeq:H.numnbrs"); - memory->create(H.jlist,m_cap,"qeq:H.jlist"); - memory->create(H.val,m_cap,"qeq:H.val"); + memory->create(H.firstnbr,n_cap,"qtpie:H.firstnbr"); + memory->create(H.numnbrs,n_cap,"qtpie:H.numnbrs"); + memory->create(H.jlist,m_cap,"qtpie:H.jlist"); + memory->create(H.val,m_cap,"qtpie:H.val"); } /* ---------------------------------------------------------------------- */ From af6efcc5145d377aa964575b2be7a045a8c4749f Mon Sep 17 00:00:00 2001 From: Navraj Lalli Date: Fri, 13 Sep 2024 16:43:13 +0100 Subject: [PATCH 28/40] Add fix qtpie/reaxff documentation --- doc/src/Commands_fix.rst | 1 + doc/src/fix.rst | 1 + doc/src/fix_qtpie_reaxff.rst | 178 ++++++++++++++++++++ doc/utils/sphinx-config/false_positives.txt | 4 + 4 files changed, 184 insertions(+) create mode 100644 doc/src/fix_qtpie_reaxff.rst diff --git a/doc/src/Commands_fix.rst b/doc/src/Commands_fix.rst index d9febcc289..6b75d6779c 100644 --- a/doc/src/Commands_fix.rst +++ b/doc/src/Commands_fix.rst @@ -186,6 +186,7 @@ OPT. * :doc:`qeq/slater ` * :doc:`qmmm ` * :doc:`qtb ` + * :doc:`qtpie/reaxff ` * :doc:`rattle ` * :doc:`reaxff/bonds (k) ` * :doc:`reaxff/species (k) ` diff --git a/doc/src/fix.rst b/doc/src/fix.rst index 4919c226fd..b17906d414 100644 --- a/doc/src/fix.rst +++ b/doc/src/fix.rst @@ -365,6 +365,7 @@ accelerated styles exist. * :doc:`qeq/slater ` - charge equilibration via Slater method * :doc:`qmmm ` - functionality to enable a quantum mechanics/molecular mechanics coupling * :doc:`qtb ` - implement quantum thermal bath scheme +* :doc:`qtpie/reaxff ` - apply QTPIE charge equilibration * :doc:`rattle ` - RATTLE constraints on bonds and/or angles * :doc:`reaxff/bonds ` - write out ReaxFF bond information * :doc:`reaxff/species ` - write out ReaxFF molecule information diff --git a/doc/src/fix_qtpie_reaxff.rst b/doc/src/fix_qtpie_reaxff.rst new file mode 100644 index 0000000000..eb5bb4e982 --- /dev/null +++ b/doc/src/fix_qtpie_reaxff.rst @@ -0,0 +1,178 @@ +.. index:: fix qtpie/reaxff + +fix qtpie/reaxff command +======================== + +Syntax +"""""" + +.. code-block:: LAMMPS + + fix ID group-ID qtpie/reaxff Nevery cutlo cuthi tolerance params gfile args + +* ID, group-ID are documented in :doc:`fix ` command +* qtpie/reaxff = style name of this fix command +* Nevery = perform QTPIE every this many steps +* cutlo,cuthi = lo and hi cutoff for Taper radius +* tolerance = precision to which charges will be equilibrated +* params = reaxff or a filename +* gfile = the name of a file containing Gaussian orbital exponents +* one or more keywords or keyword/value pairs may be appended + + .. parsed-literal:: + + keyword = *maxiter* + *maxiter* N = limit the number of iterations to *N* + +Examples +"""""""" + +.. code-block:: LAMMPS + + fix 1 all qtpie/reaxff 1 0.0 10.0 1.0e-6 reaxff exp.qtpie + fix 1 all qtpie/reaxff 1 0.0 10.0 1.0e-6 params.qtpie exp.qtpie maxiter 500 + +Description +""""""""""" + +The QTPIE charge equilibration method is an extension of the QEq charge +equilibration method. With QTPIE, the partial charges on individual atoms +are computed by minimizing the electrostatic energy of the system in the +same way as the QEq method but where the Mulliken electronegativity, +:math:`\chi_i`, of each atom in the QEq charge equilibration scheme +:ref:`(Rappe and Goddard) ` is replaced with an effective +electronegativity given by :ref:`(Chen) ` + +.. math:: + \chi_{\mathrm{eff},i} = \frac{\sum_{j=1}^{N} (\chi_i - \chi_j + \phi_j - \phi_i) S_{ij}} + {\sum_{m=1}^{N}S_{im}}, + +which acts to penalize long-range charge transfer seen with the QEq charge +equilibration scheme. In this equation, :math:`N` is the number of atoms in +the system, :math:`S_{ij}` is the overlap integral between atom :math:`i` +and atom :math:`j`, and :math:`\phi_i` and :math:`\phi_j` are the electric +potentials at the position of atom :math:`i` and :math:`j` due to +an external electric field, respectively. + +This fix is typically used in conjunction with the ReaxFF force +field model as implemented in the :doc:`pair_style reaxff ` +command, but it can be used with any potential in LAMMPS, so long as it +defines and uses charges on each atom. For more technical details about the +charge equilibration performed by `fix qtpie/reaxff`, which is the same as in +:doc:`fix qeq/reaxff ` except for the use of +:math:`\chi_{\mathrm{eff},i}`, please refer to :ref:`(Aktulga) `. +To be explicit, this fix replaces :math:`\chi_k` of eq. 3 in +:ref:`(Aktulga) ` with :math:`\chi_{\mathrm{eff},k}`. + +This fix requires the Mulliken electronegativity, :math:`\chi`, in eV, the +self-Coulomb potential, :math:`\eta`, in eV, and the shielded Coulomb +constant, :math:`\gamma`, in :math:`\AA^{-1}`. If the *params* setting above +is the word "reaxff", then these are extracted from the +:doc:`pair_style reaxff ` command and the ReaxFF force field +file it reads in. If a file name is specified for *params*, then the +parameters are taken from the specified file and the file must contain +one line for each atom type. The latter form must be used when performing +QTPIE with a non-ReaxFF potential. Each line should be formatted as follows, +ensuring that the parameters are given in units of eV, eV, and :math:`\AA^{-1}`, +respectively: + +.. parsed-literal:: + + itype chi eta gamma + +where *itype* is the atom type from 1 to Ntypes. Note that eta is +defined here as twice the eta value in the ReaxFF file. + +The overlap integrals in the equation for :math:`\chi_{\mathrm{eff},i}` +are computed by using normalized 1s Gaussian type orbitals. The Gaussian +orbital exponents, :math:`\alpha`, that are needed to compute the overlap +integrals are taken from the file given by *gfile*. +This file must contain one line for each atom type and provide the Gaussian +orbital exponent for each atom type in units of inverse square Bohr radius. +Each line should be formatted as follows: + +.. parsed-literal:: + + itype alpha + +Empty lines or any text following the pound sign (#) are ignored. An example +*gfile* for a system with two atom types is + +.. parsed-literal:: + + # An example gfile. Exponents are taken from Table 2.2 of Chen, J. (2009). + # Theory and applications of fluctuating-charge models. + # The units of the exponents are 1 / (Bohr radius)^2 . + 1 0.2240 # O + 2 0.5434 # H + +The optional *maxiter* keyword allows changing the max number +of iterations in the linear solver. The default value is 200. + +.. note:: + + In order to solve the self-consistent equations for electronegativity + equalization, LAMMPS imposes the additional constraint that all the + charges in the fix group must add up to zero. The initial charge + assignments should also satisfy this constraint. LAMMPS will print a + warning if that is not the case. + +Restart, fix_modify, output, run start/stop, minimize info +""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" + +No information about this fix is written to :doc:`binary restart files +`. This fix computes a global scalar (the number of +iterations) for access by various :doc:`output commands `. +No parameter of this fix can be used with the *start/stop* keywords of +the :doc:`run ` command. + +This fix is invoked during :doc:`energy minimization `. + +Restrictions +"""""""""""" + +This fix is part of the REAXFF package. It is only enabled if +LAMMPS was built with that package. See the :doc:`Build package +` page for more info. + +This fix does not correctly handle interactions involving multiple +periodic images of the same atom. Hence, it should not be used for +periodic cell dimensions less than 10 Angstroms. + +This fix may be used in combination with :doc:`fix efield ` +and will apply the external electric field during charge equilibration, +but there may be only one fix efield instance used and the electric field +must be applied to all atoms in the system. Consequently, `fix efield` must +be used with *group-ID* all and must not be used with the keyword *region*. +Equal-style variables can be used for electric field vector +components without any further settings. Atom-style variables can be used +for spatially-varying electric field vector components, but the resulting +electric potential must be specified as an atom-style variable using +the *potential* keyword for `fix efield`. + +Related commands +"""""""""""""""" + +:doc:`pair_style reaxff `, :doc:`fix qeq/reaxff ` + +Default +""""""" + +maxiter 200 + +---------- + +.. _Rappe3: + +**(Rappe)** Rappe and Goddard III, Journal of Physical Chemistry, 95, +3358-3363 (1991). + +.. _qtpie-Chen: + +**(Chen)** Chen, Jiahao. Theory and applications of fluctuating-charge models. +University of Illinois at Urbana-Champaign, 2009. + +.. _qeq-Aktulga2: + +**(Aktulga)** Aktulga, Fogarty, Pandit, Grama, Parallel Computing, 38, +245-259 (2012). diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index cfbddbe5f6..65c1031fcf 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -1303,6 +1303,7 @@ gettimeofday geturl gewald Gezelter +gfile Gflop gfortran ghostneigh @@ -1709,6 +1710,7 @@ Jewett jgissing ji Jiang +Jiahao Jiao jik JIK @@ -2363,6 +2365,7 @@ mui Mukherjee Mulders Müller +Mulliken mult multi multibody @@ -3069,6 +3072,7 @@ qqr qqrd Qsb qtb +qtpie quadratically quadrupolar quadrupole From 3f232caf9b7628abd6feffa768a0c88abf53c1c8 Mon Sep 17 00:00:00 2001 From: Navraj Lalli Date: Fri, 13 Sep 2024 17:13:59 +0100 Subject: [PATCH 29/40] Fix whitespace --- doc/src/fix_qtpie_reaxff.rst | 56 ++++++++++++++++----------------- src/REAXFF/fix_qtpie_reaxff.cpp | 2 +- 2 files changed, 29 insertions(+), 29 deletions(-) diff --git a/doc/src/fix_qtpie_reaxff.rst b/doc/src/fix_qtpie_reaxff.rst index eb5bb4e982..5900c3c6e7 100644 --- a/doc/src/fix_qtpie_reaxff.rst +++ b/doc/src/fix_qtpie_reaxff.rst @@ -37,41 +37,41 @@ Description The QTPIE charge equilibration method is an extension of the QEq charge equilibration method. With QTPIE, the partial charges on individual atoms -are computed by minimizing the electrostatic energy of the system in the -same way as the QEq method but where the Mulliken electronegativity, -:math:`\chi_i`, of each atom in the QEq charge equilibration scheme -:ref:`(Rappe and Goddard) ` is replaced with an effective +are computed by minimizing the electrostatic energy of the system in the +same way as the QEq method but where the Mulliken electronegativity, +:math:`\chi_i`, of each atom in the QEq charge equilibration scheme +:ref:`(Rappe and Goddard) ` is replaced with an effective electronegativity given by :ref:`(Chen) ` .. math:: \chi_{\mathrm{eff},i} = \frac{\sum_{j=1}^{N} (\chi_i - \chi_j + \phi_j - \phi_i) S_{ij}} {\sum_{m=1}^{N}S_{im}}, -which acts to penalize long-range charge transfer seen with the QEq charge +which acts to penalize long-range charge transfer seen with the QEq charge equilibration scheme. In this equation, :math:`N` is the number of atoms in -the system, :math:`S_{ij}` is the overlap integral between atom :math:`i` -and atom :math:`j`, and :math:`\phi_i` and :math:`\phi_j` are the electric +the system, :math:`S_{ij}` is the overlap integral between atom :math:`i` +and atom :math:`j`, and :math:`\phi_i` and :math:`\phi_j` are the electric potentials at the position of atom :math:`i` and :math:`j` due to -an external electric field, respectively. +an external electric field, respectively. This fix is typically used in conjunction with the ReaxFF force field model as implemented in the :doc:`pair_style reaxff ` command, but it can be used with any potential in LAMMPS, so long as it -defines and uses charges on each atom. For more technical details about the +defines and uses charges on each atom. For more technical details about the charge equilibration performed by `fix qtpie/reaxff`, which is the same as in -:doc:`fix qeq/reaxff ` except for the use of +:doc:`fix qeq/reaxff ` except for the use of :math:`\chi_{\mathrm{eff},i}`, please refer to :ref:`(Aktulga) `. -To be explicit, this fix replaces :math:`\chi_k` of eq. 3 in +To be explicit, this fix replaces :math:`\chi_k` of eq. 3 in :ref:`(Aktulga) ` with :math:`\chi_{\mathrm{eff},k}`. This fix requires the Mulliken electronegativity, :math:`\chi`, in eV, the self-Coulomb potential, :math:`\eta`, in eV, and the shielded Coulomb -constant, :math:`\gamma`, in :math:`\AA^{-1}`. If the *params* setting above -is the word "reaxff", then these are extracted from the -:doc:`pair_style reaxff ` command and the ReaxFF force field -file it reads in. If a file name is specified for *params*, then the -parameters are taken from the specified file and the file must contain -one line for each atom type. The latter form must be used when performing +constant, :math:`\gamma`, in :math:`\AA^{-1}`. If the *params* setting above +is the word "reaxff", then these are extracted from the +:doc:`pair_style reaxff ` command and the ReaxFF force field +file it reads in. If a file name is specified for *params*, then the +parameters are taken from the specified file and the file must contain +one line for each atom type. The latter form must be used when performing QTPIE with a non-ReaxFF potential. Each line should be formatted as follows, ensuring that the parameters are given in units of eV, eV, and :math:`\AA^{-1}`, respectively: @@ -80,15 +80,15 @@ respectively: itype chi eta gamma -where *itype* is the atom type from 1 to Ntypes. Note that eta is +where *itype* is the atom type from 1 to Ntypes. Note that eta is defined here as twice the eta value in the ReaxFF file. The overlap integrals in the equation for :math:`\chi_{\mathrm{eff},i}` are computed by using normalized 1s Gaussian type orbitals. The Gaussian -orbital exponents, :math:`\alpha`, that are needed to compute the overlap -integrals are taken from the file given by *gfile*. +orbital exponents, :math:`\alpha`, that are needed to compute the overlap +integrals are taken from the file given by *gfile*. This file must contain one line for each atom type and provide the Gaussian -orbital exponent for each atom type in units of inverse square Bohr radius. +orbital exponent for each atom type in units of inverse square Bohr radius. Each line should be formatted as follows: .. parsed-literal:: @@ -100,7 +100,7 @@ Empty lines or any text following the pound sign (#) are ignored. An example .. parsed-literal:: - # An example gfile. Exponents are taken from Table 2.2 of Chen, J. (2009). + # An example gfile. Exponents are taken from Table 2.2 of Chen, J. (2009). # Theory and applications of fluctuating-charge models. # The units of the exponents are 1 / (Bohr radius)^2 . 1 0.2240 # O @@ -142,12 +142,12 @@ periodic cell dimensions less than 10 Angstroms. This fix may be used in combination with :doc:`fix efield ` and will apply the external electric field during charge equilibration, but there may be only one fix efield instance used and the electric field -must be applied to all atoms in the system. Consequently, `fix efield` must +must be applied to all atoms in the system. Consequently, `fix efield` must be used with *group-ID* all and must not be used with the keyword *region*. -Equal-style variables can be used for electric field vector -components without any further settings. Atom-style variables can be used -for spatially-varying electric field vector components, but the resulting -electric potential must be specified as an atom-style variable using +Equal-style variables can be used for electric field vector +components without any further settings. Atom-style variables can be used +for spatially-varying electric field vector components, but the resulting +electric potential must be specified as an atom-style variable using the *potential* keyword for `fix efield`. Related commands @@ -169,7 +169,7 @@ maxiter 200 .. _qtpie-Chen: -**(Chen)** Chen, Jiahao. Theory and applications of fluctuating-charge models. +**(Chen)** Chen, Jiahao. Theory and applications of fluctuating-charge models. University of Illinois at Urbana-Champaign, 2009. .. _qeq-Aktulga2: diff --git a/src/REAXFF/fix_qtpie_reaxff.cpp b/src/REAXFF/fix_qtpie_reaxff.cpp index c26c77c882..b08c6808ac 100644 --- a/src/REAXFF/fix_qtpie_reaxff.cpp +++ b/src/REAXFF/fix_qtpie_reaxff.cpp @@ -14,7 +14,7 @@ /* ---------------------------------------------------------------------- Contributing authors: - Efstratios M Kritikos, California Institute of Technology + Efstratios M Kritikos, California Institute of Technology (Implemented original version in LAMMMPS Aug 2019) Navraj S Lalli, Imperial College London (Reimplemented QTPIE as a new fix in LAMMPS Aug 2024 and extended functionality) From 96c776c51f6547a66431f7986213dc457959b9fb Mon Sep 17 00:00:00 2001 From: Navraj Lalli Date: Fri, 13 Sep 2024 18:08:14 +0100 Subject: [PATCH 30/40] Add log files for qtpie/reaxff examples --- ...log.29Aug24.reaxff.water-qtpie-field.g++.1 | 127 ++++++++++++++++++ ...log.29Aug24.reaxff.water-qtpie-field.g++.4 | 127 ++++++++++++++++++ .../log.29Aug24.reaxff.water-qtpie.g++.1 | 126 +++++++++++++++++ .../log.29Aug24.reaxff.water-qtpie.g++.4 | 126 +++++++++++++++++ 4 files changed, 506 insertions(+) create mode 100644 examples/reaxff/water/log.29Aug24.reaxff.water-qtpie-field.g++.1 create mode 100644 examples/reaxff/water/log.29Aug24.reaxff.water-qtpie-field.g++.4 create mode 100644 examples/reaxff/water/log.29Aug24.reaxff.water-qtpie.g++.1 create mode 100644 examples/reaxff/water/log.29Aug24.reaxff.water-qtpie.g++.4 diff --git a/examples/reaxff/water/log.29Aug24.reaxff.water-qtpie-field.g++.1 b/examples/reaxff/water/log.29Aug24.reaxff.water-qtpie-field.g++.1 new file mode 100644 index 0000000000..33221ff080 --- /dev/null +++ b/examples/reaxff/water/log.29Aug24.reaxff.water-qtpie-field.g++.1 @@ -0,0 +1,127 @@ +LAMMPS (29 Aug 2024 - Development - patch_29Aug2024-76-g3f232caf9b) + using 1 OpenMP thread(s) per MPI task +# QTPIE Water + +boundary p p p +units real +atom_style charge + +read_data data.water +Reading data file ... + orthogonal box = (0 0 0) to (31.043046 31.043046 31.043046) + 1 by 1 by 1 MPI processor grid + reading atoms ... + 3000 atoms + read_data CPU = 0.056 seconds + +variable x index 1 +variable y index 1 +variable z index 1 + +replicate $x $y $z +replicate 1 $y $z +replicate 1 1 $z +replicate 1 1 1 +Replication is creating a 1x1x1 = 1 times larger system... + orthogonal box = (0 0 0) to (31.043046 31.043046 31.043046) + 1 by 1 by 1 MPI processor grid + 3000 atoms + replicate CPU = 0.001 seconds + +pair_style reaxff NULL safezone 3.0 mincap 150 +pair_coeff * * qeq_ff.water O H +WARNING: Changed valency_val to valency_boc for X (src/REAXFF/reaxff_ffield.cpp:294) +neighbor 0.5 bin +neigh_modify every 1 delay 0 check yes + +velocity all create 300.0 4928459 rot yes dist gaussian + +fix 1 all qtpie/reaxff 1 0.0 10.0 1.0e-6 reaxff gauss_exp.txt +fix 2 all nvt temp 300 300 50.0 +fix 3 all efield 0.0 0.0 0.05 + +timestep 0.5 +thermo 10 +thermo_style custom step temp press density vol + +run 20 + +CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE + +Your simulation uses code contributions which should be cited: + +- Type Label Framework: https://doi.org/10.1021/acs.jpcb.3c08419 + +@Article{Gissinger24, + author = {Jacob R. Gissinger, Ilia Nikiforov, Yaser Afshar, Brendon Waters, Moon-ki Choi, Daniel S. Karls, Alexander Stukowski, Wonpil Im, Hendrik Heinz, Axel Kohlmeyer, and Ellad B. Tadmor}, + title = {Type Label Framework for Bonded Force Fields in LAMMPS}, + journal = {J. Phys. Chem. B}, + year = 2024, + volume = 128, + number = 13, + pages = {3282–-3297} +} + +- pair reaxff command: doi:10.1016/j.parco.2011.08.005 + +@Article{Aktulga12, + author = {H. M. Aktulga and J. C. Fogarty and S. A. Pandit and A. Y. Grama}, + title = {Parallel Reactive Molecular Dynamics: {N}umerical Methods and Algorithmic Techniques}, + journal = {Parallel Computing}, + year = 2012, + volume = 38, + number = {4--5}, + pages = {245--259} +} + +CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE + +Neighbor list info ... + update: every = 1 steps, delay = 0 steps, check = yes + max neighbors/atom: 2000, page size: 100000 + master list distance cutoff = 10.5 + ghost atom cutoff = 10.5 + binsize = 5.25, bins = 6 6 6 + 2 neighbor lists, perpetual/occasional/extra = 2 0 0 + (1) pair reaxff, perpetual + attributes: half, newton off, ghost + pair build: half/bin/ghost/newtoff + stencil: full/ghost/bin/3d + bin: standard + (2) fix qtpie/reaxff, perpetual, copy from (1) + attributes: half, newton off + pair build: copy + stencil: none + bin: none +Per MPI rank memory allocation (min/avg/max) = 539.2 | 539.2 | 539.2 Mbytes + Step Temp Press Density Volume + 0 300 10137.041 1 29915.273 + 10 296.09128 3564.7969 1 29915.273 + 20 293.04308 10299.201 1 29915.273 +Loop time of 10.7863 on 1 procs for 20 steps with 3000 atoms + +Performance: 0.080 ns/day, 299.620 hours/ns, 1.854 timesteps/s, 5.563 katom-step/s +100.0% CPU use with 1 MPI tasks x 1 OpenMP threads + +MPI task timing breakdown: +Section | min time | avg time | max time |%varavg| %total +--------------------------------------------------------------- +Pair | 4.7275 | 4.7275 | 4.7275 | 0.0 | 43.83 +Neigh | 0.17533 | 0.17533 | 0.17533 | 0.0 | 1.63 +Comm | 0.0017376 | 0.0017376 | 0.0017376 | 0.0 | 0.02 +Output | 8.2065e-05 | 8.2065e-05 | 8.2065e-05 | 0.0 | 0.00 +Modify | 5.8812 | 5.8812 | 5.8812 | 0.0 | 54.52 +Other | | 0.0005226 | | | 0.00 + +Nlocal: 3000 ave 3000 max 3000 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Nghost: 11077 ave 11077 max 11077 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Neighs: 971775 ave 971775 max 971775 min +Histogram: 1 0 0 0 0 0 0 0 0 0 + +Total # of neighbors = 971775 +Ave neighs/atom = 323.925 +Neighbor list builds = 2 +Dangerous builds = 0 +Total wall time: 0:00:12 diff --git a/examples/reaxff/water/log.29Aug24.reaxff.water-qtpie-field.g++.4 b/examples/reaxff/water/log.29Aug24.reaxff.water-qtpie-field.g++.4 new file mode 100644 index 0000000000..07a348604e --- /dev/null +++ b/examples/reaxff/water/log.29Aug24.reaxff.water-qtpie-field.g++.4 @@ -0,0 +1,127 @@ +LAMMPS (29 Aug 2024 - Development - patch_29Aug2024-76-g3f232caf9b) + using 1 OpenMP thread(s) per MPI task +# QTPIE Water + +boundary p p p +units real +atom_style charge + +read_data data.water +Reading data file ... + orthogonal box = (0 0 0) to (31.043046 31.043046 31.043046) + 1 by 2 by 2 MPI processor grid + reading atoms ... + 3000 atoms + read_data CPU = 0.053 seconds + +variable x index 1 +variable y index 1 +variable z index 1 + +replicate $x $y $z +replicate 1 $y $z +replicate 1 1 $z +replicate 1 1 1 +Replication is creating a 1x1x1 = 1 times larger system... + orthogonal box = (0 0 0) to (31.043046 31.043046 31.043046) + 1 by 2 by 2 MPI processor grid + 3000 atoms + replicate CPU = 0.002 seconds + +pair_style reaxff NULL safezone 3.0 mincap 150 +pair_coeff * * qeq_ff.water O H +WARNING: Changed valency_val to valency_boc for X (src/REAXFF/reaxff_ffield.cpp:294) +neighbor 0.5 bin +neigh_modify every 1 delay 0 check yes + +velocity all create 300.0 4928459 rot yes dist gaussian + +fix 1 all qtpie/reaxff 1 0.0 10.0 1.0e-6 reaxff gauss_exp.txt +fix 2 all nvt temp 300 300 50.0 +fix 3 all efield 0.0 0.0 0.05 + +timestep 0.5 +thermo 10 +thermo_style custom step temp press density vol + +run 20 + +CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE + +Your simulation uses code contributions which should be cited: + +- Type Label Framework: https://doi.org/10.1021/acs.jpcb.3c08419 + +@Article{Gissinger24, + author = {Jacob R. Gissinger, Ilia Nikiforov, Yaser Afshar, Brendon Waters, Moon-ki Choi, Daniel S. Karls, Alexander Stukowski, Wonpil Im, Hendrik Heinz, Axel Kohlmeyer, and Ellad B. Tadmor}, + title = {Type Label Framework for Bonded Force Fields in LAMMPS}, + journal = {J. Phys. Chem. B}, + year = 2024, + volume = 128, + number = 13, + pages = {3282–-3297} +} + +- pair reaxff command: doi:10.1016/j.parco.2011.08.005 + +@Article{Aktulga12, + author = {H. M. Aktulga and J. C. Fogarty and S. A. Pandit and A. Y. Grama}, + title = {Parallel Reactive Molecular Dynamics: {N}umerical Methods and Algorithmic Techniques}, + journal = {Parallel Computing}, + year = 2012, + volume = 38, + number = {4--5}, + pages = {245--259} +} + +CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE + +Neighbor list info ... + update: every = 1 steps, delay = 0 steps, check = yes + max neighbors/atom: 2000, page size: 100000 + master list distance cutoff = 10.5 + ghost atom cutoff = 10.5 + binsize = 5.25, bins = 6 6 6 + 2 neighbor lists, perpetual/occasional/extra = 2 0 0 + (1) pair reaxff, perpetual + attributes: half, newton off, ghost + pair build: half/bin/ghost/newtoff + stencil: full/ghost/bin/3d + bin: standard + (2) fix qtpie/reaxff, perpetual, copy from (1) + attributes: half, newton off + pair build: copy + stencil: none + bin: none +Per MPI rank memory allocation (min/avg/max) = 260.5 | 262.2 | 263.6 Mbytes + Step Temp Press Density Volume + 0 300 10137.041 1 29915.273 + 10 296.09128 3564.7969 1 29915.273 + 20 293.04308 10299.201 1 29915.273 +Loop time of 3.14492 on 4 procs for 20 steps with 3000 atoms + +Performance: 0.275 ns/day, 87.359 hours/ns, 6.359 timesteps/s, 19.078 katom-step/s +99.6% CPU use with 4 MPI tasks x 1 OpenMP threads + +MPI task timing breakdown: +Section | min time | avg time | max time |%varavg| %total +--------------------------------------------------------------- +Pair | 1.6557 | 1.6847 | 1.7281 | 2.1 | 53.57 +Neigh | 0.086503 | 0.086968 | 0.087627 | 0.2 | 2.77 +Comm | 0.003309 | 0.046699 | 0.075729 | 12.4 | 1.48 +Output | 5.0156e-05 | 5.483e-05 | 6.8111e-05 | 0.0 | 0.00 +Modify | 1.3254 | 1.3261 | 1.3266 | 0.0 | 42.16 +Other | | 0.0004552 | | | 0.01 + +Nlocal: 750 ave 760 max 735 min +Histogram: 1 0 0 0 1 0 0 0 0 2 +Nghost: 6230.5 ave 6253 max 6193 min +Histogram: 1 0 0 0 0 0 1 0 1 1 +Neighs: 276995 ave 280886 max 271360 min +Histogram: 1 0 0 0 1 0 0 0 1 1 + +Total # of neighbors = 1107981 +Ave neighs/atom = 369.327 +Neighbor list builds = 2 +Dangerous builds = 0 +Total wall time: 0:00:03 diff --git a/examples/reaxff/water/log.29Aug24.reaxff.water-qtpie.g++.1 b/examples/reaxff/water/log.29Aug24.reaxff.water-qtpie.g++.1 new file mode 100644 index 0000000000..1187a755ee --- /dev/null +++ b/examples/reaxff/water/log.29Aug24.reaxff.water-qtpie.g++.1 @@ -0,0 +1,126 @@ +LAMMPS (29 Aug 2024 - Development - patch_29Aug2024-76-g3f232caf9b) + using 1 OpenMP thread(s) per MPI task +# QTPIE Water + +boundary p p p +units real +atom_style charge + +read_data data.water +Reading data file ... + orthogonal box = (0 0 0) to (31.043046 31.043046 31.043046) + 1 by 1 by 1 MPI processor grid + reading atoms ... + 3000 atoms + read_data CPU = 0.055 seconds + +variable x index 1 +variable y index 1 +variable z index 1 + +replicate $x $y $z +replicate 1 $y $z +replicate 1 1 $z +replicate 1 1 1 +Replication is creating a 1x1x1 = 1 times larger system... + orthogonal box = (0 0 0) to (31.043046 31.043046 31.043046) + 1 by 1 by 1 MPI processor grid + 3000 atoms + replicate CPU = 0.001 seconds + +pair_style reaxff NULL safezone 3.0 mincap 150 +pair_coeff * * qeq_ff.water O H +WARNING: Changed valency_val to valency_boc for X (src/REAXFF/reaxff_ffield.cpp:294) +neighbor 0.5 bin +neigh_modify every 1 delay 0 check yes + +velocity all create 300.0 4928459 rot yes dist gaussian + +fix 1 all qtpie/reaxff 1 0.0 10.0 1.0e-6 reaxff gauss_exp.txt +fix 2 all nvt temp 300 300 50.0 + +timestep 0.5 +thermo 10 +thermo_style custom step temp press density vol + +run 20 + +CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE + +Your simulation uses code contributions which should be cited: + +- Type Label Framework: https://doi.org/10.1021/acs.jpcb.3c08419 + +@Article{Gissinger24, + author = {Jacob R. Gissinger, Ilia Nikiforov, Yaser Afshar, Brendon Waters, Moon-ki Choi, Daniel S. Karls, Alexander Stukowski, Wonpil Im, Hendrik Heinz, Axel Kohlmeyer, and Ellad B. Tadmor}, + title = {Type Label Framework for Bonded Force Fields in LAMMPS}, + journal = {J. Phys. Chem. B}, + year = 2024, + volume = 128, + number = 13, + pages = {3282–-3297} +} + +- pair reaxff command: doi:10.1016/j.parco.2011.08.005 + +@Article{Aktulga12, + author = {H. M. Aktulga and J. C. Fogarty and S. A. Pandit and A. Y. Grama}, + title = {Parallel Reactive Molecular Dynamics: {N}umerical Methods and Algorithmic Techniques}, + journal = {Parallel Computing}, + year = 2012, + volume = 38, + number = {4--5}, + pages = {245--259} +} + +CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE + +Neighbor list info ... + update: every = 1 steps, delay = 0 steps, check = yes + max neighbors/atom: 2000, page size: 100000 + master list distance cutoff = 10.5 + ghost atom cutoff = 10.5 + binsize = 5.25, bins = 6 6 6 + 2 neighbor lists, perpetual/occasional/extra = 2 0 0 + (1) pair reaxff, perpetual + attributes: half, newton off, ghost + pair build: half/bin/ghost/newtoff + stencil: full/ghost/bin/3d + bin: standard + (2) fix qtpie/reaxff, perpetual, copy from (1) + attributes: half, newton off + pair build: copy + stencil: none + bin: none +Per MPI rank memory allocation (min/avg/max) = 539.2 | 539.2 | 539.2 Mbytes + Step Temp Press Density Volume + 0 300 10138.375 1 29915.273 + 10 295.97879 3575.2769 1 29915.273 + 20 292.76583 10309.128 1 29915.273 +Loop time of 10.8138 on 1 procs for 20 steps with 3000 atoms + +Performance: 0.080 ns/day, 300.383 hours/ns, 1.849 timesteps/s, 5.548 katom-step/s +99.9% CPU use with 1 MPI tasks x 1 OpenMP threads + +MPI task timing breakdown: +Section | min time | avg time | max time |%varavg| %total +--------------------------------------------------------------- +Pair | 4.7177 | 4.7177 | 4.7177 | 0.0 | 43.63 +Neigh | 0.17607 | 0.17607 | 0.17607 | 0.0 | 1.63 +Comm | 0.0017295 | 0.0017295 | 0.0017295 | 0.0 | 0.02 +Output | 8.5431e-05 | 8.5431e-05 | 8.5431e-05 | 0.0 | 0.00 +Modify | 5.9177 | 5.9177 | 5.9177 | 0.0 | 54.72 +Other | | 0.0004911 | | | 0.00 + +Nlocal: 3000 ave 3000 max 3000 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Nghost: 11077 ave 11077 max 11077 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Neighs: 971830 ave 971830 max 971830 min +Histogram: 1 0 0 0 0 0 0 0 0 0 + +Total # of neighbors = 971830 +Ave neighs/atom = 323.94333 +Neighbor list builds = 2 +Dangerous builds = 0 +Total wall time: 0:00:12 diff --git a/examples/reaxff/water/log.29Aug24.reaxff.water-qtpie.g++.4 b/examples/reaxff/water/log.29Aug24.reaxff.water-qtpie.g++.4 new file mode 100644 index 0000000000..372156b6a2 --- /dev/null +++ b/examples/reaxff/water/log.29Aug24.reaxff.water-qtpie.g++.4 @@ -0,0 +1,126 @@ +LAMMPS (29 Aug 2024 - Development - patch_29Aug2024-76-g3f232caf9b) + using 1 OpenMP thread(s) per MPI task +# QTPIE Water + +boundary p p p +units real +atom_style charge + +read_data data.water +Reading data file ... + orthogonal box = (0 0 0) to (31.043046 31.043046 31.043046) + 1 by 2 by 2 MPI processor grid + reading atoms ... + 3000 atoms + read_data CPU = 0.053 seconds + +variable x index 1 +variable y index 1 +variable z index 1 + +replicate $x $y $z +replicate 1 $y $z +replicate 1 1 $z +replicate 1 1 1 +Replication is creating a 1x1x1 = 1 times larger system... + orthogonal box = (0 0 0) to (31.043046 31.043046 31.043046) + 1 by 2 by 2 MPI processor grid + 3000 atoms + replicate CPU = 0.002 seconds + +pair_style reaxff NULL safezone 3.0 mincap 150 +pair_coeff * * qeq_ff.water O H +WARNING: Changed valency_val to valency_boc for X (src/REAXFF/reaxff_ffield.cpp:294) +neighbor 0.5 bin +neigh_modify every 1 delay 0 check yes + +velocity all create 300.0 4928459 rot yes dist gaussian + +fix 1 all qtpie/reaxff 1 0.0 10.0 1.0e-6 reaxff gauss_exp.txt +fix 2 all nvt temp 300 300 50.0 + +timestep 0.5 +thermo 10 +thermo_style custom step temp press density vol + +run 20 + +CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE + +Your simulation uses code contributions which should be cited: + +- Type Label Framework: https://doi.org/10.1021/acs.jpcb.3c08419 + +@Article{Gissinger24, + author = {Jacob R. Gissinger, Ilia Nikiforov, Yaser Afshar, Brendon Waters, Moon-ki Choi, Daniel S. Karls, Alexander Stukowski, Wonpil Im, Hendrik Heinz, Axel Kohlmeyer, and Ellad B. Tadmor}, + title = {Type Label Framework for Bonded Force Fields in LAMMPS}, + journal = {J. Phys. Chem. B}, + year = 2024, + volume = 128, + number = 13, + pages = {3282–-3297} +} + +- pair reaxff command: doi:10.1016/j.parco.2011.08.005 + +@Article{Aktulga12, + author = {H. M. Aktulga and J. C. Fogarty and S. A. Pandit and A. Y. Grama}, + title = {Parallel Reactive Molecular Dynamics: {N}umerical Methods and Algorithmic Techniques}, + journal = {Parallel Computing}, + year = 2012, + volume = 38, + number = {4--5}, + pages = {245--259} +} + +CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE + +Neighbor list info ... + update: every = 1 steps, delay = 0 steps, check = yes + max neighbors/atom: 2000, page size: 100000 + master list distance cutoff = 10.5 + ghost atom cutoff = 10.5 + binsize = 5.25, bins = 6 6 6 + 2 neighbor lists, perpetual/occasional/extra = 2 0 0 + (1) pair reaxff, perpetual + attributes: half, newton off, ghost + pair build: half/bin/ghost/newtoff + stencil: full/ghost/bin/3d + bin: standard + (2) fix qtpie/reaxff, perpetual, copy from (1) + attributes: half, newton off + pair build: copy + stencil: none + bin: none +Per MPI rank memory allocation (min/avg/max) = 260.5 | 262.2 | 263.6 Mbytes + Step Temp Press Density Volume + 0 300 10138.375 1 29915.273 + 10 295.97879 3575.2769 1 29915.273 + 20 292.76583 10309.128 1 29915.273 +Loop time of 3.13598 on 4 procs for 20 steps with 3000 atoms + +Performance: 0.276 ns/day, 87.111 hours/ns, 6.378 timesteps/s, 19.133 katom-step/s +99.6% CPU use with 4 MPI tasks x 1 OpenMP threads + +MPI task timing breakdown: +Section | min time | avg time | max time |%varavg| %total +--------------------------------------------------------------- +Pair | 1.6622 | 1.695 | 1.7252 | 2.2 | 54.05 +Neigh | 0.086543 | 0.087117 | 0.087848 | 0.2 | 2.78 +Comm | 0.0048192 | 0.035002 | 0.067754 | 15.4 | 1.12 +Output | 4.8033e-05 | 5.3375e-05 | 6.6893e-05 | 0.0 | 0.00 +Modify | 1.3176 | 1.3183 | 1.3189 | 0.0 | 42.04 +Other | | 0.0004753 | | | 0.02 + +Nlocal: 750 ave 760 max 735 min +Histogram: 1 0 0 0 1 0 0 0 0 2 +Nghost: 6229.5 ave 6253 max 6191 min +Histogram: 1 0 0 0 0 0 1 0 1 1 +Neighs: 277011 ave 280900 max 271380 min +Histogram: 1 0 0 0 1 0 0 0 1 1 + +Total # of neighbors = 1108044 +Ave neighs/atom = 369.348 +Neighbor list builds = 2 +Dangerous builds = 0 +Total wall time: 0:00:03 From d64be895e638ffb1ac6ec389e3c66ce87c68eae4 Mon Sep 17 00:00:00 2001 From: Navraj Lalli Date: Thu, 26 Sep 2024 16:42:01 +0100 Subject: [PATCH 31/40] Allow for output of effective electronegativities --- doc/src/fix_qtpie_reaxff.rst | 3 ++- src/REAXFF/fix_qtpie_reaxff.cpp | 7 +++++++ 2 files changed, 9 insertions(+), 1 deletion(-) diff --git a/doc/src/fix_qtpie_reaxff.rst b/doc/src/fix_qtpie_reaxff.rst index 5900c3c6e7..e0c2de1432 100644 --- a/doc/src/fix_qtpie_reaxff.rst +++ b/doc/src/fix_qtpie_reaxff.rst @@ -122,7 +122,8 @@ Restart, fix_modify, output, run start/stop, minimize info No information about this fix is written to :doc:`binary restart files `. This fix computes a global scalar (the number of -iterations) for access by various :doc:`output commands `. +iterations) and a per-atom vector (the effective electronegativity), which +can be accessed by various :doc:`output commands `. No parameter of this fix can be used with the *start/stop* keywords of the :doc:`run ` command. diff --git a/src/REAXFF/fix_qtpie_reaxff.cpp b/src/REAXFF/fix_qtpie_reaxff.cpp index b08c6808ac..30a7b7f71f 100644 --- a/src/REAXFF/fix_qtpie_reaxff.cpp +++ b/src/REAXFF/fix_qtpie_reaxff.cpp @@ -58,8 +58,14 @@ static constexpr double ANGSTROM_TO_BOHRRADIUS = 1.8897261259; FixQtpieReaxFF::FixQtpieReaxFF(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), matvecs(0), pertype_option(nullptr), gauss_file(nullptr) { + // this fix returns a global scalar (the number of iterations) scalar_flag = 1; extscalar = 0; + + // this fix returns a per-atom vector (the effective electronegativity) + peratom_flag = 1; + size_peratom_cols = 0; + imax = 200; maxwarn = 1; @@ -312,6 +318,7 @@ void FixQtpieReaxFF::allocate_storage() memory->create(Hdia_inv,nmax,"qtpie:Hdia_inv"); memory->create(b_s,nmax,"qtpie:b_s"); memory->create(chi_eff,nmax,"qtpie:chi_eff"); + vector_atom = chi_eff; memory->create(b_t,nmax,"qtpie:b_t"); memory->create(b_prc,nmax,"qtpie:b_prc"); memory->create(b_prm,nmax,"qtpie:b_prm"); From 5d0f1aeeafee9f831e4f67fba7ec978d960953d8 Mon Sep 17 00:00:00 2001 From: Navraj Lalli Date: Thu, 26 Sep 2024 17:21:16 +0100 Subject: [PATCH 32/40] Expand documentation --- doc/src/fix_qtpie_reaxff.rst | 31 ++++++++++++++++----- doc/utils/sphinx-config/false_positives.txt | 4 +++ 2 files changed, 28 insertions(+), 7 deletions(-) diff --git a/doc/src/fix_qtpie_reaxff.rst b/doc/src/fix_qtpie_reaxff.rst index e0c2de1432..2d2b183491 100644 --- a/doc/src/fix_qtpie_reaxff.rst +++ b/doc/src/fix_qtpie_reaxff.rst @@ -38,21 +38,33 @@ Description The QTPIE charge equilibration method is an extension of the QEq charge equilibration method. With QTPIE, the partial charges on individual atoms are computed by minimizing the electrostatic energy of the system in the -same way as the QEq method but where the Mulliken electronegativity, +same way as the QEq method but where the absolute electronegativity, :math:`\chi_i`, of each atom in the QEq charge equilibration scheme :ref:`(Rappe and Goddard) ` is replaced with an effective electronegativity given by :ref:`(Chen) ` .. math:: - \chi_{\mathrm{eff},i} = \frac{\sum_{j=1}^{N} (\chi_i - \chi_j + \phi_j - \phi_i) S_{ij}} + \chi_{\mathrm{eff},i} = \frac{\sum_{j=1}^{N} (\chi_i - \chi_j) S_{ij}} {\sum_{m=1}^{N}S_{im}}, which acts to penalize long-range charge transfer seen with the QEq charge equilibration scheme. In this equation, :math:`N` is the number of atoms in -the system, :math:`S_{ij}` is the overlap integral between atom :math:`i` -and atom :math:`j`, and :math:`\phi_i` and :math:`\phi_j` are the electric -potentials at the position of atom :math:`i` and :math:`j` due to -an external electric field, respectively. +the system and :math:`S_{ij}` is the overlap integral between atom :math:`i` +and atom :math:`j`. + +The effect of an external electric field can be incorporated into the QTPIE +method by modifying the absolute or effective electronegativities of each +atom :ref:`(Chen) `. This fix models the effect of an external +electric field by using the effective electronegativity given in +:ref:`(Gergs) `: + +.. math:: + \chi_{\mathrm{eff},i} = \frac{\sum_{j=1}^{N} (\chi_i - \chi_j + \phi_j - \phi_i) S_{ij}} + {\sum_{m=1}^{N}S_{im}}, + +where :math:`\phi_i` and :math:`\phi_j` are the electric +potentials at the positions of atom :math:`i` and :math:`j` +due to the external electric field. This fix is typically used in conjunction with the ReaxFF force field model as implemented in the :doc:`pair_style reaxff ` @@ -64,7 +76,7 @@ charge equilibration performed by `fix qtpie/reaxff`, which is the same as in To be explicit, this fix replaces :math:`\chi_k` of eq. 3 in :ref:`(Aktulga) ` with :math:`\chi_{\mathrm{eff},k}`. -This fix requires the Mulliken electronegativity, :math:`\chi`, in eV, the +This fix requires the absolute electronegativity, :math:`\chi`, in eV, the self-Coulomb potential, :math:`\eta`, in eV, and the shielded Coulomb constant, :math:`\gamma`, in :math:`\AA^{-1}`. If the *params* setting above is the word "reaxff", then these are extracted from the @@ -173,6 +185,11 @@ maxiter 200 **(Chen)** Chen, Jiahao. Theory and applications of fluctuating-charge models. University of Illinois at Urbana-Champaign, 2009. +.. _Gergs: + +**(Gergs)** Gergs, Dirkmann and Mussenbrock. +Journal of Applied Physics 123.24 (2018). + .. _qeq-Aktulga2: **(Aktulga)** Aktulga, Fogarty, Pandit, Grama, Parallel Computing, 38, diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index 65c1031fcf..12e952d551 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -814,6 +814,7 @@ dipoleflag dir Direc dirname +Dirkmann discoverable discretization discretized @@ -975,6 +976,7 @@ elaplong elastance Electroneg electronegative +electronegativities electronegativity electroneutral electroneutrality @@ -1291,6 +1293,7 @@ Geocomputing georg Georg Geotechnica +Gergs germain Germann Germano @@ -2390,6 +2393,7 @@ Murdick Murtola Murty Muser +Mussenbrock mutexes Muto muVT From 350551ecac5f07aef7a08b3cb24e3c8c8e64edaf Mon Sep 17 00:00:00 2001 From: Navraj Lalli Date: Thu, 26 Sep 2024 21:27:53 +0100 Subject: [PATCH 33/40] Fix whitespace --- doc/src/fix_qtpie_reaxff.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/src/fix_qtpie_reaxff.rst b/doc/src/fix_qtpie_reaxff.rst index 2d2b183491..0075c573e1 100644 --- a/doc/src/fix_qtpie_reaxff.rst +++ b/doc/src/fix_qtpie_reaxff.rst @@ -187,7 +187,7 @@ University of Illinois at Urbana-Champaign, 2009. .. _Gergs: -**(Gergs)** Gergs, Dirkmann and Mussenbrock. +**(Gergs)** Gergs, Dirkmann and Mussenbrock. Journal of Applied Physics 123.24 (2018). .. _qeq-Aktulga2: From d86de2862b66e2749b5c1592970ca60825e8db6d Mon Sep 17 00:00:00 2001 From: Navraj Lalli Date: Mon, 30 Sep 2024 12:10:33 +0100 Subject: [PATCH 34/40] Make signs consistent with efield = -grad(phi) --- doc/src/fix_qtpie_reaxff.rst | 2 +- src/REAXFF/fix_qtpie_reaxff.cpp | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/doc/src/fix_qtpie_reaxff.rst b/doc/src/fix_qtpie_reaxff.rst index 0075c573e1..b7faa772af 100644 --- a/doc/src/fix_qtpie_reaxff.rst +++ b/doc/src/fix_qtpie_reaxff.rst @@ -59,7 +59,7 @@ electric field by using the effective electronegativity given in :ref:`(Gergs) `: .. math:: - \chi_{\mathrm{eff},i} = \frac{\sum_{j=1}^{N} (\chi_i - \chi_j + \phi_j - \phi_i) S_{ij}} + \chi_{\mathrm{eff},i} = \frac{\sum_{j=1}^{N} (\chi_i - \chi_j + \phi_i - \phi_j) S_{ij}} {\sum_{m=1}^{N}S_{im}}, where :math:`\phi_i` and :math:`\phi_j` are the electric diff --git a/src/REAXFF/fix_qtpie_reaxff.cpp b/src/REAXFF/fix_qtpie_reaxff.cpp index 30a7b7f71f..946457a4da 100644 --- a/src/REAXFF/fix_qtpie_reaxff.cpp +++ b/src/REAXFF/fix_qtpie_reaxff.cpp @@ -1166,7 +1166,7 @@ void FixQtpieReaxFF::calc_chi_eff() chia = chi[type[i]]; if (efield) { if (efield->varflag != FixEfield::ATOM) { - phia = factor*(x[i][0]*efield->ex + x[i][1]*efield->ey + x[i][2]*efield->ez); + phia = -factor*(x[i][0]*efield->ex + x[i][1]*efield->ey + x[i][2]*efield->ez); } else { // atom-style potential from FixEfield phia = efield->efield[i][3]; } @@ -1189,11 +1189,11 @@ void FixQtpieReaxFF::calc_chi_eff() if (efield) { if (efield->varflag != FixEfield::ATOM) { - phib = factor*(x[j][0]*efield->ex + x[j][1]*efield->ey + x[j][2]*efield->ez); + 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 - chib + phib - phia) * overlap; + sum_n += (chia - chib + phia - phib) * overlap; } else { sum_n += (chia - chib) * overlap; } From d5f630db6c034fa568bd0e35b4f2ba24e2fc6744 Mon Sep 17 00:00:00 2001 From: Navraj Lalli Date: Mon, 30 Sep 2024 12:28:16 +0100 Subject: [PATCH 35/40] Fix sign used for atom-style potential A positive sign in front of the electric potential is consistent with E = -grad(electric potential). --- src/REAXFF/fix_qeq_reaxff.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/REAXFF/fix_qeq_reaxff.cpp b/src/REAXFF/fix_qeq_reaxff.cpp index adaf5be031..37e90f582e 100644 --- a/src/REAXFF/fix_qeq_reaxff.cpp +++ b/src/REAXFF/fix_qeq_reaxff.cpp @@ -1158,7 +1158,7 @@ void FixQEqReaxFF::get_chi_field() for (int i = 0; i < nlocal; i++) { if (mask[i] & efgroupbit) { if (region && !region->match(x[i][0],x[i][1],x[i][2])) continue; - chi_field[i] = -efield->efield[i][3]; + chi_field[i] = efield->efield[i][3]; } } } From a5ab8be0a24a18fd3114d21c510f8629cb692f2a Mon Sep 17 00:00:00 2001 From: Navraj Lalli Date: Fri, 25 Oct 2024 14:16:25 +0100 Subject: [PATCH 36/40] Clarify restriction on periodic cell dimensions --- doc/src/fix_acks2_reaxff.rst | 3 ++- doc/src/fix_qeq_reaxff.rst | 3 ++- doc/src/fix_qtpie_reaxff.rst | 6 ++++-- 3 files changed, 8 insertions(+), 4 deletions(-) diff --git a/doc/src/fix_acks2_reaxff.rst b/doc/src/fix_acks2_reaxff.rst index ebb1b48051..566e17a330 100644 --- a/doc/src/fix_acks2_reaxff.rst +++ b/doc/src/fix_acks2_reaxff.rst @@ -111,7 +111,8 @@ LAMMPS was built with that package. See the :doc:`Build package This fix does not correctly handle interactions involving multiple periodic images of the same atom. Hence, it should not be used for -periodic cell dimensions less than :math:`10~\AA`. +periodic cell dimensions smaller than the non-bonded cutoff radius, +which is typically :math:`10~\AA` for ReaxFF simulations. This fix may be used in combination with :doc:`fix efield ` and will apply the external electric field during charge equilibration, diff --git a/doc/src/fix_qeq_reaxff.rst b/doc/src/fix_qeq_reaxff.rst index e90842ea6a..c449c8cda9 100644 --- a/doc/src/fix_qeq_reaxff.rst +++ b/doc/src/fix_qeq_reaxff.rst @@ -124,7 +124,8 @@ LAMMPS was built with that package. See the :doc:`Build package This fix does not correctly handle interactions involving multiple periodic images of the same atom. Hence, it should not be used for -periodic cell dimensions less than 10 Angstroms. +periodic cell dimensions smaller than the non-bonded cutoff radius, +which is typically :math:`10~\AA` for ReaxFF simulations. This fix may be used in combination with :doc:`fix efield ` and will apply the external electric field during charge equilibration, diff --git a/doc/src/fix_qtpie_reaxff.rst b/doc/src/fix_qtpie_reaxff.rst index b7faa772af..cf59e4701a 100644 --- a/doc/src/fix_qtpie_reaxff.rst +++ b/doc/src/fix_qtpie_reaxff.rst @@ -150,7 +150,8 @@ LAMMPS was built with that package. See the :doc:`Build package This fix does not correctly handle interactions involving multiple periodic images of the same atom. Hence, it should not be used for -periodic cell dimensions less than 10 Angstroms. +periodic cell dimensions smaller than the non-bonded cutoff radius, +which is typically :math:`10~\AA` for ReaxFF simulations. This fix may be used in combination with :doc:`fix efield ` and will apply the external electric field during charge equilibration, @@ -166,7 +167,8 @@ the *potential* keyword for `fix efield`. Related commands """""""""""""""" -:doc:`pair_style reaxff `, :doc:`fix qeq/reaxff ` +:doc:`pair_style reaxff `, :doc:`fix qeq/reaxff `, +:doc:`fix acks2/reaxff ` Default """"""" From 1c48d201b41765a81a4820fbee38df09fd5b10ad Mon Sep 17 00:00:00 2001 From: Navraj Lalli Date: Fri, 8 Nov 2024 10:07:16 +0000 Subject: [PATCH 37/40] Remove unused pack_flag = 5 options --- src/REAXFF/fix_qtpie_reaxff.cpp | 44 +++------------------------------ 1 file changed, 3 insertions(+), 41 deletions(-) diff --git a/src/REAXFF/fix_qtpie_reaxff.cpp b/src/REAXFF/fix_qtpie_reaxff.cpp index 946457a4da..48c1109178 100644 --- a/src/REAXFF/fix_qtpie_reaxff.cpp +++ b/src/REAXFF/fix_qtpie_reaxff.cpp @@ -903,15 +903,6 @@ int FixQtpieReaxFF::pack_forward_comm(int n, int *list, double *buf, for (m = 0; m < n; m++) buf[m] = t[list[m]]; else if (pack_flag == 4) for (m = 0; m < n; m++) buf[m] = atom->q[list[m]]; - else if (pack_flag == 5) { - m = 0; - for (int i = 0; i < n; i++) { - int j = 2 * list[i]; - buf[m++] = d[j]; - buf[m++] = d[j+1]; - } - return m; - } return n; } @@ -929,15 +920,6 @@ void FixQtpieReaxFF::unpack_forward_comm(int n, int first, double *buf) for (m = 0, i = first; m < n; m++, i++) t[i] = buf[m]; else if (pack_flag == 4) for (m = 0, i = first; m < n; m++, i++) atom->q[i] = buf[m]; - else if (pack_flag == 5) { - int last = first + n; - m = 0; - for (i = first; i < last; i++) { - int j = 2 * i; - d[j] = buf[m++]; - d[j+1] = buf[m++]; - } - } } /* ---------------------------------------------------------------------- */ @@ -945,35 +927,15 @@ void FixQtpieReaxFF::unpack_forward_comm(int n, int first, double *buf) int FixQtpieReaxFF::pack_reverse_comm(int n, int first, double *buf) { int i, m; - if (pack_flag == 5) { - m = 0; - int last = first + n; - for (i = first; i < last; i++) { - int indxI = 2 * i; - buf[m++] = q[indxI]; - buf[m++] = q[indxI+1]; - } - return m; - } else { - for (m = 0, i = first; m < n; m++, i++) buf[m] = q[i]; - return n; - } + for (m = 0, i = first; m < n; m++, i++) buf[m] = q[i]; + return n; } /* ---------------------------------------------------------------------- */ void FixQtpieReaxFF::unpack_reverse_comm(int n, int *list, double *buf) { - if (pack_flag == 5) { - int m = 0; - for (int i = 0; i < n; i++) { - int indxI = 2 * list[i]; - q[indxI] += buf[m++]; - q[indxI+1] += buf[m++]; - } - } else { - for (int m = 0; m < n; m++) q[list[m]] += buf[m]; - } + for (int m = 0; m < n; m++) q[list[m]] += buf[m]; } /* ---------------------------------------------------------------------- From 95899b53b86b68249dd26ee4c32b18bc0c8ab398 Mon Sep 17 00:00:00 2001 From: Navraj Lalli Date: Fri, 8 Nov 2024 10:26:47 +0000 Subject: [PATCH 38/40] Add fix qtpie/reaxff to pair_style reaxff docs --- doc/src/pair_reaxff.rst | 26 ++++++++++++++------------ 1 file changed, 14 insertions(+), 12 deletions(-) diff --git a/doc/src/pair_reaxff.rst b/doc/src/pair_reaxff.rst index 03d53d1ff4..84403c15d1 100644 --- a/doc/src/pair_reaxff.rst +++ b/doc/src/pair_reaxff.rst @@ -20,7 +20,7 @@ Syntax .. parsed-literal:: keyword = *checkqeq* or *lgvdw* or *safezone* or *mincap* or *minhbonds* or *tabulate* or *list/blocking* - *checkqeq* value = *yes* or *no* = whether or not to require qeq/reaxff or acks2/reaxff fix + *checkqeq* value = *yes* or *no* = whether or not to require one of fix qeq/reaxff, fix acks2/reaxff or fix qtpie/reaxff *enobonds* value = *yes* or *no* = whether or not to tally energy of atoms with no bonds *lgvdw* value = *yes* or *no* = whether or not to use a low gradient vdW correction *safezone* = factor used for array allocation @@ -120,20 +120,22 @@ up that process. The ReaxFF parameter files provided were created using a charge equilibration (QEq) model for handling the electrostatic interactions. -Therefore, by default, LAMMPS requires that either the -:doc:`fix qeq/reaxff ` or the -:doc:`fix qeq/shielded ` or :doc:`fix acks2/reaxff ` -command be used with -*pair_style reaxff* when simulating a ReaxFF model, to equilibrate -the charges each timestep. +Therefore, by default, LAMMPS requires that +:doc:`fix qeq/reaxff ` or :doc:`fix qeq/shielded ` +or :doc:`fix acks2/reaxff ` +or :doc:`fix qtpie/reaxff ` +is used with *pair_style reaxff* when simulating a ReaxFF model, +to equilibrate the charges at each timestep. +See the :doc:`fix qeq/reaxff ` or :doc:`fix qeq/shielded ` +or :doc:`fix acks2/reaxff ` +or :doc:`fix qtpie/reaxff ` +command documentation for more details. Using the keyword *checkqeq* with the value *no* turns off the check for the QEq fixes, allowing a simulation to be run without charge equilibration. In this case, the static charges you assign to each atom will be used for computing the electrostatic interactions in -the system. See the :doc:`fix qeq/reaxff ` or -:doc:`fix qeq/shielded ` or :doc:`fix acks2/reaxff ` -command documentation for more details. +the system. Using the optional keyword *lgvdw* with the value *yes* turns on the low-gradient correction of ReaxFF for long-range London Dispersion, @@ -372,8 +374,8 @@ Related commands """""""""""""""" :doc:`pair_coeff `, :doc:`fix qeq/reaxff `, -:doc:`fix acks2/reaxff `, :doc:`fix reaxff/bonds `, -:doc:`fix reaxff/species `, +:doc:`fix acks2/reaxff `, :doc:`fix qtpie/reaxff `, +:doc:`fix reaxff/bonds `, :doc:`fix reaxff/species `, :doc:`compute reaxff/atom ` Default From e84c45c6e7d61a57bcdba143435fda23a0d6d67d Mon Sep 17 00:00:00 2001 From: Navraj Lalli Date: Fri, 8 Nov 2024 10:38:49 +0000 Subject: [PATCH 39/40] Fix whitespace --- doc/src/pair_reaxff.rst | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/doc/src/pair_reaxff.rst b/doc/src/pair_reaxff.rst index 84403c15d1..495572dc0e 100644 --- a/doc/src/pair_reaxff.rst +++ b/doc/src/pair_reaxff.rst @@ -124,9 +124,9 @@ Therefore, by default, LAMMPS requires that :doc:`fix qeq/reaxff ` or :doc:`fix qeq/shielded ` or :doc:`fix acks2/reaxff ` or :doc:`fix qtpie/reaxff ` -is used with *pair_style reaxff* when simulating a ReaxFF model, -to equilibrate the charges at each timestep. -See the :doc:`fix qeq/reaxff ` or :doc:`fix qeq/shielded ` +is used with *pair_style reaxff* when simulating a ReaxFF model, +to equilibrate the charges at each timestep. +See the :doc:`fix qeq/reaxff ` or :doc:`fix qeq/shielded ` or :doc:`fix acks2/reaxff ` or :doc:`fix qtpie/reaxff ` command documentation for more details. @@ -135,7 +135,7 @@ Using the keyword *checkqeq* with the value *no* turns off the check for the QEq fixes, allowing a simulation to be run without charge equilibration. In this case, the static charges you assign to each atom will be used for computing the electrostatic interactions in -the system. +the system. Using the optional keyword *lgvdw* with the value *yes* turns on the low-gradient correction of ReaxFF for long-range London Dispersion, From 5673375d2153c85e24ebfc371c0c6a603cfa15e3 Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Tue, 12 Nov 2024 12:32:37 -0700 Subject: [PATCH 40/40] Add more related commands to docs --- doc/src/fix_acks2_reaxff.rst | 3 ++- doc/src/fix_qeq_reaxff.rst | 3 ++- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/doc/src/fix_acks2_reaxff.rst b/doc/src/fix_acks2_reaxff.rst index 566e17a330..79a9cf8ea6 100644 --- a/doc/src/fix_acks2_reaxff.rst +++ b/doc/src/fix_acks2_reaxff.rst @@ -123,7 +123,8 @@ components in non-periodic directions. Related commands """""""""""""""" -:doc:`pair_style reaxff `, :doc:`fix qeq/reaxff ` +:doc:`pair_style reaxff `, :doc:`fix qeq/reaxff `, +:doc:`fix qtpi/reaxff ` Default """"""" diff --git a/doc/src/fix_qeq_reaxff.rst b/doc/src/fix_qeq_reaxff.rst index c449c8cda9..e1a09c4fc3 100644 --- a/doc/src/fix_qeq_reaxff.rst +++ b/doc/src/fix_qeq_reaxff.rst @@ -139,7 +139,8 @@ as an atom-style variable using the *potential* keyword for `fix efield`. Related commands """""""""""""""" -:doc:`pair_style reaxff `, :doc:`fix qeq/shielded ` +:doc:`pair_style reaxff `, :doc:`fix qeq/shielded `, +:doc:`fix acks2/reaxff `, :doc:`fix qtpie/reaxff ` Default """""""