jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

remove REAX and MEAM from src folder

This commit is contained in:
Axel Kohlmeyer 2018-12-10 12:07:51 -05:00
parent 8b5887bfb1
commit 054c525842
10 changed files with 0 additions and 3028 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

67
src/MEAM/Install.sh
View File

@ -1,67 +0,0 @@
# Install/unInstall package files in LAMMPS
# mode = 0/1/2 for uninstall/install/update
mode=$1
# enforce using portable C locale
LC_ALL=C
export LC_ALL
# arg1 = file, arg2 = file it depends on
action () {
if (test $mode = 0) then
rm -f ../$1
elif (! cmp -s $1 ../$1) then
if (test -z "$2" || test -e ../$2) then
cp $1 ..
if (test $mode = 2) then
echo " updating src/$1"
fi
fi
elif (test -n "$2") then
if (test ! -e ../$2) then
rm -f ../$1
fi
fi
}
# all package files with no dependencies
for file in *.cpp *.h; do
test -f ${file} && action $file
done
# edit 2 Makefile.package files to include/exclude package info
if (test $1 = 1) then
if (test -e ../Makefile.package) then
sed -i -e 's/[^ \t]*meam[^ \t]* //' ../Makefile.package
sed -i -e 's|^PKG_INC =[ \t]*|&-I../../lib/meam |' ../Makefile.package
sed -i -e 's|^PKG_PATH =[ \t]*|&-L../../lib/meam |' ../Makefile.package
sed -i -e 's|^PKG_LIB =[ \t]*|&-lmeam |' ../Makefile.package
sed -i -e 's|^PKG_SYSINC =[ \t]*|&$(meam_SYSINC) |' ../Makefile.package
sed -i -e 's|^PKG_SYSLIB =[ \t]*|&$(meam_SYSLIB) |' ../Makefile.package
sed -i -e 's|^PKG_SYSPATH =[ \t]*|&$(meam_SYSPATH) |' ../Makefile.package
fi
if (test -e ../Makefile.package.settings) then
sed -i -e '/^include.*meam.*$/d' ../Makefile.package.settings
# multiline form needed for BSD sed on Macs
sed -i -e '4 i \
include ..\/..\/lib\/meam\/Makefile.lammps
' ../Makefile.package.settings
fi
elif (test $1 = 0) then
if (test -e ../Makefile.package) then
sed -i -e 's/[^ \t]*meam[^ \t]* //' ../Makefile.package
fi
if (test -e ../Makefile.package.settings) then
sed -i -e '/^include.*meam.*$/d' ../Makefile.package.settings
fi
fi

11
src/MEAM/README
View File

@ -1,11 +0,0 @@
This package implements the MEAM potential as a LAMMPS pair style.
This package uses an external library in lib/meam which must be
compiled before making LAMMPS. See the lib/meam/README file and the
LAMMPS manual for information on building LAMMPS with external
libraries. The settings in the Makefile.lammps file in that directory
must be correct for LAMMPS to build correctly with this package
installed.
The person who created the MEAM library and this pair style is Greg
Wagner (while at Sandia, now at Northwestern U).

949
src/MEAM/pair_meam.cpp
View File

