diff --git a/src/.gitignore b/src/.gitignore index a2f7da0400..cdd3c99603 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -996,6 +996,8 @@ /fix_qeq_fire.h /fix_qeq_reaxff.cpp /fix_qeq_reaxff.h +/fix_qeqr_reaxff.cpp +/fix_qeqr_reaxff.h /fix_qmmm.cpp /fix_qmmm.h /fix_qtpie_reaxff.cpp diff --git a/src/OPENMP/pair_reaxff_omp.cpp b/src/OPENMP/pair_reaxff_omp.cpp index 85369cc7bf..df28b64e98 100644 --- a/src/OPENMP/pair_reaxff_omp.cpp +++ b/src/OPENMP/pair_reaxff_omp.cpp @@ -106,6 +106,7 @@ void PairReaxFFOMP::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("^qeqr/reax").size() + modify->get_fix_by_style("^qeq/shielded").size() + acks2_fixes.size() + modify->get_fix_by_style("^qtpie/reax").size(); diff --git a/src/REAXFF/fix_qeqr_reaxff.cpp b/src/REAXFF/fix_qeqr_reaxff.cpp new file mode 100644 index 0000000000..98f13da2f0 --- /dev/null +++ b/src/REAXFF/fix_qeqr_reaxff.cpp @@ -0,0 +1,125 @@ +// 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 authors: + Navraj S Lalli, Imperial College London (navrajsinghlalli@gmail.com) + +------------------------------------------------------------------------- */ + +#include "fix_qeqr_reaxff.h" + +#include "atom.h" +#include "comm.h" +#include "domain.h" +#include "error.h" +#include "fix_efield.h" +#include "force.h" +#include "group.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 + +using namespace LAMMPS_NS; +using namespace FixConst; + +static constexpr double ANGSTROM_TO_BOHRRADIUS = 1.8897261259; + +/* ---------------------------------------------------------------------- */ + +FixQEqrReaxFF::FixQEqrReaxFF(LAMMPS *lmp, int narg, char **arg) : + FixQtpieReaxFF(lmp, narg, arg) +{ +} + +/* ---------------------------------------------------------------------- */ + +void FixQEqrReaxFF::calc_chi_eff() +{ + memset(&chi_eff[0],0,atom->nmax*sizeof(double)); + + const auto x = (const double * const *)atom->x; + const int *type = atom->type; + + double dist,overlap,sum_n,sum_d,expa,expb,chia,phia,phib,p,m; + 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; + + if (efield) { + if (efield->varflag != FixEfield::CONSTANT) + efield->update_efield_variables(); + + // compute chi_eff for each local atom + for (i = 0; i < nn; i++) { + expa = gauss_exp[type[i]]; + chia = chi[type[i]]; + if (efield->varflag != FixEfield::ATOM) { + phia = -factor*(x[i][0]*efield->ex + x[i][1]*efield->ey + x[i][2]*efield->ez); + } else { // atom-style potential from FixEfield + phia = efield->efield[i][3]; + } + + sum_n = 0.0; + sum_d = 0.0; + + for (j = 0; j < nt; j++) { + dist = distance(x[i],x[j])*ANGSTROM_TO_BOHRRADIUS; // in atomic units + + if (dist < dist_cutoff) { + expb = gauss_exp[type[j]]; + + // overlap integral of two normalised 1s Gaussian type orbitals + p = expa + expb; + m = expa * expb / p; + overlap = pow((4.0*m/p),0.75) * exp(-m*dist*dist); + + if (efield->varflag != FixEfield::ATOM) { + phib = -factor*(x[j][0]*efield->ex + x[j][1]*efield->ey + x[j][2]*efield->ez); + } else { // atom-style potential from FixEfield + phib = efield->efield[j][3]; + } + sum_n += (chia + scale * (phia - phib)) * overlap; + sum_d += overlap; + } + } + chi_eff[i] = sum_n / sum_d; + } + } else { + for (i = 0; i < nn; i++) { + chi_eff[i] = chi[type[i]]; + } + } +} diff --git a/src/REAXFF/fix_qeqr_reaxff.h b/src/REAXFF/fix_qeqr_reaxff.h new file mode 100644 index 0000000000..a48a6ee632 --- /dev/null +++ b/src/REAXFF/fix_qeqr_reaxff.h @@ -0,0 +1,38 @@ +/* -*- 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(qeqr/reaxff,FixQEqrReaxFF); +// clang-format on +#else + +#ifndef LMP_FIX_QEQR_REAXFF_H +#define LMP_FIX_QEQR_REAXFF_H + +#include "fix_qtpie_reaxff.h" + +namespace LAMMPS_NS { + +class FixQEqrReaxFF : public FixQtpieReaxFF { + public: + FixQEqrReaxFF(class LAMMPS *, int, char **); + + protected: + void calc_chi_eff() override; +}; + +} // namespace LAMMPS_NS + +#endif +#endif diff --git a/src/REAXFF/fix_qtpie_reaxff.cpp b/src/REAXFF/fix_qtpie_reaxff.cpp index 3d5289ac6d..86331b2335 100644 --- a/src/REAXFF/fix_qtpie_reaxff.cpp +++ b/src/REAXFF/fix_qtpie_reaxff.cpp @@ -81,8 +81,9 @@ FixQtpieReaxFF::FixQtpieReaxFF(LAMMPS *lmp, int narg, char **arg) : imax = 200; maxwarn = 1; + scale = 1.0; - if ((narg < 9) || (narg > 12)) error->all(FLERR,"Illegal fix {} command", style); + if ((narg < 9) || (narg > 14)) 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); @@ -101,6 +102,11 @@ FixQtpieReaxFF::FixQtpieReaxFF(LAMMPS *lmp, int narg, char **arg) : error->all(FLERR,"Illegal fix {} command", style); imax = utils::numeric(FLERR,arg[iarg+1],false,lmp); iarg++; + } else if (strcmp(arg[iarg],"scale") == 0) { + if (iarg+1 > narg-1) + error->all(FLERR,"Illegal fix {} command", style); + scale = utils::numeric(FLERR,arg[iarg+1],false,lmp); + iarg++; } else error->all(FLERR,"Illegal fix {} command", style); iarg++; } @@ -1172,7 +1178,7 @@ void FixQtpieReaxFF::calc_chi_eff() } else { // atom-style potential from FixEfield phib = efield->efield[j][3]; } - sum_n += (chia - chib + phia - phib) * overlap; + sum_n += (chia - chib + scale * (phia - phib)) * overlap; } else { sum_n += (chia - chib) * overlap; } diff --git a/src/REAXFF/fix_qtpie_reaxff.h b/src/REAXFF/fix_qtpie_reaxff.h index 2f86e27a7a..a0724912ff 100644 --- a/src/REAXFF/fix_qtpie_reaxff.h +++ b/src/REAXFF/fix_qtpie_reaxff.h @@ -91,6 +91,7 @@ class FixQtpieReaxFF : public Fix { char *gauss_file; // input file for gaussian orbital exponents double *gauss_exp; // array of gaussian orbital exponents for each atom type double dist_cutoff; // separation distance beyond which to neglect overlap integrals + double scale; // scaling factor for electric polarization effects void pertype_parameters(char *); void init_shielding(); @@ -128,7 +129,7 @@ class FixQtpieReaxFF : public Fix { void vector_sum(double *, double, double *, double, double *, int); void vector_add(double *, double, double *, int); - void calc_chi_eff(); + virtual void calc_chi_eff(); double find_min_exp(const double*, const int); double distance(const double*, const double*); diff --git a/src/REAXFF/pair_reaxff.cpp b/src/REAXFF/pair_reaxff.cpp index b36d9e79aa..56a2535108 100644 --- a/src/REAXFF/pair_reaxff.cpp +++ b/src/REAXFF/pair_reaxff.cpp @@ -341,12 +341,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("^qeqr/reax").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 or " - "fix qtpie/reaxff commands"); + "fix qtpie/reaxff commands or fix qeqr/reaxff"); 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 108395cc2c..83b5cc8d7e 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 FixQEqrReaxFF; friend class FixQtpieReaxFF; public: