Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F90325715
pair_eim.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
Thu, Oct 31, 13:09
Size
33 KB
Mime Type
text/x-c
Expires
Sat, Nov 2, 13:09 (2 d)
Engine
blob
Format
Raw Data
Handle
21991304
Attached To
rLAMMPS lammps
pair_eim.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 author: Xiaowang Zhou (SNL)
------------------------------------------------------------------------- */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "pair_eim.h"
#include "atom.h"
#include "force.h"
#include "comm.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
#define MAXLINE 1024
/* ---------------------------------------------------------------------- */
PairEIM::PairEIM(LAMMPS *lmp) : Pair(lmp)
{
single_enable = 0;
restartinfo = 0;
one_coeff = 1;
manybody_flag = 1;
setfl = NULL;
nmax = 0;
rho = NULL;
fp = NULL;
map = NULL;
nelements = 0;
elements = NULL;
negativity = NULL;
q0 = NULL;
cutforcesq = NULL;
Fij = NULL;
Gij = NULL;
phiij = NULL;
Fij_spline = NULL;
Gij_spline = NULL;
phiij_spline = NULL;
// set comm size needed by this Pair
comm_forward = 1;
comm_reverse = 1;
}
/* ----------------------------------------------------------------------
check if allocated, since class can be destructed when incomplete
------------------------------------------------------------------------- */
PairEIM::~PairEIM()
{
memory->destroy(rho);
memory->destroy(fp);
if (allocated) {
memory->destroy(setflag);
memory->destroy(cutsq);
delete [] map;
memory->destroy(type2Fij);
memory->destroy(type2Gij);
memory->destroy(type2phiij);
}
for (int i = 0; i < nelements; i++) delete [] elements[i];
delete [] elements;
deallocate_setfl();
delete [] negativity;
delete [] q0;
memory->destroy(cutforcesq);
memory->destroy(Fij);
memory->destroy(Gij);
memory->destroy(phiij);
memory->destroy(Fij_spline);
memory->destroy(Gij_spline);
memory->destroy(phiij_spline);
}
/* ---------------------------------------------------------------------- */
void PairEIM::compute(int eflag, int vflag)
{
int i,j,ii,jj,m,inum,jnum,itype,jtype;
double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
double rsq,r,p,rhoip,rhojp,phip,phi,coul,coulp,recip,psip;
double *coeff;
int *ilist,*jlist,*numneigh,**firstneigh;
evdwl = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = eflag_global = eflag_atom = 0;
// grow energy array if necessary
if (atom->nmax > nmax) {
memory->destroy(rho);
memory->destroy(fp);
nmax = atom->nmax;
memory->create(rho,nmax,"pair:rho");
memory->create(fp,nmax,"pair:fp");
}
double **x = atom->x;
double **f = atom->f;
int *type = atom->type;
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// zero out density
if (newton_pair) {
m = nlocal + atom->nghost;
for (i = 0; i < m; i++) {
rho[i] = 0.0;
fp[i] = 0.0;
}
} else {
for (i = 0; i < nlocal; i++) {
rho[i] = 0.0;
fp[i] = 0.0;
}
}
// rho = density at each atom
// loop over neighbors of my atoms
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
jtype = type[j];
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
if (rsq < cutforcesq[itype][jtype]) {
p = sqrt(rsq)*rdr + 1.0;
m = static_cast<int> (p);
m = MIN(m,nr-1);
p -= m;
p = MIN(p,1.0);
coeff = Fij_spline[type2Fij[itype][jtype]][m];
rho[i] += ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6];
if (newton_pair || j < nlocal) {
coeff = Fij_spline[type2Fij[jtype][itype]][m];
rho[j] += ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6];
}
}
}
}
// communicate and sum densities
rhofp = 1;
if (newton_pair) comm->reverse_comm_pair(this);
comm->forward_comm_pair(this);
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
jtype = type[j];
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
if (rsq < cutforcesq[itype][jtype]) {
p = sqrt(rsq)*rdr + 1.0;
m = static_cast<int> (p);
m = MIN(m,nr-1);
p -= m;
p = MIN(p,1.0);
coeff = Gij_spline[type2Gij[itype][jtype]][m];
fp[i] += rho[j]*(((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6]);
if (newton_pair || j < nlocal) {
fp[j] += rho[i]*(((coeff[3]*p + coeff[4])*p + coeff[5])*p +
coeff[6]);
}
}
}
}
// communicate and sum modified densities
rhofp = 2;
if (newton_pair) comm->reverse_comm_pair(this);
comm->forward_comm_pair(this);
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
itype = type[i];
if (eflag) {
phi = 0.5*rho[i]*fp[i];
if (eflag_global) eng_vdwl += phi;
if (eflag_atom) eatom[i] += phi;
}
}
// compute forces on each atom
// loop over neighbors of my atoms
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
jtype = type[j];
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
if (rsq < cutforcesq[itype][jtype]) {
r = sqrt(rsq);
p = r*rdr + 1.0;
m = static_cast<int> (p);
m = MIN(m,nr-1);
p -= m;
p = MIN(p,1.0);
// rhoip = derivative of (density at atom j due to atom i)
// rhojp = derivative of (density at atom i due to atom j)
// phi = pair potential energy
// phip = phi'
coeff = Fij_spline[type2Fij[jtype][itype]][m];
rhoip = (coeff[0]*p + coeff[1])*p + coeff[2];
coeff = Fij_spline[type2Fij[itype][jtype]][m];
rhojp = (coeff[0]*p + coeff[1])*p + coeff[2];
coeff = phiij_spline[type2phiij[itype][jtype]][m];
phip = (coeff[0]*p + coeff[1])*p + coeff[2];
phi = ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6];
coeff = Gij_spline[type2Gij[itype][jtype]][m];
coul = ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6];
coulp = (coeff[0]*p + coeff[1])*p + coeff[2];
psip = phip + (rho[i]*rho[j]-q0[itype]*q0[jtype])*coulp +
fp[i]*rhojp + fp[j]*rhoip;
recip = 1.0/r;
fpair = -psip*recip;
f[i][0] += delx*fpair;
f[i][1] += dely*fpair;
f[i][2] += delz*fpair;
if (newton_pair || j < nlocal) {
f[j][0] -= delx*fpair;
f[j][1] -= dely*fpair;
f[j][2] -= delz*fpair;
}
if (eflag) evdwl = phi-q0[itype]*q0[jtype]*coul;
if (evflag) ev_tally(i,j,nlocal,newton_pair,
evdwl,0.0,fpair,delx,dely,delz);
}
}
}
if (vflag_fdotr) virial_fdotr_compute();
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
void PairEIM::allocate()
{
allocated = 1;
int n = atom->ntypes;
memory->create(setflag,n+1,n+1,"pair:setflag");
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
setflag[i][j] = 0;
memory->create(cutsq,n+1,n+1,"pair:cutsq");
map = new int[n+1];
for (int i = 1; i <= n; i++) map[i] = -1;
memory->create(type2Fij,n+1,n+1,"pair:type2Fij");
memory->create(type2Gij,n+1,n+1,"pair:type2Gij");
memory->create(type2phiij,n+1,n+1,"pair:type2phiij");
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairEIM::settings(int narg, char **arg)
{
if (narg > 0) error->all(FLERR,"Illegal pair_style command");
}
/* ----------------------------------------------------------------------
set coeffs from set file
------------------------------------------------------------------------- */
void PairEIM::coeff(int narg, char **arg)
{
int i,j,m,n;
if (!allocated) allocate();
if (narg < 5) 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");
// read EIM element names before filename
// nelements = # of EIM elements to read from file
// elements = list of unique element names
if (nelements) {
for (i = 0; i < nelements; i++) delete [] elements[i];
delete [] elements;
}
nelements = narg - 3 - atom->ntypes;
if (nelements < 1) error->all(FLERR,"Incorrect args for pair coefficients");
elements = new char*[nelements];
for (i = 0; i < nelements; i++) {
n = strlen(arg[i+2]) + 1;
elements[i] = new char[n];
strcpy(elements[i],arg[i+2]);
}
// read EIM file
deallocate_setfl();
setfl = new Setfl();
read_file(arg[2+nelements]);
// read args that map atom types to elements in potential file
// map[i] = which element the Ith atom type is, -1 if NULL
for (i = 3 + nelements; i < narg; i++) {
m = i - (3+nelements) + 1;
for (j = 0; j < nelements; j++)
if (strcmp(arg[i],elements[j]) == 0) break;
if (j < nelements) map[m] = j;
else if (strcmp(arg[i],"NULL") == 0) map[m] = -1;
else error->all(FLERR,"Incorrect args for pair coefficients");
}
// clear setflag since coeff() called once with I,J = * *
n = atom->ntypes;
for (i = 1; i <= n; i++)
for (j = i; j <= n; j++)
setflag[i][j] = 0;
// set setflag i,j for type pairs where both are mapped to elements
// set mass of atom type if i = j
int count = 0;
for (i = 1; i <= n; i++)
for (j = i; j <= n; j++)
if (map[i] >= 0 && map[j] >= 0) {
setflag[i][j] = 1;
if (i == j) atom->set_mass(i,setfl->mass[map[i]]);
count++;
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void PairEIM::init_style()
{
// convert read-in file(s) to arrays and spline them
file2array();
array2spline();
neighbor->request(this,instance_me);
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairEIM::init_one(int i, int j)
{
cutmax = sqrt(cutforcesq[i][j]);
return cutmax;
}
/* ----------------------------------------------------------------------
read potential values from a set file
------------------------------------------------------------------------- */
void PairEIM::read_file(char *filename)
{
// open potential file
int me = comm->me;
FILE *fptr;
if (me == 0) {
fptr = force->open_potential(filename);
if (fptr == NULL) {
char str[128];
sprintf(str,"Cannot open EIM potential file %s",filename);
error->one(FLERR,str);
}
}
int npair = nelements*(nelements+1)/2;
setfl->ielement = new int[nelements];
setfl->mass = new double[nelements];
setfl->negativity = new double[nelements];
setfl->ra = new double[nelements];
setfl->ri = new double[nelements];
setfl->Ec = new double[nelements];
setfl->q0 = new double[nelements];
setfl->rcutphiA = new double[npair];
setfl->rcutphiR = new double[npair];
setfl->Eb = new double[npair];
setfl->r0 = new double[npair];
setfl->alpha = new double[npair];
setfl->beta = new double[npair];
setfl->rcutq = new double[npair];
setfl->Asigma = new double[npair];
setfl->rq = new double[npair];
setfl->rcutsigma = new double[npair];
setfl->Ac = new double[npair];
setfl->zeta = new double[npair];
setfl->rs = new double[npair];
setfl->tp = new int[npair];
if (me == 0)
if (!grabglobal(fptr))
error->one(FLERR,"Could not grab global entry from EIM potential file");
MPI_Bcast(&setfl->division,1,MPI_DOUBLE,0,world);
MPI_Bcast(&setfl->rbig,1,MPI_DOUBLE,0,world);
MPI_Bcast(&setfl->rsmall,1,MPI_DOUBLE,0,world);
for (int i = 0; i < nelements; i++) {
if (me == 0)
if (!grabsingle(fptr,i))
error->one(FLERR,"Could not grab element entry from EIM potential file");
MPI_Bcast(&setfl->ielement[i],1,MPI_INT,0,world);
MPI_Bcast(&setfl->mass[i],1,MPI_DOUBLE,0,world);
MPI_Bcast(&setfl->negativity[i],1,MPI_DOUBLE,0,world);
MPI_Bcast(&setfl->ra[i],1,MPI_DOUBLE,0,world);
MPI_Bcast(&setfl->ri[i],1,MPI_DOUBLE,0,world);
MPI_Bcast(&setfl->Ec[i],1,MPI_DOUBLE,0,world);
MPI_Bcast(&setfl->q0[i],1,MPI_DOUBLE,0,world);
}
for (int i = 0; i < nelements; i++) {
for (int j = i; j < nelements; j++) {
int ij;
if (i == j) ij = i;
else if (i < j) ij = nelements*(i+1) - (i+1)*(i+2)/2 + j;
else ij = nelements*(j+1) - (j+1)*(j+2)/2 + i;
if (me == 0)
if (grabpair(fptr,i,j) == 0)
error->one(FLERR,"Could not grab pair entry from EIM potential file");
MPI_Bcast(&setfl->rcutphiA[ij],1,MPI_DOUBLE,0,world);
MPI_Bcast(&setfl->rcutphiR[ij],1,MPI_DOUBLE,0,world);
MPI_Bcast(&setfl->Eb[ij],1,MPI_DOUBLE,0,world);
MPI_Bcast(&setfl->r0[ij],1,MPI_DOUBLE,0,world);
MPI_Bcast(&setfl->alpha[ij],1,MPI_DOUBLE,0,world);
MPI_Bcast(&setfl->beta[ij],1,MPI_DOUBLE,0,world);
MPI_Bcast(&setfl->rcutq[ij],1,MPI_DOUBLE,0,world);
MPI_Bcast(&setfl->Asigma[ij],1,MPI_DOUBLE,0,world);
MPI_Bcast(&setfl->rq[ij],1,MPI_DOUBLE,0,world);
MPI_Bcast(&setfl->rcutsigma[ij],1,MPI_DOUBLE,0,world);
MPI_Bcast(&setfl->Ac[ij],1,MPI_DOUBLE,0,world);
MPI_Bcast(&setfl->zeta[ij],1,MPI_DOUBLE,0,world);
MPI_Bcast(&setfl->rs[ij],1,MPI_DOUBLE,0,world);
MPI_Bcast(&setfl->tp[ij],1,MPI_INT,0,world);
}
}
setfl->nr = 5000;
setfl->cut = 0.0;
for (int i = 0; i < npair; i++) {
if (setfl->cut < setfl->rcutphiA[i]) setfl->cut = setfl->rcutphiA[i];
if (setfl->cut < setfl->rcutphiR[i]) setfl->cut = setfl->rcutphiR[i];
if (setfl->cut < setfl->rcutq[i]) setfl->cut = setfl->rcutq[i];
if (setfl->cut < setfl->rcutsigma[i]) setfl->cut = setfl->rcutsigma[i];
}
setfl->dr = setfl->cut/(setfl->nr-1.0);
memory->create(setfl->cuts,nelements,nelements,"pair:cuts");
for (int i = 0; i < nelements; i++) {
for (int j = 0; j < nelements; j++) {
if (i > j) {
setfl->cuts[i][j] = setfl->cuts[j][i];
} else {
int ij;
if (i == j) {
ij = i;
} else {
ij = nelements*(i+1) - (i+1)*(i+2)/2 + j;
}
setfl->cuts[i][j] = setfl->rcutphiA[ij];
if (setfl->cuts[i][j] < setfl->rcutphiR[ij])
setfl->cuts[i][j] = setfl->rcutphiR[ij];
if (setfl->cuts[i][j] < setfl->rcutq[ij])
setfl->cuts[i][j] = setfl->rcutq[ij];
if (setfl->cuts[i][j] < setfl->rcutsigma[ij])
setfl->cuts[i][j] = setfl->rcutsigma[ij];
}
}
}
memory->create(setfl->Fij,nelements,nelements,setfl->nr+1,"pair:Fij");
memory->create(setfl->Gij,nelements,nelements,setfl->nr+1,"pair:Gij");
memory->create(setfl->phiij,nelements,nelements,setfl->nr+1,"pair:phiij");
for (int i = 0; i < nelements; i++)
for (int j = 0; j < nelements; j++) {
for (int k = 0; k < setfl->nr; k++) {
if (i > j) {
setfl->phiij[i][j][k+1] = setfl->phiij[j][i][k+1];
} else {
double r = k*setfl->dr;
setfl->phiij[i][j][k+1] = funcphi(i,j,r);
}
}
}
for (int i = 0; i < nelements; i++)
for (int j = 0; j < nelements; j++) {
for (int k = 0; k < setfl->nr; k++) {
double r = k*setfl->dr;
setfl->Fij[i][j][k+1] = funcsigma(i,j,r);
}
}
for (int i = 0; i < nelements; i++)
for (int j = 0; j < nelements; j++) {
for (int k = 0; k < setfl->nr; k++) {
if (i > j) {
setfl->Gij[i][j][k+1] = setfl->Gij[j][i][k+1];
} else {
double r = k*setfl->dr;
setfl->Gij[i][j][k+1] = funccoul(i,j,r);
}
}
}
// close the potential file
if (me == 0) fclose(fptr);
}
/* ----------------------------------------------------------------------
deallocate data associated with setfl file
------------------------------------------------------------------------- */
void PairEIM::deallocate_setfl()
{
if (!setfl) return;
delete [] setfl->ielement;
delete [] setfl->mass;
delete [] setfl->negativity;
delete [] setfl->ra;
delete [] setfl->ri;
delete [] setfl->Ec;
delete [] setfl->q0;
delete [] setfl->rcutphiA;
delete [] setfl->rcutphiR;
delete [] setfl->Eb;
delete [] setfl->r0;
delete [] setfl->alpha;
delete [] setfl->beta;
delete [] setfl->rcutq;
delete [] setfl->Asigma;
delete [] setfl->rq;
delete [] setfl->rcutsigma;
delete [] setfl->Ac;
delete [] setfl->zeta;
delete [] setfl->rs;
delete [] setfl->tp;
memory->destroy(setfl->cuts);
memory->destroy(setfl->Fij);
memory->destroy(setfl->Gij);
memory->destroy(setfl->phiij);
delete setfl;
}
/* ----------------------------------------------------------------------
convert read-in potentials to standard array format
interpolate all file values to a single grid and cutoff
------------------------------------------------------------------------- */
void PairEIM::file2array()
{
int i,j,m,n;
int irow,icol;
int ntypes = atom->ntypes;
delete [] negativity;
delete [] q0;
delete [] cutforcesq;
negativity = new double[ntypes+1];
q0 = new double[ntypes+1];
memory->create(cutforcesq,ntypes+1,ntypes+1,"pair:cutforcesq");
for (i = 1; i <= ntypes; i++) {
if (map[i] == -1) {
negativity[i]=0.0;
q0[i]=0.0;
} else {
negativity[i]=setfl->negativity[map[i]];
q0[i]=setfl->q0[map[i]];
}
}
for (i = 1; i <= ntypes; i++)
for (j = 1; j <= ntypes; j++) {
if (map[i] == -1 || map[j] == -1) {
cutforcesq[i][j] = setfl->cut;
cutforcesq[i][j] = cutforcesq[i][j]*cutforcesq[i][j];
} else {
cutforcesq[i][j] = setfl->cuts[map[i]][map[j]];
cutforcesq[i][j] = cutforcesq[i][j]*cutforcesq[i][j];
}
}
nr = setfl->nr;
dr = setfl->dr;
// ------------------------------------------------------------------
// setup Fij arrays
// ------------------------------------------------------------------
nFij = nelements*nelements + 1;
memory->destroy(Fij);
memory->create(Fij,nFij,nr+1,"pair:Fij");
// copy each element's Fij to global Fij
n=0;
for (i = 0; i < nelements; i++)
for (j = 0; j < nelements; j++) {
for (m = 1; m <= nr; m++) Fij[n][m] = setfl->Fij[i][j][m];
n++;
}
// add extra Fij of zeroes for non-EIM types to point to (pair hybrid)
for (m = 1; m <= nr; m++) Fij[nFij-1][m] = 0.0;
// type2Fij[i][j] = which Fij array (0 to nFij-1) each type pair maps to
// setfl of Fij arrays
// value = n = sum over rows of matrix until reach irow,icol
// if atom type doesn't point to element (non-EIM atom in pair hybrid)
// then map it to last Fij array of zeroes
for (i = 1; i <= ntypes; i++) {
for (j = 1; j <= ntypes; j++) {
irow = map[i];
icol = map[j];
if (irow == -1 || icol == -1) {
type2Fij[i][j] = nFij-1;
} else {
n = 0;
for (m = 0; m < irow; m++) n += nelements;
n += icol;
type2Fij[i][j] = n;
}
}
}
// ------------------------------------------------------------------
// setup Gij arrays
// ------------------------------------------------------------------
nGij = nelements * (nelements+1) / 2 + 1;
memory->destroy(Gij);
memory->create(Gij,nGij,nr+1,"pair:Gij");
// copy each element's Gij to global Gij, only for I >= J
n=0;
for (i = 0; i < nelements; i++)
for (j = 0; j <= i; j++) {
for (m = 1; m <= nr; m++) Gij[n][m] = setfl->Gij[i][j][m];
n++;
}
// add extra Gij of zeroes for non-EIM types to point to (pair hybrid)
for (m = 1; m <= nr; m++) Gij[nGij-1][m] = 0.0;
// type2Gij[i][j] = which Gij array (0 to nGij-1) each type pair maps to
// setfl of Gij arrays only fill lower triangular Nelement matrix
// value = n = sum over rows of lower-triangular matrix until reach irow,icol
// swap indices when irow < icol to stay lower triangular
// if atom type doesn't point to element (non-EIM atom in pair hybrid)
// then map it to last Gij array of zeroes
for (i = 1; i <= ntypes; i++) {
for (j = 1; j <= ntypes; j++) {
irow = map[i];
icol = map[j];
if (irow == -1 || icol == -1) {
type2Gij[i][j] = nGij-1;
} else {
if (irow < icol) {
irow = map[j];
icol = map[i];
}
n = 0;
for (m = 0; m < irow; m++) n += m + 1;
n += icol;
type2Gij[i][j] = n;
}
}
}
// ------------------------------------------------------------------
// setup phiij arrays
// ------------------------------------------------------------------
nphiij = nelements * (nelements+1) / 2 + 1;
memory->destroy(phiij);
memory->create(phiij,nphiij,nr+1,"pair:phiij");
// copy each element pair phiij to global phiij, only for I >= J
n = 0;
for (i = 0; i < nelements; i++)
for (j = 0; j <= i; j++) {
for (m = 1; m <= nr; m++) phiij[n][m] = setfl->phiij[i][j][m];
n++;
}
// add extra phiij of zeroes for non-EIM types to point to (pair hybrid)
for (m = 1; m <= nr; m++) phiij[nphiij-1][m] = 0.0;
// type2phiij[i][j] = which phiij array (0 to nphiij-1)
// each type pair maps to
// setfl of phiij arrays only fill lower triangular Nelement matrix
// value = n = sum over rows of lower-triangular matrix until reach irow,icol
// swap indices when irow < icol to stay lower triangular
// if atom type doesn't point to element (non-EIM atom in pair hybrid)
// then map it to last phiij array of zeroes
for (i = 1; i <= ntypes; i++) {
for (j = 1; j <= ntypes; j++) {
irow = map[i];
icol = map[j];
if (irow == -1 || icol == -1) {
type2phiij[i][j] = nphiij-1;
} else {
if (irow < icol) {
irow = map[j];
icol = map[i];
}
n = 0;
for (m = 0; m < irow; m++) n += m + 1;
n += icol;
type2phiij[i][j] = n;
}
}
}
}
/* ---------------------------------------------------------------------- */
void PairEIM::array2spline()
{
rdr = 1.0/dr;
memory->destroy(Fij_spline);
memory->destroy(Gij_spline);
memory->destroy(phiij_spline);
memory->create(Fij_spline,nFij,nr+1,7,"pair:Fij");
memory->create(Gij_spline,nGij,nr+1,7,"pair:Gij");
memory->create(phiij_spline,nphiij,nr+1,7,"pair:phiij");
for (int i = 0; i < nFij; i++)
interpolate(nr,dr,Fij[i],Fij_spline[i],0.0);
for (int i = 0; i < nGij; i++)
interpolate(nr,dr,Gij[i],Gij_spline[i],0.0);
for (int i = 0; i < nphiij; i++)
interpolate(nr,dr,phiij[i],phiij_spline[i],0.0);
}
/* ---------------------------------------------------------------------- */
void PairEIM::interpolate(int n, double delta, double *f,
double **spline, double origin)
{
for (int m = 1; m <= n; m++) spline[m][6] = f[m];
spline[1][5] = spline[2][6] - spline[1][6];
spline[2][5] = 0.5 * (spline[3][6]-spline[1][6]);
spline[n-1][5] = 0.5 * (spline[n][6]-spline[n-2][6]);
spline[n][5] = 0.0;
for (int m = 3; m <= n-2; m++)
spline[m][5] = ((spline[m-2][6]-spline[m+2][6]) +
8.0*(spline[m+1][6]-spline[m-1][6])) / 12.0;
for (int m = 1; m <= n-1; m++) {
spline[m][4] = 3.0*(spline[m+1][6]-spline[m][6]) -
2.0*spline[m][5] - spline[m+1][5];
spline[m][3] = spline[m][5] + spline[m+1][5] -
2.0*(spline[m+1][6]-spline[m][6]);
}
spline[n][4] = 0.0;
spline[n][3] = 0.0;
for (int m = 1; m <= n; m++) {
spline[m][2] = spline[m][5]/delta;
spline[m][1] = 2.0*spline[m][4]/delta;
spline[m][0] = 3.0*spline[m][3]/delta;
}
}
/* ----------------------------------------------------------------------
grab global line from file and store info in setfl
return 0 if error
------------------------------------------------------------------------- */
int PairEIM::grabglobal(FILE *fptr)
{
char line[MAXLINE];
char *pch = NULL, *data = NULL;
while (pch == NULL) {
if (fgets(line,MAXLINE,fptr) == NULL) break;
pch = strstr(line,"global");
if (pch != NULL) {
data = strtok (line," \t\n\r\f");
data = strtok (NULL,"?");
sscanf(data,"%lg %lg %lg",&setfl->division,&setfl->rbig,&setfl->rsmall);
}
}
if (pch == NULL) return 0;
return 1;
}
/* ----------------------------------------------------------------------
grab elemental line from file and store info in setfl
return 0 if error
------------------------------------------------------------------------- */
int PairEIM::grabsingle(FILE *fptr, int i)
{
char line[MAXLINE];
rewind(fptr);
char *pch1 = NULL, *pch2 = NULL, *data = NULL;
while (pch1 == NULL || pch2 == NULL) {
if (fgets(line,MAXLINE,fptr) == NULL) break;
pch1 = strtok (line," \t\n\r\f");
pch1 = strstr(pch1,"element:");
if (pch1 != NULL) {
pch2 = strtok(NULL, " \t\n\r\f");
if (pch2 != NULL) {
data = strtok (NULL, "?");
if (strcmp(pch2,elements[i]) == 0) {
sscanf(data,"%d %lg %lg %lg %lg %lg %lg",&setfl->ielement[i],
&setfl->mass[i],&setfl->negativity[i],&setfl->ra[i],
&setfl->ri[i],&setfl->Ec[i],&setfl->q0[i]);
} else pch2 = NULL;
}
}
}
if (pch1 == NULL || pch2 == NULL) return 0;
return 1;
}
/* ----------------------------------------------------------------------
grab pair line from file and store info in setfl
return 0 if error
------------------------------------------------------------------------- */
int PairEIM::grabpair(FILE *fptr, int i, int j)
{
char line[MAXLINE];
rewind(fptr);
int ij;
if (i == j) ij = i;
else if (i < j) ij = nelements*(i+1) - (i+1)*(i+2)/2 + j;
else ij = nelements*(j+1) - (j+1)*(j+2)/2 + i;
char *pch1 = NULL, *pch2 = NULL, *pch3 = NULL, *data = NULL;
while (pch1 == NULL || pch2 == NULL || pch3 == NULL) {
if (fgets(line,MAXLINE,fptr) == NULL) break;
pch1 = strtok (line," \t\n\r\f");
pch1 = strstr(pch1,"pair:");
if (pch1 != NULL) {
pch2 = strtok (NULL, " \t\n\r\f");
if (pch2 != NULL) pch3 = strtok (NULL, " \t\n\r\f");
if (pch3 != NULL) data = strtok (NULL, "?");
if ((pch2 != NULL) && (pch3 != NULL)) {
if ((strcmp(pch2,elements[i]) == 0 &&
strcmp(pch3,elements[j]) == 0) ||
(strcmp(pch2,elements[j]) == 0 &&
strcmp(pch3,elements[i]) == 0)) {
sscanf(data,"%lg %lg %lg %lg %lg",
&setfl->rcutphiA[ij],&setfl->rcutphiR[ij],
&setfl->Eb[ij],&setfl->r0[ij],&setfl->alpha[ij]);
fgets(line,MAXLINE,fptr);
sscanf(line,"%lg %lg %lg %lg %lg",
&setfl->beta[ij],&setfl->rcutq[ij],&setfl->Asigma[ij],
&setfl->rq[ij],&setfl->rcutsigma[ij]);
fgets(line,MAXLINE,fptr);
sscanf(line,"%lg %lg %lg %d",
&setfl->Ac[ij],&setfl->zeta[ij],&setfl->rs[ij],
&setfl->tp[ij]);
} else {
pch1 = NULL;
pch2 = NULL;
pch3 = NULL;
}
}
}
}
if (pch1 == NULL || pch2 == NULL || pch3 == NULL) return 0;
return 1;
}
/* ----------------------------------------------------------------------
cutoff function
------------------------------------------------------------------------- */
double PairEIM::funccutoff(double rp, double rc, double r)
{
double rbig = setfl->rbig;
double rsmall = setfl->rsmall;
double a = (rsmall-rbig)/(rc-rp)*(r-rp)+rbig;
a = erfc(a);
double b = erfc(rbig);
double c = erfc(rsmall);
return (a-c)/(b-c);
}
/* ----------------------------------------------------------------------
pair interaction function phi
------------------------------------------------------------------------- */
double PairEIM::funcphi(int i, int j, double r)
{
int ij;
double value = 0.0;
if (i == j) ij = i;
else if (i < j) ij = nelements*(i+1) - (i+1)*(i+2)/2 + j;
else ij = nelements*(j+1) - (j+1)*(j+2)/2 + i;
if (r < 0.2) r = 0.2;
if (setfl->tp[ij] == 1) {
double a = setfl->Eb[ij]*setfl->alpha[ij] /
(setfl->beta[ij]-setfl->alpha[ij]);
double b = setfl->Eb[ij]*setfl->beta[ij] /
(setfl->beta[ij]-setfl->alpha[ij]);
if (r < setfl->rcutphiA[ij]) {
value -= a*exp(-setfl->beta[ij]*(r/setfl->r0[ij]-1.0))*
funccutoff(setfl->r0[ij],setfl->rcutphiA[ij],r);
}
if (r < setfl-> rcutphiR[ij]) {
value += b*exp(-setfl->alpha[ij]*(r/setfl->r0[ij]-1.0))*
funccutoff(setfl->r0[ij],setfl->rcutphiR[ij],r);
}
} else if (setfl->tp[ij] == 2) {
double a=setfl->Eb[ij]*setfl->alpha[ij]*pow(setfl->r0[ij],setfl->beta[ij])/
(setfl->beta[ij]-setfl->alpha[ij]);
double b=a*setfl->beta[ij]/setfl->alpha[ij]*
pow(setfl->r0[ij],setfl->alpha[ij]-setfl->beta[ij]);
if (r < setfl->rcutphiA[ij]) {
value -= a/pow(r,setfl->beta[ij])*
funccutoff(setfl->r0[ij],setfl->rcutphiA[ij],r);
}
if (r < setfl-> rcutphiR[ij]) {
value += b/pow(r,setfl->alpha[ij])*
funccutoff(setfl->r0[ij],setfl->rcutphiR[ij],r);
}
}
return value;
}
/* ----------------------------------------------------------------------
ion propensity function sigma
------------------------------------------------------------------------- */
double PairEIM::funcsigma(int i, int j, double r)
{
int ij;
double value = 0.0;
if (i == j) ij = i;
else if (i < j) ij = nelements*(i+1) - (i+1)*(i+2)/2 + j;
else ij = nelements*(j+1) - (j+1)*(j+2)/2 + i;
if (r < 0.2) r = 0.2;
if (r < setfl->rcutq[ij]) {
value = setfl->Asigma[ij]*(setfl->negativity[j]-setfl->negativity[i]) *
funccutoff(setfl->rq[ij],setfl->rcutq[ij],r);
}
return value;
}
/* ----------------------------------------------------------------------
charge-charge interaction function sigma
------------------------------------------------------------------------- */
double PairEIM::funccoul(int i, int j, double r)
{
int ij;
double value = 0.0;
if (i == j) ij = i;
else if (i < j) ij = nelements*(i+1) - (i+1)*(i+2)/2 + j;
else ij = nelements*(j+1) - (j+1)*(j+2)/2 + i;
if (r < 0.2) r = 0.2;
if (r < setfl->rcutsigma[ij]) {
value = setfl->Ac[ij]*exp(-setfl->zeta[ij]*r)*
funccutoff(setfl->rs[ij],setfl->rcutsigma[ij],r);
}
return value;
}
/* ---------------------------------------------------------------------- */
int PairEIM::pack_forward_comm(int n, int *list, double *buf,
int pbc_flag, int *pbc)
{
int i,j,m;
m = 0;
if (rhofp == 1) {
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = rho[j];
}
}
if (rhofp == 2) {
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = fp[j];
}
}
return m;
}
/* ---------------------------------------------------------------------- */
void PairEIM::unpack_forward_comm(int n, int first, double *buf)
{
int i,m,last;
m = 0;
last = first + n;
if (rhofp == 1) {
for (i = first; i < last; i++) rho[i] = buf[m++];
}
if (rhofp == 2) {
for (i = first; i < last; i++) fp[i] = buf[m++];
}
}
/* ---------------------------------------------------------------------- */
int PairEIM::pack_reverse_comm(int n, int first, double *buf)
{
int i,m,last;
m = 0;
last = first + n;
if (rhofp == 1) {
for (i = first; i < last; i++) buf[m++] = rho[i];
}
if (rhofp == 2) {
for (i = first; i < last; i++) buf[m++] = fp[i];
}
return m;
}
/* ---------------------------------------------------------------------- */
void PairEIM::unpack_reverse_comm(int n, int *list, double *buf)
{
int i,j,m;
m = 0;
if (rhofp == 1) {
for (i = 0; i < n; i++) {
j = list[i];
rho[j] += buf[m++];
}
}
if (rhofp == 2) {
for (i = 0; i < n; i++) {
j = list[i];
fp[j] += buf[m++];
}
}
}
/* ----------------------------------------------------------------------
memory usage of local atom-based arrays
------------------------------------------------------------------------- */
double PairEIM::memory_usage()
{
double bytes = maxeatom * sizeof(double);
bytes += maxvatom*6 * sizeof(double);
bytes += 2 * nmax * sizeof(double);
return bytes;
}
Event Timeline
Log In to Comment