@ -1,949 +0,0 @@
/* ----------------------------------------------------------------------
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: Greg Wagner (SNL)
------------------------------------------------------------------------- */
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include "pair_meam.h"
#include "atom.h"
#include "force.h"
#include "comm.h"
#include "memory.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
#define MAXLINE 1024
enum{FCC,BCC,HCP,DIM,DIAMOND,B1,C11,L12,B2};
static const int nkeywords = 21;
static const char *keywords[] = {
"Ec","alpha","rho0","delta","lattce",
"attrac","repuls","nn2","Cmin","Cmax","rc","delr",
"augt1","gsmooth_factor","re","ialloy",
"mixture_ref_t","erose_form","zbl",
"emb_lin_neg","bkgd_dyn"};
/* ---------------------------------------------------------------------- */
PairMEAM::PairMEAM(LAMMPS *lmp) : Pair(lmp)
{
if (comm->me == 0)
error->warning(FLERR,"THE pair_style meam COMMAND IS OBSOLETE AND "
"WILL BE REMOVED VERY SOON. PLEASE USE meam/c");
single_enable = 0;
restartinfo = 0;
one_coeff = 1;
manybody_flag = 1;
nmax = 0;
rho = rho0 = rho1 = rho2 = rho3 = frhop = NULL;
gamma = dgamma1 = dgamma2 = dgamma3 = arho2b = NULL;
arho1 = arho2 = arho3 = arho3b = t_ave = tsq_ave = NULL;
maxneigh = 0;
allocated = 0;
scrfcn = dscrfcn = fcpair = NULL;
nelements = 0;
elements = NULL;
mass = NULL;
// set comm size needed by this Pair
comm_forward = 38;
comm_reverse = 30;
}
/* ----------------------------------------------------------------------
free all arrays
check if allocated, since class can be destructed when incomplete
------------------------------------------------------------------------- */
PairMEAM::~PairMEAM()
{
meam_cleanup_();
memory->destroy(rho);
memory->destroy(rho0);
memory->destroy(rho1);
memory->destroy(rho2);
memory->destroy(rho3);
memory->destroy(frhop);
memory->destroy(gamma);
memory->destroy(dgamma1);
memory->destroy(dgamma2);
memory->destroy(dgamma3);
memory->destroy(arho2b);
memory->destroy(arho1);
memory->destroy(arho2);
memory->destroy(arho3);
memory->destroy(arho3b);
memory->destroy(t_ave);
memory->destroy(tsq_ave);
memory->destroy(scrfcn);
memory->destroy(dscrfcn);
memory->destroy(fcpair);
for (int i = 0; i < nelements; i++) delete [] elements[i];
delete [] elements;
delete [] mass;
if (allocated) {
memory->destroy(setflag);
memory->destroy(cutsq);
delete [] map;
delete [] fmap;
}
}
/* ---------------------------------------------------------------------- */
void PairMEAM::compute(int eflag, int vflag)
{
int i,j,ii,n,inum_half,errorflag;
int *ilist_half,*numneigh_half,**firstneigh_half;
int *numneigh_full,**firstneigh_full;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = eflag_global = vflag_global =
eflag_atom = vflag_atom = 0;
// grow local arrays if necessary
if (atom->nmax > nmax) {
memory->destroy(rho);
memory->destroy(rho0);
memory->destroy(rho1);
memory->destroy(rho2);
memory->destroy(rho3);
memory->destroy(frhop);
memory->destroy(gamma);
memory->destroy(dgamma1);
memory->destroy(dgamma2);
memory->destroy(dgamma3);
memory->destroy(arho2b);
memory->destroy(arho1);
memory->destroy(arho2);
memory->destroy(arho3);
memory->destroy(arho3b);
memory->destroy(t_ave);
memory->destroy(tsq_ave);
nmax = atom->nmax;
memory->create(rho,nmax,"pair:rho");
memory->create(rho0,nmax,"pair:rho0");
memory->create(rho1,nmax,"pair:rho1");
memory->create(rho2,nmax,"pair:rho2");
memory->create(rho3,nmax,"pair:rho3");
memory->create(frhop,nmax,"pair:frhop");
memory->create(gamma,nmax,"pair:gamma");
memory->create(dgamma1,nmax,"pair:dgamma1");
memory->create(dgamma2,nmax,"pair:dgamma2");
memory->create(dgamma3,nmax,"pair:dgamma3");
memory->create(arho2b,nmax,"pair:arho2b");
memory->create(arho1,nmax,3,"pair:arho1");
memory->create(arho2,nmax,6,"pair:arho2");
memory->create(arho3,nmax,10,"pair:arho3");
memory->create(arho3b,nmax,3,"pair:arho3b");
memory->create(t_ave,nmax,3,"pair:t_ave");
memory->create(tsq_ave,nmax,3,"pair:tsq_ave");
}
// neighbor list info
inum_half = listhalf->inum;
ilist_half = listhalf->ilist;
numneigh_half = listhalf->numneigh;
firstneigh_half = listhalf->firstneigh;
numneigh_full = listfull->numneigh;
firstneigh_full = listfull->firstneigh;
// strip neighbor lists of any special bond flags before using with MEAM
// necessary before doing neigh_f2c and neigh_c2f conversions each step
if (neighbor->ago == 0) {
neigh_strip(inum_half,ilist_half,numneigh_half,firstneigh_half);
neigh_strip(inum_half,ilist_half,numneigh_full,firstneigh_full);
}
// check size of scrfcn based on half neighbor list
int nlocal = atom->nlocal;
int nall = nlocal + atom->nghost;
n = 0;
for (ii = 0; ii < inum_half; ii++) n += numneigh_half[ilist_half[ii]];
if (n > maxneigh) {
memory->destroy(scrfcn);
memory->destroy(dscrfcn);
memory->destroy(fcpair);
maxneigh = n;
memory->create(scrfcn,maxneigh,"pair:scrfcn");
memory->create(dscrfcn,maxneigh,"pair:dscrfcn");
memory->create(fcpair,maxneigh,"pair:fcpair");
}
// zero out local arrays
for (i = 0; i < nall; i++) {
rho0[i] = 0.0;
arho2b[i] = 0.0;
arho1[i][0] = arho1[i][1] = arho1[i][2] = 0.0;
for (j = 0; j < 6; j++) arho2[i][j] = 0.0;
for (j = 0; j < 10; j++) arho3[i][j] = 0.0;
arho3b[i][0] = arho3b[i][1] = arho3b[i][2] = 0.0;
t_ave[i][0] = t_ave[i][1] = t_ave[i][2] = 0.0;
tsq_ave[i][0] = tsq_ave[i][1] = tsq_ave[i][2] = 0.0;
}
double **x = atom->x;
double **f = atom->f;
int *type = atom->type;
int ntype = atom->ntypes;
// change neighbor list indices to Fortran indexing
neigh_c2f(inum_half,ilist_half,numneigh_half,firstneigh_half);
neigh_c2f(inum_half,ilist_half,numneigh_full,firstneigh_full);
// 3 stages of MEAM calculation
// loop over my atoms followed by communication
int ifort;
int offset = 0;
errorflag = 0;
for (ii = 0; ii < inum_half; ii++) {
i = ilist_half[ii];
ifort = i+1;
meam_dens_init_(&ifort,&nmax,&ntype,type,fmap,&x[0][0],
&numneigh_half[i],firstneigh_half[i],
&numneigh_full[i],firstneigh_full[i],
&scrfcn[offset],&dscrfcn[offset],&fcpair[offset],
rho0,&arho1[0][0],&arho2[0][0],arho2b,
&arho3[0][0],&arho3b[0][0],&t_ave[0][0],&tsq_ave[0][0],
&errorflag);
if (errorflag) {
char str[128];
sprintf(str,"MEAM library error %d",errorflag);
error->one(FLERR,str);
}
offset += numneigh_half[i];
}
comm->reverse_comm_pair(this);
meam_dens_final_(&nlocal,&nmax,&eflag_either,&eflag_global,&eflag_atom,
&eng_vdwl,eatom,&ntype,type,fmap,
&arho1[0][0],&arho2[0][0],arho2b,&arho3[0][0],
&arho3b[0][0],&t_ave[0][0],&tsq_ave[0][0],gamma,dgamma1,
dgamma2,dgamma3,rho,rho0,rho1,rho2,rho3,frhop,&errorflag);
if (errorflag) {
char str[128];
sprintf(str,"MEAM library error %d",errorflag);
error->one(FLERR,str);
}
comm->forward_comm_pair(this);
offset = 0;
// vptr is first value in vatom if it will be used by meam_force()
// else vatom may not exist, so pass dummy ptr
double *vptr;
if (vflag_atom) vptr = &vatom[0][0];
else vptr = &cutmax;
for (ii = 0; ii < inum_half; ii++) {
i = ilist_half[ii];
ifort = i+1;
meam_force_(&ifort,&nmax,&eflag_either,&eflag_global,&eflag_atom,
&vflag_atom,&eng_vdwl,eatom,&ntype,type,fmap,&x[0][0],
&numneigh_half[i],firstneigh_half[i],
&numneigh_full[i],firstneigh_full[i],
&scrfcn[offset],&dscrfcn[offset],&fcpair[offset],
dgamma1,dgamma2,dgamma3,rho0,rho1,rho2,rho3,frhop,
&arho1[0][0],&arho2[0][0],arho2b,&arho3[0][0],&arho3b[0][0],
&t_ave[0][0],&tsq_ave[0][0],&f[0][0],vptr,&errorflag);
if (errorflag) {
char str[128];
sprintf(str,"MEAM library error %d",errorflag);
error->one(FLERR,str);
}
offset += numneigh_half[i];
}
// change neighbor list indices back to C indexing
neigh_f2c(inum_half,ilist_half,numneigh_half,firstneigh_half);
neigh_f2c(inum_half,ilist_half,numneigh_full,firstneigh_full);
if (vflag_fdotr) virial_fdotr_compute();
}
/* ---------------------------------------------------------------------- */
void PairMEAM::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");
map = new int[n+1];
fmap = new int[n];
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairMEAM::settings(int narg, char **/*arg*/)
{
if (narg != 0) error->all(FLERR,"Illegal pair_style command");
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairMEAM::coeff(int narg, char **arg)
{
int i,j,m,n;
if (!allocated) allocate();
if (narg < 6) 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 MEAM element names between 2 filenames
// nelements = # of MEAM elements
// elements = list of unique element names
if (nelements) {
for (i = 0; i < nelements; i++) delete [] elements[i];
delete [] elements;
delete [] mass;
}
nelements = narg - 4 - atom->ntypes;
if (nelements < 1) error->all(FLERR,"Incorrect args for pair coefficients");
elements = new char*[nelements];
mass = new double[nelements];
for (i = 0; i < nelements; i++) {
n = strlen(arg[i+3]) + 1;
elements[i] = new char[n];
strcpy(elements[i],arg[i+3]);
}
// read MEAM library and parameter files
// pass all parameters to MEAM package
// tell MEAM package that setup is done
read_files(arg[2],arg[2+nelements+1]);
meam_setup_done_(&cutmax);
// read args that map atom types to MEAM elements
// map[i] = which element the Ith atom type is, -1 if not mapped
for (i = 4 + nelements; i < narg; i++) {
m = i - (4+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 (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
setflag[i][j] = 0;
// set setflag i,j for type pairs where both are mapped to elements
// set mass for i,i in atom class
int count = 0;
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
if (map[i] >= 0 && map[j] >= 0) {
setflag[i][j] = 1;
if (i == j) atom->set_mass(FLERR,i,mass[map[i]]);
count++;
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void PairMEAM::init_style()
{
if (force->newton_pair == 0)
error->all(FLERR,"Pair style MEAM requires newton pair on");
// need full and half neighbor list
int irequest_full = neighbor->request(this,instance_me);
neighbor->requests[irequest_full]->id = 1;
neighbor->requests[irequest_full]->half = 0;
neighbor->requests[irequest_full]->full = 1;
int irequest_half = neighbor->request(this,instance_me);
neighbor->requests[irequest_half]->id = 2;
// setup Fortran-style mapping array needed by MEAM package
// fmap is indexed from 1:ntypes by Fortran and stores a Fortran index
// if type I is not a MEAM atom, fmap stores a 0
for (int i = 1; i <= atom->ntypes; i++) fmap[i-1] = map[i] + 1;
}
/* ----------------------------------------------------------------------
neighbor callback to inform pair style of neighbor list to use
half or full
------------------------------------------------------------------------- */
void PairMEAM::init_list(int id, NeighList *ptr)
{
if (id == 1) listfull = ptr;
else if (id == 2) listhalf = ptr;
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairMEAM::init_one(int /*i*/, int /*j*/)
{
return cutmax;
}
/* ---------------------------------------------------------------------- */
void PairMEAM::read_files(char *globalfile, char *userfile)
{
// open global meamf file on proc 0
FILE *fp;
if (comm->me == 0) {
fp = force->open_potential(globalfile);
if (fp == NULL) {
char str[128];
snprintf(str,128,"Cannot open MEAM potential file %s",globalfile);
error->one(FLERR,str);
}
}
// allocate parameter arrays
int params_per_line = 19;
int *lat = new int[nelements];
int *ielement = new int[nelements];
int *ibar = new int[nelements];
double *z = new double[nelements];
double *atwt = new double[nelements];
double *alpha = new double[nelements];
double *b0 = new double[nelements];
double *b1 = new double[nelements];
double *b2 = new double[nelements];
double *b3 = new double[nelements];
double *alat = new double[nelements];
double *esub = new double[nelements];
double *asub = new double[nelements];
double *t0 = new double[nelements];
double *t1 = new double[nelements];
double *t2 = new double[nelements];
double *t3 = new double[nelements];
double *rozero = new double[nelements];
bool *found = new bool[nelements];
for (int i = 0; i < nelements; i++) found[i] = false;
// read each set of params from global MEAM file
// one set of params can span multiple lines
// store params if element name is in element list
// if element name appears multiple times, only store 1st entry
int i,n,nwords;
char **words = new char*[params_per_line+1];
char line[MAXLINE],*ptr;
int eof = 0;
int nset = 0;
while (1) {
if (comm->me == 0) {
ptr = fgets(line,MAXLINE,fp);
if (ptr == NULL) {
eof = 1;
fclose(fp);
} else n = strlen(line) + 1;
}
MPI_Bcast(&eof,1,MPI_INT,0,world);
if (eof) break;
MPI_Bcast(&n,1,MPI_INT,0,world);
MPI_Bcast(line,n,MPI_CHAR,0,world);
// strip comment, skip line if blank
if ((ptr = strchr(line,'#'))) *ptr = '\0';
nwords = atom->count_words(line);
if (nwords == 0) continue;
// concatenate additional lines until have params_per_line words
while (nwords < params_per_line) {
n = strlen(line);
if (comm->me == 0) {
ptr = fgets(&line[n],MAXLINE-n,fp);
if (ptr == NULL) {
eof = 1;
fclose(fp);
} else n = strlen(line) + 1;
}
MPI_Bcast(&eof,1,MPI_INT,0,world);
if (eof) break;
MPI_Bcast(&n,1,MPI_INT,0,world);
MPI_Bcast(line,n,MPI_CHAR,0,world);
if ((ptr = strchr(line,'#'))) *ptr = '\0';
nwords = atom->count_words(line);
}
if (nwords != params_per_line)
error->all(FLERR,"Incorrect format in MEAM potential file");
// words = ptrs to all words in line
// strip single and double quotes from words
nwords = 0;
words[nwords++] = strtok(line,"' \t\n\r\f");
while ((words[nwords++] = strtok(NULL,"' \t\n\r\f"))) continue;
// skip if element name isn't in element list
for (i = 0; i < nelements; i++)
if (strcmp(words[0],elements[i]) == 0) break;
if (i >= nelements) continue;
// skip if element already appeared
if (found[i] == true) continue;
found[i] = true;
// map lat string to an integer
if (strcmp(words[1],"fcc") == 0) lat[i] = FCC;
else if (strcmp(words[1],"bcc") == 0) lat[i] = BCC;
else if (strcmp(words[1],"hcp") == 0) lat[i] = HCP;
else if (strcmp(words[1],"dim") == 0) lat[i] = DIM;
else if (strcmp(words[1],"dia") == 0) lat[i] = DIAMOND;
else error->all(FLERR,"Unrecognized lattice type in MEAM file 1");
// store parameters
z[i] = atof(words[2]);
ielement[i] = atoi(words[3]);
atwt[i] = atof(words[4]);
alpha[i] = atof(words[5]);
b0[i] = atof(words[6]);
b1[i] = atof(words[7]);
b2[i] = atof(words[8]);
b3[i] = atof(words[9]);
alat[i] = atof(words[10]);
esub[i] = atof(words[11]);
asub[i] = atof(words[12]);
t0[i] = atof(words[13]);
t1[i] = atof(words[14]);
t2[i] = atof(words[15]);
t3[i] = atof(words[16]);
rozero[i] = atof(words[17]);
ibar[i] = atoi(words[18]);
nset++;
}
// error if didn't find all elements in file
if (nset != nelements)
error->all(FLERR,"Did not find all elements in MEAM library file");
// pass element parameters to MEAM package
meam_setup_global_(&nelements,lat,z,ielement,atwt,alpha,b0,b1,b2,b3,
alat,esub,asub,t0,t1,t2,t3,rozero,ibar);
// set element masses
for (i = 0; i < nelements; i++) mass[i] = atwt[i];
// clean-up memory
delete [] words;
delete [] lat;
delete [] ielement;
delete [] ibar;
delete [] z;
delete [] atwt;
delete [] alpha;
delete [] b0;
delete [] b1;
delete [] b2;
delete [] b3;
delete [] alat;
delete [] esub;
delete [] asub;
delete [] t0;
delete [] t1;
delete [] t2;
delete [] t3;
delete [] rozero;
delete [] found;
// done if user param file is NULL
if (strcmp(userfile,"NULL") == 0) return;
// open user param file on proc 0
if (comm->me == 0) {
fp = force->open_potential(userfile);
if (fp == NULL) {
char str[128];
snprintf(str,128,"Cannot open MEAM potential file %s",userfile);
error->one(FLERR,str);
}
}
// read settings
// pass them one at a time to MEAM package
// match strings to list of corresponding ints
int which;
double value;
int nindex,index[3];
int maxparams = 6;
char **params = new char*[maxparams];
int nparams;
eof = 0;
while (1) {
if (comm->me == 0) {
ptr = fgets(line,MAXLINE,fp);
if (ptr == NULL) {
eof = 1;
fclose(fp);
} else n = strlen(line) + 1;
}
MPI_Bcast(&eof,1,MPI_INT,0,world);
if (eof) break;
MPI_Bcast(&n,1,MPI_INT,0,world);
MPI_Bcast(line,n,MPI_CHAR,0,world);
// strip comment, skip line if blank
if ((ptr = strchr(line,'#'))) *ptr = '\0';
nparams = atom->count_words(line);
if (nparams == 0) continue;
// words = ptrs to all words in line
nparams = 0;
params[nparams++] = strtok(line,"=(), '\t\n\r\f");
while (nparams < maxparams &&
(params[nparams++] = strtok(NULL,"=(), '\t\n\r\f")))
continue;
nparams--;
for (which = 0; which < nkeywords; which++)
if (strcmp(params[0],keywords[which]) == 0) break;
if (which == nkeywords) {
char str[128];
snprintf(str,128,"Keyword %s in MEAM parameter file not recognized",
params[0]);
error->all(FLERR,str);
}
nindex = nparams - 2;
for (i = 0; i < nindex; i++) index[i] = atoi(params[i+1]);
// map lattce_meam value to an integer
if (which == 4) {
if (strcmp(params[nparams-1],"fcc") == 0) value = FCC;
else if (strcmp(params[nparams-1],"bcc") == 0) value = BCC;
else if (strcmp(params[nparams-1],"hcp") == 0) value = HCP;
else if (strcmp(params[nparams-1],"dim") == 0) value = DIM;
else if (strcmp(params[nparams-1],"dia") == 0) value = DIAMOND;
else if (strcmp(params[nparams-1],"b1") == 0) value = B1;
else if (strcmp(params[nparams-1],"c11") == 0) value = C11;
else if (strcmp(params[nparams-1],"l12") == 0) value = L12;
else if (strcmp(params[nparams-1],"b2") == 0) value = B2;
else error->all(FLERR,"Unrecognized lattice type in MEAM file 2");
}
else value = atof(params[nparams-1]);
// pass single setting to MEAM package
int errorflag = 0;
meam_setup_param_(&which,&value,&nindex,index,&errorflag);
if (errorflag) {
char str[128];
sprintf(str,"MEAM library error %d",errorflag);
error->all(FLERR,str);
}
}
delete [] params;
}
/* ---------------------------------------------------------------------- */
int PairMEAM::pack_forward_comm(int n, int *list, double *buf,
int /*pbc_flag*/, int * /*pbc*/)
{
int i,j,k,m;
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = rho0[j];
buf[m++] = rho1[j];
buf[m++] = rho2[j];
buf[m++] = rho3[j];
buf[m++] = frhop[j];
buf[m++] = gamma[j];
buf[m++] = dgamma1[j];
buf[m++] = dgamma2[j];
buf[m++] = dgamma3[j];
buf[m++] = arho2b[j];
buf[m++] = arho1[j][0];
buf[m++] = arho1[j][1];
buf[m++] = arho1[j][2];
buf[m++] = arho2[j][0];
buf[m++] = arho2[j][1];
buf[m++] = arho2[j][2];
buf[m++] = arho2[j][3];
buf[m++] = arho2[j][4];
buf[m++] = arho2[j][5];
for (k = 0; k < 10; k++) buf[m++] = arho3[j][k];
buf[m++] = arho3b[j][0];
buf[m++] = arho3b[j][1];
buf[m++] = arho3b[j][2];
buf[m++] = t_ave[j][0];
buf[m++] = t_ave[j][1];
buf[m++] = t_ave[j][2];
buf[m++] = tsq_ave[j][0];
buf[m++] = tsq_ave[j][1];
buf[m++] = tsq_ave[j][2];
}
return m;
}
/* ---------------------------------------------------------------------- */
void PairMEAM::unpack_forward_comm(int n, int first, double *buf)
{
int i,k,m,last;
m = 0;
last = first + n;
for (i = first; i < last; i++) {
rho0[i] = buf[m++];
rho1[i] = buf[m++];
rho2[i] = buf[m++];
rho3[i] = buf[m++];
frhop[i] = buf[m++];
gamma[i] = buf[m++];
dgamma1[i] = buf[m++];
dgamma2[i] = buf[m++];
dgamma3[i] = buf[m++];
arho2b[i] = buf[m++];
arho1[i][0] = buf[m++];
arho1[i][1] = buf[m++];
arho1[i][2] = buf[m++];
arho2[i][0] = buf[m++];
arho2[i][1] = buf[m++];
arho2[i][2] = buf[m++];
arho2[i][3] = buf[m++];
arho2[i][4] = buf[m++];
arho2[i][5] = buf[m++];
for (k = 0; k < 10; k++) arho3[i][k] = buf[m++];
arho3b[i][0] = buf[m++];
arho3b[i][1] = buf[m++];
arho3b[i][2] = buf[m++];
t_ave[i][0] = buf[m++];
t_ave[i][1] = buf[m++];
t_ave[i][2] = buf[m++];
tsq_ave[i][0] = buf[m++];
tsq_ave[i][1] = buf[m++];
tsq_ave[i][2] = buf[m++];
}
}
/* ---------------------------------------------------------------------- */
int PairMEAM::pack_reverse_comm(int n, int first, double *buf)
{
int i,k,m,last;
m = 0;
last = first + n;
for (i = first; i < last; i++) {
buf[m++] = rho0[i];
buf[m++] = arho2b[i];
buf[m++] = arho1[i][0];
buf[m++] = arho1[i][1];
buf[m++] = arho1[i][2];
buf[m++] = arho2[i][0];
buf[m++] = arho2[i][1];
buf[m++] = arho2[i][2];
buf[m++] = arho2[i][3];
buf[m++] = arho2[i][4];
buf[m++] = arho2[i][5];
for (k = 0; k < 10; k++) buf[m++] = arho3[i][k];
buf[m++] = arho3b[i][0];
buf[m++] = arho3b[i][1];
buf[m++] = arho3b[i][2];
buf[m++] = t_ave[i][0];
buf[m++] = t_ave[i][1];
buf[m++] = t_ave[i][2];
buf[m++] = tsq_ave[i][0];
buf[m++] = tsq_ave[i][1];
buf[m++] = tsq_ave[i][2];
}
return m;
}
/* ---------------------------------------------------------------------- */
void PairMEAM::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];
rho0[j] += buf[m++];
arho2b[j] += buf[m++];
arho1[j][0] += buf[m++];
arho1[j][1] += buf[m++];
arho1[j][2] += buf[m++];
arho2[j][0] += buf[m++];
arho2[j][1] += buf[m++];
arho2[j][2] += buf[m++];
arho2[j][3] += buf[m++];
arho2[j][4] += buf[m++];
arho2[j][5] += buf[m++];
for (k = 0; k < 10; k++) arho3[j][k] += buf[m++];
arho3b[j][0] += buf[m++];
arho3b[j][1] += buf[m++];
arho3b[j][2] += buf[m++];
t_ave[j][0] += buf[m++];
t_ave[j][1] += buf[m++];
t_ave[j][2] += buf[m++];
tsq_ave[j][0] += buf[m++];
tsq_ave[j][1] += buf[m++];
tsq_ave[j][2] += buf[m++];
}
}
/* ----------------------------------------------------------------------
memory usage of local atom-based arrays
------------------------------------------------------------------------- */
double PairMEAM::memory_usage()
{
double bytes = 11 * nmax * sizeof(double);
bytes += (3 + 6 + 10 + 3 + 3 + 3) * nmax * sizeof(double);
bytes += 3 * maxneigh * sizeof(double);
return bytes;
}
/* ----------------------------------------------------------------------
strip special bond flags from neighbor list entries
are not used with MEAM
need to do here so Fortran lib doesn't see them
done once per reneighbor so that neigh_f2c and neigh_c2f don't see them
------------------------------------------------------------------------- */
void PairMEAM::neigh_strip(int inum, int *ilist,
int *numneigh, int **firstneigh)
{
int i,j,ii,jnum;
int *jlist;
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
jlist = firstneigh[i];
jnum = numneigh[i];
for (j = 0; j < jnum; j++) jlist[j] &= NEIGHMASK;
}
}
/* ----------------------------------------------------------------------
toggle neighbor list indices between zero- and one-based values
needed for access by MEAM Fortran library
------------------------------------------------------------------------- */
void PairMEAM::neigh_f2c(int inum, int *ilist, int *numneigh, int **firstneigh)
{
int i,j,ii,jnum;
int *jlist;
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
jlist = firstneigh[i];
jnum = numneigh[i];
for (j = 0; j < jnum; j++) jlist[j]--;
}
}
void PairMEAM::neigh_c2f(int inum, int *ilist, int *numneigh, int **firstneigh)
{
int i,j,ii,jnum;
int *jlist;
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
jlist = firstneigh[i];
jnum = numneigh[i];
for (j = 0; j < jnum; j++) jlist[j]++;
}
}

