Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F85654236
pair_reax.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Mon, Sep 30, 19:11
Size
27 KB
Mime Type
text/x-c
Expires
Wed, Oct 2, 19:11 (1 d, 21 h)
Engine
blob
Format
Raw Data
Handle
21231744
Attached To
rLAMMPS lammps
pair_reax.cpp
View Options
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
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: Aidan Thompson (Sandia, athomps@sandia.gov)
Hansohl Cho (MIT, hansohl@mit.edu)
LAMMPS implementation of the Reactive Force Field (ReaxFF) is based on
Aidan Thompson's GRASP code
(General Reactive Atomistic Simulation Program)
and Ardi Van Duin's original ReaxFF code
------------------------------------------------------------------------- */
#include "mpi.h"
#include "math.h"
#include "stdio.h"
#include "stdlib.h"
#include "string.h"
#include "pair_reax.h"
#include "pair_reax_fortran.h"
#include "atom.h"
#include "update.h"
#include "force.h"
#include "comm.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
#define SMALL 0.0001
/* ---------------------------------------------------------------------- */
PairREAX::PairREAX(LAMMPS *lmp) : Pair(lmp)
{
single_enable = 0;
restartinfo = 0;
one_coeff = 1;
no_virial_fdotr_compute = 1;
nextra = 14;
pvector = new double[nextra];
cutmax = 0.0;
hbcut = 6.0;
ihbnew = 1;
itripstaball = 1;
iprune = 4;
ihb = 1;
chpot = 0;
nmax = 0;
arow_ptr = NULL;
ch = NULL;
elcvec = NULL;
rcg = NULL;
wcg = NULL;
pcg = NULL;
poldcg = NULL;
qcg = NULL;
matmax = 0;
aval = NULL;
acol_ind = NULL;
comm_forward = 1;
comm_reverse = 1;
precision = 1.0e-6;
}
/* ----------------------------------------------------------------------
free all arrays
check if allocated, since class can be destructed when incomplete
------------------------------------------------------------------------- */
PairREAX::~PairREAX()
{
delete [] pvector;
if (allocated) {
memory->destroy(setflag);
memory->destroy(cutsq);
for (int i = 1; i <= atom->ntypes; i++)
delete [] param_list[i].params;
delete [] param_list;
delete [] map;
}
memory->destroy(arow_ptr);
memory->destroy(ch);
memory->destroy(elcvec);
memory->destroy(rcg);
memory->destroy(wcg);
memory->destroy(pcg);
memory->destroy(poldcg);
memory->destroy(qcg);
memory->destroy(aval);
memory->destroy(acol_ind);
}
/* ---------------------------------------------------------------------- */
void PairREAX::compute(int eflag, int vflag)
{
int i,j;
double evdwl,ecoul;
double energy_charge_equilibration;
evdwl = ecoul = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = eflag_global = vflag_global =
eflag_atom = vflag_atom = 0;
if (vflag_global) FORTRAN(cbkvirial, CBKVIRIAL).Lvirial = 1;
else FORTRAN(cbkvirial, CBKVIRIAL).Lvirial = 0;
if (vflag_atom) FORTRAN(cbkvirial, CBKVIRIAL).Latomvirial = 1;
else FORTRAN(cbkvirial, CBKVIRIAL).Latomvirial = 0;
// reallocate charge equilibration and CG arrays if necessary
if (atom->nmax > nmax) {
memory->destroy(rcg);
memory->destroy(wcg);
memory->destroy(pcg);
memory->destroy(poldcg);
memory->destroy(qcg);
nmax = atom->nmax;
int n = nmax+1;
memory->create(arow_ptr,n,"reax:arow_ptr");
memory->create(ch,n,"reax:ch");
memory->create(elcvec,n,"reax:elcvec");
memory->create(rcg,n,"reax:rcg");
memory->create(wcg,n,"reax:wcg");
memory->create(pcg,n,"reax:pcg");
memory->create(poldcg,n,"reax:poldcg");
memory->create(qcg,n,"reax:qcg");
}
// calculate the atomic charge distribution
compute_charge(energy_charge_equilibration);
// transfer LAMMPS positions and neighbor lists to REAX
write_reax_positions();
write_reax_vlist();
// determine whether this bond is owned by the processor or not
FORTRAN(srtbon1, SRTBON1)(&iprune, &ihb, &hbcut, &ihbnew, &itripstaball);
// communicate with other processors for the atomic bond order calculations
FORTRAN(cbkabo, CBKABO).abo;
// communicate local atomic bond order to ghost atomic bond order
packflag = 0;
comm->forward_comm_pair(this);
FORTRAN(molec, MOLEC)();
FORTRAN(encalc, ENCALC)();
FORTRAN(mdsav, MDSAV)(&comm->me);
// read forces from ReaxFF Fortran
read_reax_forces();
// extract global and per-atom energy from ReaxFF Fortran
// compute_charge already contributed to eatom
if (eflag_global) {
evdwl += FORTRAN(cbkenergies, CBKENERGIES).eb;
evdwl += FORTRAN(cbkenergies, CBKENERGIES).ea;
evdwl += FORTRAN(cbkenergies, CBKENERGIES).elp;
evdwl += FORTRAN(cbkenergies, CBKENERGIES).emol;
evdwl += FORTRAN(cbkenergies, CBKENERGIES).ev;
evdwl += FORTRAN(cbkenergies, CBKENERGIES).epen;
evdwl += FORTRAN(cbkenergies, CBKENERGIES).ecoa;
evdwl += FORTRAN(cbkenergies, CBKENERGIES).ehb;
evdwl += FORTRAN(cbkenergies, CBKENERGIES).et;
evdwl += FORTRAN(cbkenergies, CBKENERGIES).eco;
evdwl += FORTRAN(cbkenergies, CBKENERGIES).ew;
evdwl += FORTRAN(cbkenergies, CBKENERGIES).efi;
ecoul += FORTRAN(cbkenergies, CBKENERGIES).ep;
ecoul += energy_charge_equilibration;
eng_vdwl += evdwl;
eng_coul += ecoul;
// Store the different parts of the energy
// in a list for output by compute pair command
pvector[0] = FORTRAN(cbkenergies, CBKENERGIES).eb;
pvector[1] = FORTRAN(cbkenergies, CBKENERGIES).ea;
pvector[2] = FORTRAN(cbkenergies, CBKENERGIES).elp;
pvector[3] = FORTRAN(cbkenergies, CBKENERGIES).emol;
pvector[4] = FORTRAN(cbkenergies, CBKENERGIES).ev;
pvector[5] = FORTRAN(cbkenergies, CBKENERGIES).epen;
pvector[6] = FORTRAN(cbkenergies, CBKENERGIES).ecoa;
pvector[7] = FORTRAN(cbkenergies, CBKENERGIES).ehb;
pvector[8] = FORTRAN(cbkenergies, CBKENERGIES).et;
pvector[9] = FORTRAN(cbkenergies, CBKENERGIES).eco;
pvector[10] = FORTRAN(cbkenergies, CBKENERGIES).ew;
pvector[11] = FORTRAN(cbkenergies, CBKENERGIES).ep;
pvector[12] = FORTRAN(cbkenergies, CBKENERGIES).efi;
pvector[13] = energy_charge_equilibration;
}
if (eflag_atom) {
int ntotal = atom->nlocal + atom->nghost;
for (i = 0; i < ntotal; i++)
eatom[i] += FORTRAN(cbkd,CBKD).estrain[i];
}
// extract global and per-atom virial from ReaxFF Fortran
if (vflag_global) {
virial[0] = -FORTRAN(cbkvirial, CBKVIRIAL).virial[0];
virial[1] = -FORTRAN(cbkvirial, CBKVIRIAL).virial[1];
virial[2] = -FORTRAN(cbkvirial, CBKVIRIAL).virial[2];
virial[3] = -FORTRAN(cbkvirial, CBKVIRIAL).virial[3];
virial[4] = -FORTRAN(cbkvirial, CBKVIRIAL).virial[4];
virial[5] = -FORTRAN(cbkvirial, CBKVIRIAL).virial[5];
}
if (vflag_atom) {
int ntotal = atom->nlocal + atom->nghost;
j = 0;
for (i = 0; i < ntotal; i++) {
vatom[i][0] = -FORTRAN(cbkvirial, CBKVIRIAL).atomvirial[j+0];
vatom[i][1] = -FORTRAN(cbkvirial, CBKVIRIAL).atomvirial[j+1];
vatom[i][2] = -FORTRAN(cbkvirial, CBKVIRIAL).atomvirial[j+2];
vatom[i][3] = -FORTRAN(cbkvirial, CBKVIRIAL).atomvirial[j+3];
vatom[i][4] = -FORTRAN(cbkvirial, CBKVIRIAL).atomvirial[j+4];
vatom[i][5] = -FORTRAN(cbkvirial, CBKVIRIAL).atomvirial[j+5];
j += 6;
}
}
}
/* ---------------------------------------------------------------------- */
void PairREAX::write_reax_positions()
{
int j, jx, jy, jz, jia;
double **x = atom->x;
double *q = atom->q;
int *type = atom->type;
int *tag = atom->tag;
int nlocal = atom->nlocal;
int nghost = atom->nghost;
FORTRAN(rsmall, RSMALL).na = nlocal+nghost;
FORTRAN(rsmall, RSMALL).na_local = nlocal;
if (nlocal+nghost > ReaxParams::nat)
error->one(FLERR,"Reax_defs.h setting for NATDEF is too small");
jx = 0;
jy = ReaxParams::nat;
jz = 2*ReaxParams::nat;
jia = 0;
j = 0;
for (int i = 0; i < nlocal+nghost; i++, j++) {
FORTRAN(cbkc, CBKC).c[j+jx] = x[i][0];
FORTRAN(cbkc, CBKC).c[j+jy] = x[i][1];
FORTRAN(cbkc, CBKC).c[j+jz] = x[i][2];
FORTRAN(cbkch, CBKCH).ch[j] = q[i];
FORTRAN(cbkia, CBKIA).ia[j+jia] = map[type[i]];
FORTRAN(cbkia, CBKIA).iag[j+jia] = map[type[i]];
FORTRAN(cbkc, CBKC).itag[j] = tag[i];
}
}
/* ---------------------------------------------------------------------- */
void PairREAX::write_reax_vlist()
{
int ii, jj, i, j, iii, jjj;
double xitmp, yitmp, zitmp;
double xjtmp, yjtmp, zjtmp;
int itag,jtag;
int nvpair, nvlself, nvpairmax;
int nbond;
int inum,jnum;
int *ilist,*jlist,*numneigh,**firstneigh;
double delr2;
double delx, dely, delz;
double **x = atom->x;
int *tag = atom->tag;
int nlocal = atom->nlocal;
int nghost = atom->nghost;
nvpairmax = ReaxParams::nneighmax * ReaxParams::nat;
nvpair = 0;
nvlself =0;
nbond = 0;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
xitmp = x[i][0];
yitmp = x[i][1];
zitmp = x[i][2];
itag = tag[i];
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
xjtmp = x[j][0];
yjtmp = x[j][1];
zjtmp = x[j][2];
jtag = tag[j];
delx = xitmp - xjtmp;
dely = yitmp - yjtmp;
delz = zitmp - zjtmp;
delr2 = delx*delx+dely*dely+delz*delz;
if (delr2 <= rcutvsq) {
if (i < j) {
iii = i+1;
jjj = j+1;
} else {
iii = j+1;
jjj = i+1;
}
if (nvpair >= nvpairmax)
error->one(FLERR,"Reax_defs.h setting for NNEIGHMAXDEF is too small");
FORTRAN(cbkpairs, CBKPAIRS).nvl1[nvpair] = iii;
FORTRAN(cbkpairs, CBKPAIRS).nvl2[nvpair] = jjj;
FORTRAN(cbknvlbo, CBKNVLBO).nvlbo[nvpair] = 0;
if (delr2 <= rcutbsq) {
FORTRAN(cbknvlbo, CBKNVLBO).nvlbo[nvpair] = 1;
nbond++;
}
FORTRAN(cbknvlown, CBKNVLOWN).nvlown[nvpair] = 0;
if (j < nlocal)
FORTRAN(cbknvlown, CBKNVLOWN).nvlown[nvpair] = 1;
else if (itag < jtag)
FORTRAN(cbknvlown, CBKNVLOWN).nvlown[nvpair] = 1;
else if (itag == jtag) {
if (delz > SMALL)
FORTRAN(cbknvlown, CBKNVLOWN).nvlown[nvpair] = 1;
else if (fabs(delz) < SMALL) {
if (dely > SMALL)
FORTRAN(cbknvlown, CBKNVLOWN).nvlown[nvpair] = 1;
else if (fabs(dely) < SMALL && delx > SMALL)
FORTRAN(cbknvlown, CBKNVLOWN).nvlown[nvpair] = 1;
}
}
nvpair++;
}
}
}
int ntotal = nlocal + nghost;
for (int i = nlocal; i < ntotal; i++) {
xitmp = x[i][0];
yitmp = x[i][1];
zitmp = x[i][2];
itag = tag[i];
for (int j = i+1; j < ntotal; j++) {
xjtmp = x[j][0];
yjtmp = x[j][1];
zjtmp = x[j][2];
jtag = tag[j];
delx = xitmp - xjtmp;
dely = yitmp - yjtmp;
delz = zitmp - zjtmp;
delr2 = delx*delx+dely*dely+delz*delz;
// don't need to check the double count since i < j in the ghost region
if (delr2 <= rcutvsq) {
iii = i+1;
jjj = j+1;
if (nvpair >= nvpairmax)
error->one(FLERR,"Reax_defs.h setting for NNEIGHMAXDEF is too small");
FORTRAN(cbkpairs, CBKPAIRS).nvl1[nvpair] = iii;
FORTRAN(cbkpairs, CBKPAIRS).nvl2[nvpair] = jjj;
FORTRAN(cbknvlbo, CBKNVLBO).nvlbo[nvpair] = 0;
if (delr2 <= rcutbsq) {
FORTRAN(cbknvlbo, CBKNVLBO).nvlbo[nvpair] = 1;
nbond++;
}
FORTRAN(cbknvlown, CBKNVLOWN).nvlown[nvpair] = 0;
nvpair++;
}
}
}
FORTRAN(cbkpairs, CBKPAIRS).nvpair = nvpair;
FORTRAN(cbkpairs, CBKPAIRS).nvlself = nvlself;
}
/* ---------------------------------------------------------------------- */
void PairREAX::read_reax_forces()
{
double ftmp[3];
double **f = atom->f;
int ntotal = atom->nlocal + atom->nghost;
int j = 0;
for (int i = 0; i < ntotal; i++) {
ftmp[0] = -FORTRAN(cbkd, CBKD).d[j];
ftmp[1] = -FORTRAN(cbkd, CBKD).d[j+1];
ftmp[2] = -FORTRAN(cbkd, CBKD).d[j+2];
f[i][0] = ftmp[0];
f[i][1] = ftmp[1];
f[i][2] = ftmp[2];
j += 3;
}
}
/* ---------------------------------------------------------------------- */
void PairREAX::allocate()
{
allocated = 1;
int n = atom->ntypes;
memory->create(setflag,n+1,n+1,"pair:setflag");
memory->create(cutsq,n+1,n+1,"pair:cutsq");
param_list = new ff_params[n+1];
for (int i = 1; i <= n; i++)
param_list[i].params = new double[5];
map = new int[n+1];
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairREAX::settings(int narg, char **arg)
{
if (narg != 0 && narg !=4) error->all(FLERR,"Illegal pair_style command");
if (narg == 4) {
hbcut = force->numeric(arg[0]);
ihbnew = static_cast<int> (force->numeric(arg[1]));
itripstaball = static_cast<int> (force->numeric(arg[2]));
precision = force->numeric(arg[3]);
if (hbcut <= 0.0 ||
(ihbnew != 0 && ihbnew != 1) ||
(itripstaball != 0 && itripstaball != 1) ||
precision <= 0.0)
error->all(FLERR,"Illegal pair_style command");
}
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairREAX::coeff(int narg, char **arg)
{
if (!allocated) allocate();
if (narg != 3 + atom->ntypes)
error->all(FLERR,"Incorrect args for pair coefficients");
// insure I,J args are * *
if (strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0)
error->all(FLERR,"Incorrect args for pair coefficients");
// insure filename is ffield.reax
if (strcmp(arg[2],"ffield.reax") != 0)
error->all(FLERR,"Incorrect args for pair coefficients");
// read args that map atom types to elements in potential file
// map[i] = which element the Ith atom type is, -1 if NULL
// NOTE: for now throw an error if NULL is used to disallow use with hybrid
// qEq atrix solver needs to be modified to exclude atoms
for (int i = 3; i < narg; i++) {
if (strcmp(arg[i],"NULL") == 0) {
map[i-2] = -1;
error->all(FLERR,"Cannot currently use pair reax with pair hybrid");
continue;
}
map[i-2] = force->inumeric(arg[i]);
}
int n = atom->ntypes;
int count = 0;
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++) {
setflag[i][j] = 1;
count++;
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void PairREAX::init_style()
{
if (atom->tag_enable == 0)
error->all(FLERR,"Pair style reax requires atom IDs");
if (force->newton_pair == 0)
error->all(FLERR,"Pair style reax requires newton pair on");
if (strcmp(update->unit_style,"real") != 0 && comm->me == 0)
error->warning(FLERR,"Not using real units with pair reax");
int irequest = neighbor->request(this);
neighbor->requests[irequest]->newton = 2;
FORTRAN(readc, READC)();
FORTRAN(reaxinit, REAXINIT)();
FORTRAN(ffinpt, FFINPT)();
FORTRAN(tap7th, TAP7TH)();
// turn off read_in by fort.3 in REAX Fortran
int ngeofor_tmp = -1;
FORTRAN(setngeofor, SETNGEOFOR)(&ngeofor_tmp);
if (comm->me == 0) FORTRAN(readgeo, READGEO)();
// initial setup for cutoff radius of VLIST and BLIST in ReaxFF
double vlbora;
FORTRAN(getswb, GETSWB)(&swb);
cutmax=MAX(swb, hbcut);
rcutvsq=cutmax*cutmax;
FORTRAN(getvlbora, GETVLBORA)(&vlbora);
rcutbsq=vlbora*vlbora;
// parameters for charge equilibration from ReaxFF input, fort.4
// verify that no LAMMPS type to REAX type mapping was invalid
int nelements;
FORTRAN(getnso, GETNSO)(&nelements);
FORTRAN(getswa, GETSWA)(&swa);
double chi, eta, gamma;
for (int itype = 1; itype <= atom->ntypes; itype++) {
if (map[itype] < 1 || map[itype] > nelements)
error->all(FLERR,"Invalid REAX atom type");
chi = FORTRAN(cbkchb, CBKCHB).chi[map[itype]-1];
eta = FORTRAN(cbkchb, CBKCHB).eta[map[itype]-1];
gamma = FORTRAN(cbkchb, CBKCHB).gam[map[itype]-1];
param_list[itype].np = 5;
param_list[itype].rcutsq = cutmax;
param_list[itype].params[0] = chi;
param_list[itype].params[1] = eta;
param_list[itype].params[2] = gamma;
param_list[itype].params[3] = swa;
param_list[itype].params[4] = swb;
}
taper_setup();
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairREAX::init_one(int i, int j)
{
return cutmax;
}
/* ---------------------------------------------------------------------- */
int PairREAX::pack_comm(int n, int *list, double *buf, int pbc_flag, int *pbc)
{
int i,j,m;
m = 0;
if (packflag == 0) {
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = FORTRAN(cbkabo, CBKABO).abo[j];
}
} else {
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = wcg[j];
}
}
return 1;
}
/* ---------------------------------------------------------------------- */
void PairREAX::unpack_comm(int n, int first, double *buf)
{
int i,m,last;
m = 0;
last = first + n;
if (packflag == 0) {
for (i = first; i < last; i++)
FORTRAN(cbkabo, CBKABO).abo[i] = buf[m++];
} else {
for (i = first; i < last; i++)
wcg[i] = buf[m++];
}
}
/* ---------------------------------------------------------------------- */
int PairREAX::pack_reverse_comm(int n, int first, double *buf)
{
int i,m,last;
m = 0;
last = first + n;
for (i = first; i < last; i++)
buf[m++] = wcg[i];
return 1;
}
/* ---------------------------------------------------------------------- */
void PairREAX::unpack_reverse_comm(int n, int *list, double *buf)
{
int i,j,m;
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
wcg[j] += buf[m++];
}
}
/* ----------------------------------------------------------------------
charge equilibration routines
------------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- */
void PairREAX::taper_setup()
{
double swb2,swa2,swb3,swa3,d1,d7;
d1=swb-swa;
d7=pow(d1,7);
swa2=swa*swa;
swa3=swa2*swa;
swb2=swb*swb;
swb3=swb2*swb;
swc7= 20.0e0/d7;
swc6= -70.0e0*(swa+swb)/d7;
swc5= 84.0e0*(swa2+3.0e0*swa*swb+swb2)/d7;
swc4= -35.0e0*(swa3+9.0e0*swa2*swb+9.0e0*swa*swb2+swb3)/d7;
swc3= 140.0e0*(swa3*swb+3.0e0*swa2*swb2+swa*swb3)/d7;
swc2=-210.0e0*(swa3*swb2+swa2*swb3)/d7;
swc1= 140.0e0*swa3*swb3/d7;
swc0=(-35.0e0*swa3*swb2*swb2+21.0e0*swa2*swb3*swb2+
7.0e0*swa*swb3*swb3+swb3*swb3*swb)/d7;
}
/* ---------------------------------------------------------------------- */
double PairREAX::taper_E(const double &r, const double &r2)
{
double r3=r2*r;
return swc7*r3*r3*r+swc6*r3*r3+swc5*r3*r2+swc4*r2*r2+swc3*r3+swc2*r2+
swc1*r+swc0;
}
/* ---------------------------------------------------------------------- */
double PairREAX::taper_F(const double &r, const double &r2)
{
double r3=r2*r;
return 7.0e0*swc7*r3*r3+6.0e0*swc6*r3*r2+5.0e0*swc5*r2*r2+
4.0e0*swc4*r3+3.0e0*swc3*r2+2.0e0*swc2*r+swc1;
}
/* ----------------------------------------------------------------------
compute current charge distributions based on the charge equilibration
------------------------------------------------------------------------- */
void PairREAX::compute_charge(double &energy_charge_equilibration)
{
double xitmp, yitmp, zitmp;
double xjtmp, yjtmp, zjtmp;
int itype, jtype, itag, jtag;
int ii, jj, i, j;
double delr2, delr_norm, gamt, hulp1, hulp2;
double delx, dely, delz;
double qsum,qi;
int nmatentries;
double sw;
int inum,jnum;
int *ilist,*jlist,*numneigh,**firstneigh;
double **x = atom->x;
double *q = atom->q;
int *type = atom->type;
int *tag = atom->tag;
int nlocal = atom->nlocal;
int nghost = atom->nghost;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// realloc neighbor based arrays if necessary
int numneigh_total = 0;
for (ii = 0; ii < inum; ii++)
numneigh_total += numneigh[ilist[ii]];
if (numneigh_total + 2*nlocal > matmax) {
memory->destroy(aval);
memory->destroy(acol_ind);
matmax = numneigh_total + 2*nlocal;
memory->create(aval,matmax,"reax:aval");
memory->create(acol_ind,matmax,"reax:acol_ind");
}
// build linear system
nmatentries = 0;
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
xitmp = x[i][0];
yitmp = x[i][1];
zitmp = x[i][2];
itype = type[i];
itag = tag[i];
jlist = firstneigh[i];
jnum = numneigh[i];
arow_ptr[i] = nmatentries;
aval[nmatentries] = 2.0*param_list[itype].params[1];
acol_ind[nmatentries] = i;
nmatentries++;
aval[nmatentries] = 1.0;
acol_ind[nmatentries] = nlocal + nghost;
nmatentries++;
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
xjtmp = x[j][0];
yjtmp = x[j][1];
zjtmp = x[j][2];
jtype = type[j];
jtag = tag[j];
delx = xitmp - xjtmp;
dely = yitmp - yjtmp;
delz = zitmp - zjtmp;
delr2 = delx*delx+dely*dely+delz*delz;
// avoid counting local-ghost pair twice since
// ReaxFF uses half neigh list with newton off
if (j >= nlocal) {
if (itag > jtag) {
if ((itag+jtag) % 2 == 0) continue;
} else if (itag < jtag) {
if ((itag+jtag) % 2 == 1) continue;
} else {
if (zjtmp < zitmp) continue;
if (zjtmp == zitmp && yjtmp < yitmp) continue;
if (zjtmp == zitmp && yjtmp == yitmp && xjtmp < xitmp) continue;
}
}
// rcutvsq = cutmax*cutmax, in ReaxFF
if (delr2 <= rcutvsq) {
gamt = sqrt(param_list[itype].params[2]*param_list[jtype].params[2]);
delr_norm = sqrt(delr2);
sw = taper_E(delr_norm, delr2);
hulp1=(delr_norm*delr2+(1.0/(gamt*gamt*gamt)));
hulp2=sw*14.40/cbrt(hulp1);
aval[nmatentries] = hulp2;
acol_ind[nmatentries] = j;
nmatentries++;
}
}
}
// in this case, we don't use Midpoint method
// so, we don't need to consider ghost-ghost interactions
// but, need to fill the arow_ptr[] arrays for the ghost atoms
for (i = nlocal; i < nlocal+nghost; i++)
arow_ptr[i] = nmatentries;
arow_ptr[nlocal+nghost] = nmatentries;
// add rhs matentries to linear system
for (ii =0; ii<inum; ii++) {
i = ilist[ii];
itype = type[i];
elcvec[i] = -param_list[itype].params[0];
}
for (i = nlocal; i < nlocal+nghost; i++) elcvec[i] = 0.0;
// assign current charges to charge vector
qsum = 0.0;
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
qi = q[i];
ch[i] = qi;
if (i < nlocal) qsum += qi;
}
for (i = nlocal; i < nlocal+nghost; i++) {
qi = q[i];
ch[i] = qi;
}
double qtot;
MPI_Allreduce(&qsum,&qtot,1,MPI_DOUBLE,MPI_SUM,world);
elcvec[nlocal+nghost] = 0.0;
ch[nlocal+nghost] = chpot;
// solve the linear system using CG sover
charge_reax(nlocal,nghost,ch,aval,acol_ind,arow_ptr,elcvec);
// calculate the charge equilibration energy
energy_charge_equilibration = 0;
// have already updated charge distributions for the current structure
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
itype = type[i];
// 23.02 is the ReaxFF conversion from eV to kcal/mol
// should really use constants.evfactor ~23.06
// but that would break consistency with serial ReaxFF code
// NOTE: this hard-wired real units
// if want other units would have to change params[] in file
qi = 23.02 * (param_list[itype].params[0]*ch[i]+
param_list[itype].params[1]*ch[i]*ch[i]);
energy_charge_equilibration += qi;
if (eflag_atom) eatom[i] += qi;
}
// copy charge vector back to particles from the calculated values
for (i = 0; i < nlocal+nghost; i++) q[i] = ch[i];
chpot = ch[nlocal+nghost];
}
/* ---------------------------------------------------------------------- */
void PairREAX::charge_reax(const int & nlocal, const int & nghost,
double ch[], double aval[], int acol_ind[],
int arow_ptr[], double elcvec[])
{
cg_solve(nlocal,nghost,aval,acol_ind,arow_ptr,ch,elcvec);
}
/* ----------------------------------------------------------------------
CG solver for linear systems
------------------------------------------------------------------------- */
void PairREAX::cg_solve(const int & nlocal, const int & nghost,
double aval[], int acol_ind[], int arow_ptr[],
double x[], double b[])
{
double one, zero, rho, rho_old, alpha, beta, gamma;
int iter, maxiter;
int n;
double sumtmp;
// parallel CG method by A. P. Thompson
// distributed (partial) vectors: b, r, q, A
// accumulated (full) vectors: x, w, p
// r = b-A.x
// w = r (ReverseComm + Comm)
double *r = rcg;
double *w = wcg;
double *p = pcg;
double *p_old = poldcg;
double *q = qcg;
n = nlocal+nghost+1;
one = 1.0;
zero = 0.0;
maxiter = 100;
for (int i = 0; i < n; i++) w[i] = 0;
// construct r = b-Ax
sparse_product(n, nlocal, nghost, aval, acol_ind, arow_ptr, x, r);
// not using BLAS library
for (int i=0; i<n; i++) {
r[i] = b[i] - r[i];
w[i] = r[i];
}
packflag = 1;
comm->reverse_comm_pair(this);
comm->forward_comm_pair(this);
MPI_Allreduce(&w[n-1], &sumtmp, 1, MPI_DOUBLE, MPI_SUM, world);
w[n-1] = sumtmp;
rho_old = one;
for (iter = 1; iter < maxiter; iter++) {
rho = 0.0;
for (int i=0; i<nlocal; i++) rho += w[i]*w[i];
MPI_Allreduce(&rho, &sumtmp, 1, MPI_DOUBLE, MPI_SUM, world);
rho = sumtmp + w[n-1]*w[n-1];
if (rho < precision) break;
for (int i = 0; i<n; i++) p[i] = w[i];
if (iter > 1) {
beta = rho/rho_old;
for (int i = 0; i<n; i++) p[i] += beta*p_old[i];
}
sparse_product(n, nlocal, nghost, aval, acol_ind, arow_ptr, p, q);
gamma = 0.0;
for (int i=0; i<n; i++) gamma += p[i]*q[i];
MPI_Allreduce(&gamma, &sumtmp, 1, MPI_DOUBLE, MPI_SUM, world);
gamma = sumtmp;
alpha = rho/gamma;
for (int i=0; i<n; i++) {
x[i] += alpha*p[i];
r[i] -= alpha*q[i];
w[i] = r[i];
}
comm->reverse_comm_pair(this);
comm->forward_comm_pair(this);
MPI_Allreduce(&w[n-1], &sumtmp, 1, MPI_DOUBLE, MPI_SUM, world);
w[n-1] = sumtmp;
for (int i=0; i<n; i++) p_old[i] = p[i];
rho_old = rho;
}
}
/* ----------------------------------------------------------------------
sparse maxtrix operations
------------------------------------------------------------------------- */
void PairREAX::sparse_product(const int &n, const int &nlocal,
const int &nghost,
double aval[], int acol_ind[], int arow_ptr[],
double *x, double *r)
{
int i,j,jj;
for (i=0; i<n; i++) r[i] = 0.0;
for (i=0; i<nlocal; i++) {
r[i] += aval[arow_ptr[i]]*x[i];
for (j=arow_ptr[i]+1; j<arow_ptr[i+1]; j++) {
jj = acol_ind[j];
r[i] += aval[j]*x[jj];
r[jj] += aval[j]*x[i];
}
}
for (i=nlocal; i<nlocal+nghost; i++)
for (j=arow_ptr[i]; j<arow_ptr[i+1]; j++) {
jj = acol_ind[j];
r[i] += aval[j]*x[jj];
r[jj] += aval[j]*x[i];
}
}
/* ----------------------------------------------------------------------
memory usage of local atom-based arrays
------------------------------------------------------------------------- */
double PairREAX::memory_usage()
{
double bytes = nmax * sizeof(int);
bytes += 7 * nmax * sizeof(double);
bytes += matmax * sizeof(int);
bytes += matmax * sizeof(double);
return bytes;
}
Event Timeline
Log In to Comment