diff --git a/src/KOKKOS/fix_acks2_reaxff_kokkos.cpp b/src/KOKKOS/fix_acks2_reaxff_kokkos.cpp index d02cd1f255..cb65667a17 100644 --- a/src/KOKKOS/fix_acks2_reaxff_kokkos.cpp +++ b/src/KOKKOS/fix_acks2_reaxff_kokkos.cpp @@ -393,8 +393,7 @@ void FixACKS2ReaxFFKokkos::pre_force(int vflag) k_X_diag.template sync(); } - if (field_flag) - get_chi_field(); + if (efield) get_chi_field(); // init_matvec @@ -521,8 +520,7 @@ void FixACKS2ReaxFFKokkos::allocate_array() d_z = k_z.template view(); } - if (field_flag) - get_chi_field(); + if (efield) get_chi_field(); // init_storage Kokkos::parallel_for(Kokkos::RangePolicy(0,NN),*this); diff --git a/src/KOKKOS/fix_qeq_reaxff_kokkos.cpp b/src/KOKKOS/fix_qeq_reaxff_kokkos.cpp index af67d25af8..b35bbc46ee 100644 --- a/src/KOKKOS/fix_qeq_reaxff_kokkos.cpp +++ b/src/KOKKOS/fix_qeq_reaxff_kokkos.cpp @@ -257,8 +257,7 @@ void FixQEqReaxFFKokkos::pre_force(int /*vflag*/) // init_matvec - if (field_flag) - get_chi_field(); + if (efield) get_chi_field(); k_s_hist.template sync(); k_t_hist.template sync(); @@ -384,8 +383,7 @@ void FixQEqReaxFFKokkos::allocate_array() // init_storage - if (field_flag) - get_chi_field(); + if (efield) get_chi_field(); FixQEqReaxFFKokkosZeroFunctor zero_functor(this); Kokkos::parallel_for(ignum,zero_functor); diff --git a/src/OPENMP/fix_qeq_reaxff_omp.cpp b/src/OPENMP/fix_qeq_reaxff_omp.cpp index 2071755fd1..0b0ba589cf 100644 --- a/src/OPENMP/fix_qeq_reaxff_omp.cpp +++ b/src/OPENMP/fix_qeq_reaxff_omp.cpp @@ -232,8 +232,7 @@ void FixQEqReaxFFOMP::compute_H() void FixQEqReaxFFOMP::init_storage() { - if (field_flag) - get_chi_field(); + if (efield) get_chi_field(); #if defined(_OPENMP) #pragma omp parallel for schedule(static) @@ -241,8 +240,7 @@ void FixQEqReaxFFOMP::init_storage() for (int i = 0; i < NN; i++) { Hdia_inv[i] = 1. / eta[atom->type[i]]; b_s[i] = -chi[atom->type[i]]; - if (field_flag) - b_s[i] -= chi_field[i]; + if (efield) b_s[i] -= chi_field[i]; b_t[i] = -1.0; b_prc[i] = 0; b_prm[i] = 0; @@ -279,8 +277,7 @@ void FixQEqReaxFFOMP::pre_force(int /* vflag */) if (n > n_cap*DANGER_ZONE || m_fill > m_cap*DANGER_ZONE) reallocate_matrix(); - if (field_flag) - get_chi_field(); + if (efield) get_chi_field(); init_matvec(); @@ -318,8 +315,7 @@ void FixQEqReaxFFOMP::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 (field_flag) - b_s[i] -= chi_field[i]; + if (efield) b_s[i] -= chi_field[i]; b_t[i] = -1.0; // Predictor Step @@ -348,8 +344,7 @@ void FixQEqReaxFFOMP::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 (field_flag) - b_s[i] -= chi_field[i]; + if (efield) b_s[i] -= chi_field[i]; b_t[i] = -1.0; /* linear extrapolation for s & t from previous solutions */ diff --git a/src/REAXFF/fix_acks2_reaxff.cpp b/src/REAXFF/fix_acks2_reaxff.cpp index aa4923d91d..1de4eaaf1c 100644 --- a/src/REAXFF/fix_acks2_reaxff.cpp +++ b/src/REAXFF/fix_acks2_reaxff.cpp @@ -316,15 +316,13 @@ void FixACKS2ReaxFF::init_bondcut() void FixACKS2ReaxFF::init_storage() { - if (field_flag) - get_chi_field(); + if (efield) get_chi_field(); for (int ii = 0; ii < NN; ii++) { int i = ilist[ii]; if (atom->mask[i] & groupbit) { b_s[i] = -chi[atom->type[i]]; - if (field_flag) - b_s[i] -= chi_field[i]; + if (efield) b_s[i] -= chi_field[i]; b_s[NN + i] = 0.0; s[i] = 0.0; s[NN + i] = 0.0; @@ -366,8 +364,7 @@ void FixACKS2ReaxFF::pre_force(int /*vflag*/) if (n > n_cap*DANGER_ZONE || m_fill > m_cap*DANGER_ZONE) reallocate_matrix(); - if (field_flag) - get_chi_field(); + if (efield) get_chi_field(); init_matvec(); @@ -402,8 +399,7 @@ void FixACKS2ReaxFF::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 (field_flag) - b_s[i] -= chi_field[i]; + if (efield) b_s[i] -= chi_field[i]; b_s[NN+i] = 0.0; /* cubic extrapolation for s from previous solutions */ diff --git a/src/REAXFF/fix_qeq_reaxff.cpp b/src/REAXFF/fix_qeq_reaxff.cpp index ac143d95c4..5b5a792c16 100644 --- a/src/REAXFF/fix_qeq_reaxff.cpp +++ b/src/REAXFF/fix_qeq_reaxff.cpp @@ -33,6 +33,7 @@ #include "neigh_request.h" #include "neighbor.h" #include "pair.h" +#include "region.h" #include "respa.h" #include "text_file_reader.h" #include "update.h" @@ -382,11 +383,14 @@ void FixQEqReaxFF::init() if (group->count(igroup) == 0) error->all(FLERR,"Fix qeq/reaxff group has no atoms"); - field_flag = 0; - for (int n = 0; n < modify->nfix; n++) - if (utils::strmatch(modify->fix[n]->style,"^efield")) - field_flag = 1; + efield = nullptr; + int ifix = modify->find_fix_by_style("^efield"); + if (ifix >= 0) efield = (FixEfield *) modify->fix[ifix]; + if (efield && (strcmp(update->unit_style,"real") != 0)) + error->all(FLERR,"Must use unit_style real with fix qeq/reax and external fields"); + if (efield && efield->varflag != FixEfield::CONSTANT) + error->all(FLERR,"Cannot yet use fix qeq/reaxff with variable efield"); // need a half neighbor list w/ Newton off and ghost neighbors // built whenever re-neighboring occurs @@ -511,16 +515,14 @@ void FixQEqReaxFF::min_setup_pre_force(int vflag) void FixQEqReaxFF::init_storage() { - if (field_flag) - get_chi_field(); + 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 (field_flag) - b_s[i] -= chi_field[i]; + if (efield) b_s[i] -= chi_field[i]; b_t[i] = -1.0; b_prc[i] = 0; b_prm[i] = 0; @@ -558,8 +560,7 @@ void FixQEqReaxFF::pre_force(int /*vflag*/) if (n > n_cap*DANGER_ZONE || m_fill > m_cap*DANGER_ZONE) reallocate_matrix(); - if (field_flag) - get_chi_field(); + if (efield) get_chi_field(); init_matvec(); @@ -600,8 +601,7 @@ void FixQEqReaxFF::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 (field_flag) - b_s[i] -= chi_field[i]; + if (efield) b_s[i] -= chi_field[i]; b_t[i] = -1.0; /* quadratic extrapolation for s & t from previous solutions */ @@ -1072,25 +1072,42 @@ void FixQEqReaxFF::vector_add(double* dest, double c, double* v, int k) void FixQEqReaxFF::get_chi_field() { - int nlocal = atom->nlocal; - memset(&chi_field[0],0.0,atom->nmax*sizeof(double)); + if (!efield) return; - if (!(strcmp(update->unit_style,"real") == 0)) - error->all(FLERR,"Must use unit_style real with fix qeq/reax and external fields"); + const double * const *x = (const double * const *)atom->x; + const int *mask = atom->mask; + const imageint *image = atom->image; - double factor = 1.0/force->qe2f; + const double factor = 1.0/force->qe2f; + const int nlocal = atom->nlocal; - // loop over all fixes, find fix efield + // update electric field region if necessary - for (int n = 0; n < modify->nfix; n++) { - if (utils::strmatch(modify->fix[n]->style,"^efield")) { + Region *region = nullptr; + if (efield->iregion >= 0) { + region = domain->regions[efield->iregion]; + region->prematch(); + } - FixEfield* fix_efield = (FixEfield*) modify->fix[n]; - double* field_energy = fix_efield->get_energy(); // Real units of kcal/mol/angstrom, need to convert to eV + // we currently only constant efield. Also atom selection is for the group of fix efield. - for (int i = 0; i < nlocal; i++) - chi_field[i] += field_energy[i]*factor; + if (efield->varflag == FixEfield::CONSTANT) { + double unwrap[3]; + const double fx = efield->ex; + const double fy = efield->ey; + const double fz = efield->ez; + const int efgroupbit = efield->groupbit; + + // charge interactions + // force = qE, potential energy = F dot x in unwrapped coords + + 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*(fx*unwrap[0] + fy*unwrap[1] + fz*unwrap[2]); + } } } } diff --git a/src/REAXFF/fix_qeq_reaxff.h b/src/REAXFF/fix_qeq_reaxff.h index cc5bd09791..ac54ff72ca 100644 --- a/src/REAXFF/fix_qeq_reaxff.h +++ b/src/REAXFF/fix_qeq_reaxff.h @@ -65,6 +65,7 @@ class FixQEqReaxFF : public Fix { 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 @@ -140,8 +141,6 @@ class FixQEqReaxFF : public Fix { // 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 - - int field_flag; // 1: field enabled }; } // namespace LAMMPS_NS diff --git a/src/fix_efield.cpp b/src/fix_efield.cpp index 229e63b74d..c6d8c36ada 100644 --- a/src/fix_efield.cpp +++ b/src/fix_efield.cpp @@ -19,30 +19,29 @@ #include "fix_efield.h" -#include #include "atom.h" -#include "update.h" -#include "domain.h" #include "comm.h" -#include "modify.h" -#include "force.h" -#include "respa.h" -#include "input.h" -#include "variable.h" -#include "region.h" -#include "memory.h" +#include "domain.h" #include "error.h" +#include "force.h" +#include "input.h" +#include "memory.h" +#include "modify.h" +#include "region.h" +#include "respa.h" +#include "update.h" +#include "variable.h" + +#include using namespace LAMMPS_NS; using namespace FixConst; -enum{NONE,CONSTANT,EQUAL,ATOM}; - /* ---------------------------------------------------------------------- */ FixEfield::FixEfield(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), xstr(nullptr), ystr(nullptr), zstr(nullptr), - estr(nullptr), idregion(nullptr), efield(nullptr), energy(nullptr) + estr(nullptr), idregion(nullptr), efield(nullptr) { if (narg < 6) error->all(FLERR,"Illegal fix efield command"); @@ -125,7 +124,6 @@ FixEfield::~FixEfield() delete [] estr; delete [] idregion; memory->destroy(efield); - memory->destroy(energy); } /* ---------------------------------------------------------------------- */ @@ -453,69 +451,3 @@ double FixEfield::compute_vector(int n) return fsum_all[n+1]; } -/* ---------------------------------------------------------------------- - get E -------------------------------------------------------------------------- */ - -double* FixEfield::get_energy() -{ - double **x = atom->x; - int *mask = atom->mask; - imageint *image = atom->image; - - int nlocal = atom->nlocal; - - // reallocate energy array if necessary - - if (atom->nmax > maxatom_energy) { - maxatom_energy = atom->nmax; - memory->destroy(energy); - memory->create(energy,maxatom_energy,"efield:energy"); - } - memset(&energy[0],0.0,maxatom_energy*sizeof(double)); - - // update region if necessary - - Region *region = NULL; - if (iregion >= 0) { - region = domain->regions[iregion]; - region->prematch(); - } - - int warn_flag_local = 0; - - // constant efield - - if (varflag == CONSTANT) { - double unwrap[3]; - - // charge interactions - // force = qE, potential energy = F dot x in unwrapped coords - - if (qflag) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - if (region && !region->match(x[i][0],x[i][1],x[i][2])) continue; - const double fx = ex; - const double fy = ey; - const double fz = ez; - - domain->unmap(x[i],image[i],unwrap); - energy[i] -= fx*unwrap[0] + fy*unwrap[1] + fz*unwrap[2]; - if (fabs(fx*(x[i][0]-unwrap[0])) + fabs(fy*(x[i][1]-unwrap[1])) + - fabs(fz*(x[i][2]-unwrap[2])) > 0.0) - warn_flag_local = 1; - } - } - } - } else { - error->all(FLERR,"Cannot yet use fix qeq/reaxff with variable efield"); - } - - int warn_flag; - MPI_Allreduce(&warn_flag_local,&warn_flag,1,MPI_INT,MPI_SUM,world); - if (warn_flag && comm->me == 0) - error->warning(FLERR,"Using non-zero image flags in field direction with fix qeq/reaxff"); - - return energy; -} diff --git a/src/fix_efield.h b/src/fix_efield.h index 86ccdae7e5..20201518ab 100644 --- a/src/fix_efield.h +++ b/src/fix_efield.h @@ -25,6 +25,7 @@ FixStyle(efield,FixEfield); namespace LAMMPS_NS { class FixEfield : public Fix { + friend class FixQEqReaxFF; public: FixEfield(class LAMMPS *, int, char **); ~FixEfield(); @@ -38,9 +39,10 @@ class FixEfield : public Fix { double memory_usage(); double compute_scalar(); double compute_vector(int); - double* get_energy(); - private: + enum { NONE, CONSTANT, EQUAL, ATOM }; + + protected: double ex, ey, ez; int varflag, iregion; char *xstr, *ystr, *zstr, *estr; @@ -50,8 +52,8 @@ class FixEfield : public Fix { double qe2f; int qflag, muflag; - int maxatom,maxatom_energy; - double **efield,*energy; + int maxatom, maxatom_energy; + double **efield; int force_flag; double fsum[4], fsum_all[4];