159
src/MEAM/pair_meam.h
View File

@ -1,159 +0,0 @@
/* -*- c++ -*- ----------------------------------------------------------
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.
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
PairStyle(meam,PairMEAM)
#else
#ifndef LMP_PAIR_MEAM_H
#define LMP_PAIR_MEAM_H
extern "C" {
void meam_setup_global_(int *, int *, double *, int *, double *, double *,
double *, double *, double *, double *, double *,
double *, double *, double *, double *, double *,
double *, double *, int *);
void meam_setup_param_(int *, double *, int *, int *, int *);
void meam_setup_done_(double *);
void meam_dens_init_(int *, int *, int *, int *, int *,
double *, int *, int *, int *, int *,
double *, double *, double *, double *,
double *, double *,
double *, double *, double *, double *, double *,
int *);
void meam_dens_final_(int *, int *, int *, int *, int *, double *, double *,
int *, int *, int *,
double *, double *, double *, double *,
double *, double *, double *,
double *, double *, double *, double *,
double *, double *,
double *, double *, double *, double *, int *);
void meam_force_(int *, int *, int *, int *, int *, int *,
double *, double *, int *, int *, int *,
double *, int *, int *, int *, int *, double *, double *,
double *, double *, double *, double *, double *, double *,
double *, double *, double *, double *, double *, double *,
double *, double *, double *, double *, double *, double *, int *);
void meam_cleanup_();
}
#include "pair.h"
namespace LAMMPS_NS {
class PairMEAM : public Pair {
public:
PairMEAM(class LAMMPS *);
~PairMEAM();
void compute(int, int);
void settings(int, char **);
void coeff(int, char **);
void init_style();
void init_list(int, class NeighList *);
double init_one(int, int);
int pack_forward_comm(int, int *, double *, int, int *);
void unpack_forward_comm(int, int, double *);
int pack_reverse_comm(int, int, double *);
void unpack_reverse_comm(int, int *, double *);
double memory_usage();
private:
double cutmax; // max cutoff for all elements
int nelements; // # of unique elements
char **elements; // names of unique elements
double *mass; // mass of each element
int *map; // mapping from atom types to elements
int *fmap; // Fortran version of map array for MEAM lib
int maxneigh;
double *scrfcn,*dscrfcn,*fcpair;
int nmax;
double *rho,*rho0,*rho1,*rho2,*rho3,*frhop;
double *gamma,*dgamma1,*dgamma2,*dgamma3,*arho2b;
double **arho1,**arho2,**arho3,**arho3b,**t_ave,**tsq_ave;
void allocate();
void read_files(char *, char *);
void neigh_strip(int, int *, int *, int **);
void neigh_f2c(int, int *, int *, int **);
void neigh_c2f(int, int *, int *, int **);
};
}
#endif
#endif
/* ERROR/WARNING messages:
W: The pair_style meam command is unsupported. Please use pair_style meam/c instead
UNDOCUMENTED
E: MEAM library error %d
A call to the MEAM Fortran library returned an error.
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
E: Pair style MEAM requires newton pair on
See the newton command. This is a restriction to use the MEAM
potential.
E: Cannot open MEAM potential file %s
The specified MEAM potential file cannot be opened. Check that the
path and name are correct.
E: Incorrect format in MEAM potential file
Incorrect number of words per line in the potential file.
E: Unrecognized lattice type in MEAM file 1
The lattice type in an entry of the MEAM library file is not
valid.
E: Did not find all elements in MEAM library file
The requested elements were not found in the MEAM file.
E: Keyword %s in MEAM parameter file not recognized
Self-explanatory.
E: Unrecognized lattice type in MEAM file 2
The lattice type in an entry of the MEAM parameter file is not
valid.
*/

