Homec4science
Diffusion lammps (lammm-master) Browse

lammps/src/REAX/pair_reax.cpp
62b27e591539lammm-master

Last Change
Raw File

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;
}
c4science · Help