Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F102272511
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
Tue, Feb 18, 23:46
Size
28 KB
Mime Type
text/x-c
Expires
Thu, Feb 20, 23:46 (2 d)
Engine
blob
Format
Raw Data
Handle
24321390
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 MIN(a,b) ((a) < (b) ? (a) : (b))
#define MAX(a,b) ((a) > (b) ? (a) : (b))
#define SMALL 0.0001
/* ---------------------------------------------------------------------- */
PairREAX::PairREAX(LAMMPS *lmp) : Pair(lmp)
{
single_enable = 0;
one_coeff = 1;
no_virial_compute = 1;
cutmax = 0.0;
hbcut = 6.0;
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()
{
if (allocated) {
memory->destroy_2d_int_array(setflag);
memory->destroy_2d_double_array(cutsq);
for (int i = 1; i <= atom->ntypes; i++)
delete [] param_list[i].params;
delete [] param_list;
delete [] map;
}
memory->sfree(arow_ptr);
memory->sfree(ch);
memory->sfree(elcvec);
memory->sfree(rcg);
memory->sfree(wcg);
memory->sfree(pcg);
memory->sfree(poldcg);
memory->sfree(qcg);
memory->sfree(aval);
memory->sfree(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 = 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->sfree(rcg);
memory->sfree(wcg);
memory->sfree(pcg);
memory->sfree(poldcg);
memory->sfree(qcg);
nmax = atom->nmax;
int n = nmax+1;
arow_ptr = (int *) memory->smalloc(n*sizeof(int),"reax:arow_ptr");
ch = (double *) memory->smalloc(n*sizeof(double),"reax:ch");
elcvec = (double *) memory->smalloc(n*sizeof(double),"reax:elcvec");
rcg = (double *) memory->smalloc(n*sizeof(double),"reax:rcg");
wcg = (double *) memory->smalloc(n*sizeof(double),"reax:wcg");
pcg = (double *) memory->smalloc(n*sizeof(double),"reax:pcg");
poldcg = (double *) memory->smalloc(n*sizeof(double),"reax:poldcg");
qcg = (double *) memory->smalloc(n*sizeof(double),"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);
// 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;
}
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()
{
double xtmp, ytmp, ztmp;
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("reax_defs.h::NATDEF 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;
int *jlist;
int *numneigh,**firstneigh;
double delr2;
double delx, dely, delz;
double rtmp[3];
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];
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("reax_defs.h::NNEIGHMAXDEF 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("reax_defs.h::NNEIGHMAXDEF 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;
setflag = memory->create_2d_int_array(n+1,n+1,"pair:setflag");
cutsq = memory->create_2d_double_array(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 !=2) error->all("Illegal pair_style command");
if (narg == 2) {
hbcut = force->numeric(arg[0]);
precision = force->numeric(arg[1]);
if (hbcut <= 0.0 || precision <= 0.0)
error->all("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("Incorrect args for pair coefficients");
// insure I,J args are * *
if (strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0)
error->all("Incorrect args for pair coefficients");
// insure filename is ffield.reax
if (strcmp(arg[2],"ffield.reax") != 0)
error->all("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
for (int i = 3; i < narg; i++) {
if (strcmp(arg[i],"NULL") == 0) {
map[i-2] = -1;
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("Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void PairREAX::init_style()
{
if (atom->tag_enable == 0)
error->all("Pair style reax requires atom IDs");
if (force->newton_pair == 0)
error->all("Pair style reax requires newton pair on");
if (strcmp(update->unit_style,"real") != 0 && comm->me == 0)
error->warning("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("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,k,m,last,size;
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,k,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;
double rtmp[3];
int inum,jnum;
int *ilist;
int *jlist;
int *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->sfree(aval);
memory->sfree(acol_ind);
matmax = numneigh_total + 2*nlocal;
aval = (double *) memory->smalloc(matmax*sizeof(double),"reax:aval");
acol_ind = (int *) memory->smalloc(matmax*sizeof(int),"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];
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[])
{
double chpottmp, suma;
double sumtmp;
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;
}
/* ----------------------------------------------------------------------
print out the itemized energies to the log file
this is not currently called, but should
be made available to the user in some way
------------------------------------------------------------------------- */
void PairREAX::output_itemized_energy(double energy_charge_equilibration)
{
double etmp[14],etmp2[14];
etmp[0] = FORTRAN(cbkenergies, CBKENERGIES).eb;
etmp[1] = FORTRAN(cbkenergies, CBKENERGIES).ea;
etmp[2] = FORTRAN(cbkenergies, CBKENERGIES).elp;
etmp[3] = FORTRAN(cbkenergies, CBKENERGIES).emol;
etmp[4] = FORTRAN(cbkenergies, CBKENERGIES).ev;
etmp[5] = FORTRAN(cbkenergies, CBKENERGIES).epen;
etmp[6] = FORTRAN(cbkenergies, CBKENERGIES).ecoa;
etmp[7] = FORTRAN(cbkenergies, CBKENERGIES).ehb;
etmp[8] = FORTRAN(cbkenergies, CBKENERGIES).et;
etmp[9] = FORTRAN(cbkenergies, CBKENERGIES).eco;
etmp[10] = FORTRAN(cbkenergies, CBKENERGIES).ew;
etmp[11] = FORTRAN(cbkenergies, CBKENERGIES).ep;
etmp[12] = FORTRAN(cbkenergies, CBKENERGIES).efi;
etmp[13] = energy_charge_equilibration;
MPI_Allreduce(etmp,etmp2,14,MPI_DOUBLE,MPI_SUM,world);
if (comm->me == 0) {
fprintf(screen,"eb = %g \n",etmp2[0]);
fprintf(screen,"ea = %g \n",etmp2[1]);
fprintf(screen,"elp = %g \n",etmp2[2]);
fprintf(screen,"emol = %g \n",etmp2[3]);
fprintf(screen,"ecoa = %g \n",etmp2[4]);
fprintf(screen,"epen = %g \n",etmp2[5]);
fprintf(screen,"ecoa = %g \n",etmp2[6]);
fprintf(screen,"ehb = %g \n",etmp2[7]);
fprintf(screen,"et = %g \n",etmp2[8]);
fprintf(screen,"eco = %g \n",etmp2[9]);
fprintf(screen,"ew = %g \n",etmp2[10]);
fprintf(screen,"ep = %g \n",etmp2[11]);
fprintf(screen,"efi = %g \n",etmp2[12]);
fprintf(screen,"eqeq = %g \n",etmp2[13]);
fprintf(logfile,"eb = %g \n",etmp2[0]);
fprintf(logfile,"ea = %g \n",etmp2[1]);
fprintf(logfile,"elp = %g \n",etmp2[2]);
fprintf(logfile,"emol = %g \n",etmp2[3]);
fprintf(logfile,"ecoa = %g \n",etmp2[4]);
fprintf(logfile,"epen = %g \n",etmp2[5]);
fprintf(logfile,"ecoa = %g \n",etmp2[6]);
fprintf(logfile,"ehb = %g \n",etmp2[7]);
fprintf(logfile,"et = %g \n",etmp2[8]);
fprintf(logfile,"eco = %g \n",etmp2[9]);
fprintf(logfile,"ew = %g \n",etmp2[10]);
fprintf(logfile,"ep = %g \n",etmp2[11]);
fprintf(logfile,"efi = %g \n",etmp2[12]);
fprintf(logfile,"eqeq = %g \n",etmp2[13]);
}
}
Event Timeline
Log In to Comment