67
src/REAX/Install.sh
View File

@ -1,67 +0,0 @@
# Install/unInstall package files in LAMMPS
# mode = 0/1/2 for uninstall/install/update
mode=$1
# arg1 = file, arg2 = file it depends on
# enforce using portable C locale
LC_ALL=C
export LC_ALL
action () {
if (test $mode = 0) then
rm -f ../$1
elif (! cmp -s $1 ../$1) then
if (test -z "$2" || test -e ../$2) then
cp $1 ..
if (test $mode = 2) then
echo " updating src/$1"
fi
fi
elif (test -n "$2") then
if (test ! -e ../$2) then
rm -f ../$1
fi
fi
}
# all package files with no dependencies
for file in *.cpp *.h; do
test -f ${file} && action $file
done
# edit 2 Makefile.package files to include/exclude package info
if (test $1 = 1) then
if (test -e ../Makefile.package) then
sed -i -e 's/[^ \t]*reax[^ \t]* //' ../Makefile.package
sed -i -e 's|^PKG_INC =[ \t]*|&-I../../lib/reax |' ../Makefile.package
sed -i -e 's|^PKG_PATH =[ \t]*|&-L../../lib/reax |' ../Makefile.package
sed -i -e 's|^PKG_LIB =[ \t]*|&-lreax |' ../Makefile.package
sed -i -e 's|^PKG_SYSINC =[ \t]*|&$(reax_SYSINC) |' ../Makefile.package
sed -i -e 's|^PKG_SYSLIB =[ \t]*|&$(reax_SYSLIB) |' ../Makefile.package
sed -i -e 's|^PKG_SYSPATH =[ \t]*|&$(reax_SYSPATH) |' ../Makefile.package
fi
if (test -e ../Makefile.package.settings) then
sed -i -e '/^include.*reax.*$/d' ../Makefile.package.settings
# multiline form needed for BSD sed on Macs
sed -i -e '4 i \
include ..\/..\/lib\/reax\/Makefile.lammps
' ../Makefile.package.settings
fi
elif (test $1 = 0) then
if (test -e ../Makefile.package) then
sed -i -e 's/[^ \t]*reax[^ \t]* //' ../Makefile.package
fi
if (test -e ../Makefile.package.settings) then
sed -i -e '/^include.*reax.*$/d' ../Makefile.package.settings
fi
fi

