ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

move chi field calculation to fix qeq/reaxff

Browse Source
This commit is contained in:
Axel Kohlmeyer
2021-10-11 14:35:23 -04:00
parent 70cbb72e42
commit 6d2b32f0b2
8 changed files with 73 additions and 136 deletions
Download Patch File Download Diff File Expand all files Collapse all files

6
src/KOKKOS/fix_acks2_reaxff_kokkos.cpp
View File

@ -393,8 +393,7 @@ void FixACKS2ReaxFFKokkos<DeviceType>::pre_force(int vflag)
k_X_diag.template sync<DeviceType>();
}
if (field_flag)
get_chi_field();
if (efield) get_chi_field();
// init_matvec
@ -521,8 +520,7 @@ void FixACKS2ReaxFFKokkos<DeviceType>::allocate_array()
d_z = k_z.template view<DeviceType>();
}
if (field_flag)
get_chi_field();
if (efield) get_chi_field();
// init_storage
Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType,TagACKS2Zero>(0,NN),*this);

6
src/KOKKOS/fix_qeq_reaxff_kokkos.cpp
View File

@ -257,8 +257,7 @@ void FixQEqReaxFFKokkos<DeviceType>::pre_force(int /*vflag*/)
// init_matvec
if (field_flag)
get_chi_field();
if (efield) get_chi_field();
k_s_hist.template sync<DeviceType>();
k_t_hist.template sync<DeviceType>();
@ -384,8 +383,7 @@ void FixQEqReaxFFKokkos<DeviceType>::allocate_array()
// init_storage
if (field_flag)
get_chi_field();
if (efield) get_chi_field();
FixQEqReaxFFKokkosZeroFunctor<DeviceType> zero_functor(this);
Kokkos::parallel_for(ignum,zero_functor);

15
src/OPENMP/fix_qeq_reaxff_omp.cpp
View File

@ -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 */

12
src/REAXFF/fix_acks2_reaxff.cpp
View File

@ -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 */

65
src/REAXFF/fix_qeq_reaxff.cpp
View File

@ -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]);
}
}
}
}

3
src/REAXFF/fix_qeq_reaxff.h
View File

@ -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

92
src/fix_efield.cpp
View File

@ -19,30 +19,29 @@
#include "fix_efield.h"
#include <cstring>
#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 <cstring>
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;
}

10
src/fix_efield.h
View File

@ -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];