267
src/REAX/fix_reax_bonds.cpp
View File

@ -1,267 +0,0 @@
/* ----------------------------------------------------------------------
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: Aidan Thompson (Sandia)
------------------------------------------------------------------------- */
#ifdef LAMMPS_BIGBIG
#error LAMMPS_BIGBIG is not supported by the REAX package
#endif
#include <cstdlib>
#include <cstring>
#include "fix_reax_bonds.h"
#include "pair_reax_fortran.h"
#include "atom.h"
#include "update.h"
#include "force.h"
#include "modify.h"
#include "compute.h"
#include "input.h"
#include "variable.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
using namespace FixConst;
/* ---------------------------------------------------------------------- */
FixReaxBonds::FixReaxBonds(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
if (narg < 5) error->all(FLERR,"Illegal fix reax/bonds command");
MPI_Comm_rank(world,&me);
nevery = force->inumeric(FLERR,arg[3]);
if (nevery < 1) error->all(FLERR,"Illegal fix reax/bonds command");
if (me == 0) {
char *suffix = strrchr(arg[4],'.');
if (suffix && strcmp(suffix,".gz") == 0) {
#ifdef LAMMPS_GZIP
char gzip[128];
snprintf(gzip,128,"gzip -6 > %s",arg[4]);
#ifdef _WIN32
fp = _popen(gzip,"wb");
#else
fp = popen(gzip,"w");
#endif
#else
error->one(FLERR,"Cannot open gzipped file");
#endif
} else fp = fopen(arg[4],"w");
if (fp == NULL) {
char str[128];
snprintf(str,128,"Cannot open fix reax/bonds file %s",arg[4]);
error->one(FLERR,str);
}
}
}
/* ---------------------------------------------------------------------- */
FixReaxBonds::~FixReaxBonds()
{
if (me == 0) fclose(fp);
}
/* ---------------------------------------------------------------------- */
int FixReaxBonds::setmask()
{
int mask = 0;
mask |= END_OF_STEP;
return mask;
}
/* ----------------------------------------------------------------------
perform initial write
------------------------------------------------------------------------- */
void FixReaxBonds::setup(int vflag)
{
end_of_step();
}
/* ---------------------------------------------------------------------- */
void FixReaxBonds::init()
{
// insure ReaxFF is defined
if (force->pair_match("reax",1) == NULL)
error->all(FLERR,"Cannot use fix reax/bonds without pair_style reax");
}
/* ---------------------------------------------------------------------- */
void FixReaxBonds::end_of_step()
{
OutputReaxBonds(update->ntimestep,fp);
if (me == 0) fflush(fp);
}
/* ---------------------------------------------------------------------- */
void FixReaxBonds::OutputReaxBonds(bigint ntimestep, FILE *fp)
{
int nparticles,nparticles_tot,nbuf,nbuf_local,most,j;
int ii,jn,mbond,numbonds,nsbmax,nsbmax_most;
int nprocs,nlocal_tmp,itmp;
int k,kk,jj,jbufknum;
double cutof3;
double *buf;
MPI_Request irequest;
MPI_Comm_size(world,&nprocs);
nparticles = atom->nlocal;
nparticles_tot = static_cast<int> (atom->natoms);
jn = ReaxParams::nat;
mbond = ReaxParams::mbond;
FORTRAN(getnsbmax,GETNSBMAX)(&nsbmax);
FORTRAN(getcutof3,GETCUTOF3)(&cutof3);
MPI_Allreduce(&nparticles,&most,1,MPI_INT,MPI_MAX,world);
MPI_Allreduce(&nsbmax,&nsbmax_most,1,MPI_INT,MPI_MAX,world);
if (me == 0) {
fprintf(fp,"# Timestep " BIGINT_FORMAT " \n",ntimestep);
fprintf(fp,"# \n");
fprintf(fp,"# Number of particles %d \n",nparticles_tot);
fprintf(fp,"# \n");
fprintf(fp,"# Max number of bonds per atom %d with "
"coarse bond order cutoff %5.3f \n",
nsbmax_most,cutof3);
fprintf(fp,"# Particle connection table and bond orders \n");
fprintf(fp,"# id type nb id_1...id_nb mol bo_1...bo_nb abo nlp q \n");
}
// allocate a temporary buffer for the snapshot info
// big enough for largest number of atoms on any one proc
// nbuf_local = size of local buffer for table of atom bonds
nbuf = 1+(2*nsbmax_most+7)*most;
memory->create(buf,nbuf,"reax/bonds:buf");
j = 0;
buf[j++] = nparticles;
for (int iparticle=0;iparticle<nparticles;iparticle++) {
buf[j++] = atom->tag[iparticle]; //atom tag
buf[j++] = FORTRAN(cbkia,CBKIA).iag[iparticle]; //atom type
jbufknum = j++;
numbonds = FORTRAN(cbkia,CBKIA).iag[iparticle+jn];
// connection table based on coarse bond order cutoff (> cutof3)
kk = 0;
for (k=0;k<numbonds;k++) {
ii = FORTRAN(cbknubon2,CBKNUBON2).nubon1[iparticle+jn*k];
if (FORTRAN(cbkbo,CBKBO).bo[ii-1] > cutof3) {
kk++;
jj = FORTRAN(cbkia,CBKIA).iag[iparticle+jn*(k+2)];
buf[j++] = FORTRAN(cbkc,CBKC).itag[jj-1];
}
}
buf[jbufknum] = kk; //no.bonds
buf[j++]=FORTRAN(cbkia,CBKIA).iag[iparticle+jn*(mbond+2)]; //molec.id
// bond orders (> cutof3)
kk = 0;
for (k=0;k<numbonds;k++) {
ii = FORTRAN(cbknubon2,CBKNUBON2).nubon1[iparticle+jn*k];
if (FORTRAN(cbkbo,CBKBO).bo[ii-1] > cutof3) {
kk++;
buf[j++] = FORTRAN(cbkbo,CBKBO).bo[ii-1];
}
}
// atom bond order (abo), no. of lone pairs (vlp), charge (ch)
buf[j++] = FORTRAN(cbkabo,CBKABO).abo[iparticle];
buf[j++] = FORTRAN(cbklonpar,CBKLONPAR).vlp[iparticle];
buf[j++] = atom->q[iparticle];
}
nbuf_local = j-1;
// node 0 pings each node, receives their buffer, writes to file
// all other nodes wait for ping, send buffer to node 0
if (me == 0) {
for (int inode = 0; inode<nprocs; inode++) {
j = 0;
if (inode == 0) {
nlocal_tmp = nparticles;
j++;
} else {
MPI_Irecv(&buf[0],nbuf,MPI_DOUBLE,inode,0,world,&irequest);
MPI_Send(&itmp,0,MPI_INT,inode,0,world);
MPI_Wait(&irequest,MPI_STATUS_IGNORE);
nlocal_tmp = nint(buf[j++]);
}
for (int iparticle=0;iparticle<nlocal_tmp;iparticle++) {
// print atom tag, atom type, no.bonds
numbonds = nint(buf[j+2]);
fprintf(fp," %d %d %d",nint(buf[j]),nint(buf[j+1]),numbonds);
j += 3;
if (numbonds > nsbmax_most)
error->one(FLERR,"Fix reax/bonds numbonds > nsbmax_most");
// print connection table
for (k=0;k<numbonds;k++)
fprintf(fp," %d",nint(buf[j++]));
// print molecule id
fprintf(fp," %d",nint(buf[j++]));
// print bond orders
for (k=0;k<numbonds;k++)
fprintf(fp,"%14.3f",buf[j++]);
// print sum of bond orders, no. of lone pairs, charge
fprintf(fp,"%14.3f%14.3f%14.3f\n",buf[j],buf[j+1],buf[j+2]);
j+=3;
}
}
} else {
MPI_Recv(&itmp,0,MPI_INT,0,0,world,MPI_STATUS_IGNORE);
MPI_Rsend(&buf[0],nbuf_local,MPI_DOUBLE,0,0,world);
}
if (me == 0) fprintf(fp,"# \n");
memory->destroy(buf);
}
/* ---------------------------------------------------------------------- */
int FixReaxBonds::nint(const double &r)
{
int i = 0;
if (r>0.0) i = static_cast<int>(r+0.5);
else if (r<0.0) i = static_cast<int>(r-0.5);
return i;
}

77
src/REAX/fix_reax_bonds.h
View File

@ -1,77 +0,0 @@
/* -*- c++ -*- ----------------------------------------------------------
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.
------------------------------------------------------------------------- */
#ifdef FIX_CLASS
FixStyle(reax/bonds,FixReaxBonds)
#else
#ifndef LMP_FIX_REAX_BONDS_H
#define LMP_FIX_REAX_BONDS_H
#include <cstdio>
#include "fix.h"
namespace LAMMPS_NS {
class FixReaxBonds : public Fix {
public:
FixReaxBonds(class LAMMPS *, int, char **);
~FixReaxBonds();
int setmask();
void init();
void setup(int);
void end_of_step();
private:
int me;
int nfreq;
FILE *fp;
void OutputReaxBonds(bigint, FILE*);
int nint(const double&);
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Cannot open gzipped file
UNDOCUMENTED
E: Cannot open fix reax/bonds file %s
The output file for the fix reax/bonds command cannot be opened.
Check that the path and name are correct.
E: Cannot use fix reax/bonds without pair_style reax
Self-explanatory.
E: Fix reax/bonds numbonds > nsbmax_most
The limit of the number of bonds expected by the ReaxFF force field
was exceeded.
*/

1067
src/REAX/pair_reax.cpp
View File

File diff suppressed because it is too large Load Diff

154
src/REAX/pair_reax.h
View File

@ -1,154 +0,0 @@
/* -*- c++ -*- ----------------------------------------------------------
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.
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
PairStyle(reax,PairREAX)
#else
#ifndef LMP_PAIR_REAX_H
#define LMP_PAIR_REAX_H
#include "pair.h"
namespace LAMMPS_NS {
class PairREAX : public Pair {
public:
PairREAX(class LAMMPS *);
~PairREAX();
void compute(int, int);
void settings(int, char **);
void coeff(int, char **);
void init_style();
double init_one(int, int);
double memory_usage();
int pack_forward_comm(int, int *, double *, int, int *);
void unpack_forward_comm(int, int, double *);
int pack_reverse_comm(int, int, double *);
void unpack_reverse_comm(int, int *, double *);
private:
double cutmax;
double rcutvsq,rcutbsq;
int iprune,ihb,ihbnew,itripstaball;
double hbcut,swb;
double swa;
double swc0, swc1, swc2, swc3, swc4, swc5, swc6, swc7;
double precision;
int packflag;
struct ff_params {
double rcutsq;
int np;
double *params;
};
ff_params *param_list;
int *map;
int nentries;
double chpot;
int *arow_ptr,*acol_ind;
double *ch,*elcvec;
double *rcg,*wcg,*pcg,*poldcg,*qcg;
double *aval;
int nmax,matmax;
void allocate();
void read_files(char *, char *);
void neigh_f2c(int, int *, int *, int **);
void neigh_c2f(int, int *, int *, int **);
void write_reax_positions();
void write_reax_vlist();
void read_reax_forces();
void read_reax_atom_virial();
void taper_setup();
double taper_E(const double &, const double &);
double taper_F(const double &, const double &);
void compute_charge(double &);
void sparse_product(const int &, const int &, const int &, double[],
int[], int[], double[], double[]);
void cg_solve(const int &, const int &, double[], int[],
int[], double[], double[]);
void charge_reax(const int &, const int &, double[],
double[], int[], int[], double[]);
};
}
#endif
#endif
/* ERROR/WARNING messages:
W: The pair_style reax command is unsupported. Please switch to pair_style reax/c instead
UNDOCUMENTED
E: Reax_defs.h setting for NATDEF is too small
Edit the setting in the ReaxFF library and re-compile the
library and re-build LAMMPS.
E: Reax_defs.h setting for NNEIGHMAXDEF is too small
Edit the setting in the ReaxFF library and re-compile the
library and re-build LAMMPS.
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
E: Cannot currently use pair reax with pair hybrid
This is not yet supported.
E: Pair style reax requires atom IDs
This is a requirement to use the ReaxFF potential.
E: Pair style reax requires newton pair on
This is a requirement to use the ReaxFF potential.
E: Pair style reax requires atom attribute q
The atom style defined does not have this attribute.
W: Not using real units with pair reax
This is most likely an error, unless you have created your own ReaxFF
parameter file in a different set of units.
E: Invalid REAX atom type
There is a mis-match between LAMMPS atom types and the elements
listed in the ReaxFF force field file.
U: pair style reax is now deprecated and will soon be retired. Users should switch to pair_style reax/c
Self-explanatory.
*/

210
src/REAX/pair_reax_fortran.h
View File

@ -1,210 +0,0 @@
/* -*- c++ -*- ----------------------------------------------------------
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.
------------------------------------------------------------------------- */
// machine-specific C++ -> Fortran calling syntax
// It defines the FORTRAN macro for converting variable and function
// names from FORTRAN to C. Different compilers do this in different
// ways. The default is add an underscore to the lower case string.
// Other definitions of the macro can be invoked by defining the
// corresponding macro at compile time using -D e.g. -D_IBM
// CONS(a,b) should return ab, the concatenation of its arguments.
// If compiler is using strict ISO C standards, the ## works.
// Otherwise try the old /**/ trick and test.
// If that fails, you will need to figure out
// a definition for the FORTRAN macro that works on your machine.
#if __STDC__
#define CONS(a,b) a##b
#elif defined(_IBM)
#define CONS(a,b) a##b
#else
#define CONS(a,b) a/**/b
#warning "The following declaration is a test of the CONS macro"
#warning "If it fails, pair_reax_fortran.h must be modified by hand"
static int my_apples_my_oranges = 1;
static int my_applesoroanges = CONS(my_apples,_my_oranges);
#endif
#ifdef _IBM
#define FORTRAN(lcname,ucname) lcname
#endif
#ifdef _F2C_LINUX
#define FORTRAN(lcname,ucname) CONS(lcname,__)
#endif
#ifndef FORTRAN
#define FORTRAN(lcname,ucname) CONS(lcname,_)
#endif
// hard-wired array sizes set in Fortran library
// accesses include file from Fortran library
#include "reax_defs.h"
class ReaxParams {
public:
enum {nneighmax=NNEIGHMAXDEF,
nat=NATDEF,
nattot=NATTOTDEF,
nsort=NSORTDEF,
mbond=MBONDDEF,
nbomax=NBOMAXDEF,
};
};
// data structures corresponding to values in Fortran library
extern "C" struct {
double abo[ReaxParams::nat];
} FORTRAN(cbkabo,CBKABO);
extern "C" struct {
double bo[ReaxParams::nbomax];
} FORTRAN(cbkbo,CBKBO);
extern "C" struct {
double c[3*ReaxParams::nat]; double cglobal[3*ReaxParams::nattot];
int itag[ReaxParams::nat];
} FORTRAN(cbkc,CBKC);
extern "C" struct {
double ch[ReaxParams::nat];
} FORTRAN(cbkch,CBKCH);
extern "C" struct {
double chi[ReaxParams::nsort];
double eta[ReaxParams::nsort];
double gam[ReaxParams::nsort];
} FORTRAN(cbkchb,CBKCHB);
extern "C" struct {
double d[3*ReaxParams::nat]; double estrain[ReaxParams::nat];
} FORTRAN(cbkd,CBKD);
extern "C" struct {
double atomvirial[6*ReaxParams::nat];
double virial[6];
int Lvirial;
int Latomvirial;
} FORTRAN(cbkvirial,CBKVIRIAL);
extern "C" struct {
int ia[ReaxParams::nat*(ReaxParams::mbond+3)];
int iag[ReaxParams::nat*(ReaxParams::mbond+3)];
} FORTRAN(cbkia,CBKIA);
extern "C" struct {
double vlp[ReaxParams::nat];
double dvlpdsbo[ReaxParams::nat];
} FORTRAN(cbklonpar,CBKLONPAR);
extern "C" struct {
int nubon1[ReaxParams::nat*(ReaxParams::mbond)];
int nubon2[ReaxParams::nat*(ReaxParams::mbond)];
} FORTRAN(cbknubon2,CBKNUBON2);
extern "C" struct {
int nvl1[ReaxParams::nneighmax * ReaxParams::nat];
int nvl2[ReaxParams::nneighmax * ReaxParams::nat];
int nvpair;
int nvlself;
} FORTRAN(cbkpairs,CBKPAIRS);
extern "C" struct {
int nvlbo[ReaxParams::nneighmax * ReaxParams::nat];
} FORTRAN(cbknvlbo,CBKNVLBO);
extern "C" struct {
int nvlown[ReaxParams::nneighmax * ReaxParams::nat];
} FORTRAN(cbknvlown,CBKNVLOWN);
extern "C" struct {
char qa[20*ReaxParams::nattot+10];
} FORTRAN(cbkqa,CBKQA);
extern "C" struct {
double eb;
double eoop;
double epen;
double estrc;
double deda[3];
double pressu;
double efi;
double elp;
double emol;
double ea;
double eres;
double et;
double eradbo;
double ev;
double eco;
double ecoa;
double ehb;
double sw;
double ew;
double ep;
double ekin;
} FORTRAN(cbkenergies,CBKENERGIES);
extern "C" struct {
double tset;
double dseed;
double tempmd;
double ts2;
double ts22;
int nmolo;
int nmolo5;
int nbon;
int na;
int namov;
int na_local;
} FORTRAN(rsmall,RSMALL);
// external routines provided by Fortran library
extern "C" void FORTRAN(readc,READC)();
extern "C" void FORTRAN(reaxinit,REAXINIT)();
extern "C" void FORTRAN(ffinpt,FFINPT)();
extern "C" void FORTRAN(tap7th,TAP7TH)();
extern "C" void FORTRAN(taper,TAPER)(double*,double*);
extern "C" void FORTRAN(readgeo,READGEO)();
extern "C" void FORTRAN(srtatom,SRTATOM)();
extern "C" void FORTRAN(vlist,VLIST) ();
extern "C" void FORTRAN(srtbon1,SRTBON1)(int*,int*,double*,int*,int*);
extern "C" void FORTRAN(molec,MOLEC)();
extern "C" void FORTRAN(encalc,ENCALC)();
extern "C" void FORTRAN(getswb,GETSWB)(double*);
extern "C" void FORTRAN(getswa,GETSWA)(double*);
extern "C" void FORTRAN(getvrange,GET_VRANGE)(double*);
extern "C" void FORTRAN(getnvlist,GET_NVLIST)(int*);
extern "C" void FORTRAN(getvlbora,GETVLBORA)(double*);
extern "C" void FORTRAN(cgsolve,CGSOLVE)
(int*,double*,int*,double*,double*,int*);
extern "C" void FORTRAN(getnval,GETNVAL)(int*);
extern "C" void FORTRAN(getntor,GETNTOR)(int*);
extern "C" void FORTRAN(getnhb,GETNHB)(int*);
extern "C" void FORTRAN(getnbonall,GETNBONALL)(int*);
extern "C" void FORTRAN(getnneighmax,GETNNEIGHMAX)(int*);
extern "C" void FORTRAN(getnat,GETNAT)(int*);
extern "C" void FORTRAN(getnattot,GETNATTOT)(int*);
extern "C" void FORTRAN(getnsort,GETNSORT)(int*);
extern "C" void FORTRAN(getmbond,GETMBOND)(int*);
extern "C" void FORTRAN(getnso,GETNSO)(int*);
extern "C" void FORTRAN(setngeofor,SETNGEOFOR)(int*);
extern "C" void FORTRAN(mdsav,MDSAV)(int*);
extern "C" void FORTRAN(getnsbmax,GETNSBMAX)(int*);
extern "C" void FORTRAN(getnsbma2,GETNSBMA2)(int*);
extern "C" void FORTRAN(getcutof3,GETCUTOF3)(double*);