From 36776f425b93993bd4e38b12e1cb3f8512d2c54f Mon Sep 17 00:00:00 2001 From: sjplimp Date: Mon, 21 Mar 2016 16:01:21 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@14756 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/GRANULAR/fix_pour.cpp | 9 +-- src/MANYBODY/pair_airebo.cpp | 101 ++++++++++++++++++++++++--- src/MANYBODY/pair_airebo.h | 6 +- src/MANYBODY/pair_rebo.cpp | 1 - src/USER-OMP/pair_airebo_omp.cpp | 17 +++-- src/USER-REAXC/fix_qeq_reax.cpp | 3 +- src/USER-REAXC/fix_reaxc_species.cpp | 10 ++- src/VORONOI/README | 2 +- src/domain.cpp | 13 ++-- 9 files changed, 132 insertions(+), 30 deletions(-) diff --git a/src/GRANULAR/fix_pour.cpp b/src/GRANULAR/fix_pour.cpp index 662d46667e..65b2223e1d 100644 --- a/src/GRANULAR/fix_pour.cpp +++ b/src/GRANULAR/fix_pour.cpp @@ -42,6 +42,7 @@ enum{ONE,RANGE,POLY}; enum{LAYOUT_UNIFORM,LAYOUT_NONUNIFORM,LAYOUT_TILED}; // several files #define EPSILON 0.001 +#define SMALL 1.0e-10 /* ---------------------------------------------------------------------- */ @@ -900,13 +901,13 @@ void FixPour::options(int narg, char **arg) molfrac[nmol-1] = 1.0; iarg += 2; } else if (strcmp(arg[iarg],"molfrac") == 0) { - if (mode != MOLECULE) error->all(FLERR,"Illegal fix deposit command"); - if (iarg+nmol+1 > narg) error->all(FLERR,"Illegal fix deposit command"); + if (mode != MOLECULE) error->all(FLERR,"Illegal fix pour command"); + if (iarg+nmol+1 > narg) error->all(FLERR,"Illegal fix pour command"); molfrac[0] = force->numeric(FLERR,arg[iarg+1]); for (int i = 1; i < nmol; i++) molfrac[i] = molfrac[i-1] + force->numeric(FLERR,arg[iarg+i+1]); if (molfrac[nmol-1] < 1.0-EPSILON || molfrac[nmol-1] > 1.0+EPSILON) - error->all(FLERR,"Illegal fix deposit command"); + error->all(FLERR,"Illegal fix pour command"); molfrac[nmol-1] = 1.0; iarg += nmol+1; @@ -976,7 +977,7 @@ void FixPour::options(int narg, char **arg) } double sum = 0.0; for (int i = 0; i < npoly; i++) sum += frac_poly[i]; - if (sum != 1.0) + if (fabs(sum - 1.0) > SMALL) error->all(FLERR,"Fix pour polydisperse fractions do not sum to 1.0"); } else error->all(FLERR,"Illegal fix pour command"); diff --git a/src/MANYBODY/pair_airebo.cpp b/src/MANYBODY/pair_airebo.cpp index d38440a1ff..076848ddf0 100644 --- a/src/MANYBODY/pair_airebo.cpp +++ b/src/MANYBODY/pair_airebo.cpp @@ -15,6 +15,8 @@ Contributing author: Ase Henry (MIT) Bugfixes and optimizations: Marcel Fallet & Steve Stuart (Clemson), Axel Kohlmeyer (Temple U) + AIREBO-M modification to optionally replace LJ with Morse potentials. + Thomas C. O'Connor (JHU) 2014 ------------------------------------------------------------------------- */ #include @@ -51,6 +53,8 @@ PairAIREBO::PairAIREBO(LAMMPS *lmp) : Pair(lmp) single_enable = 0; one_coeff = 1; ghostneigh = 1; + ljflag = torflag = 1; + morseflag = 0; maxlocal = 0; REBO_numneigh = NULL; @@ -142,7 +146,6 @@ void PairAIREBO::settings(int narg, char **arg) cutlj = force->numeric(FLERR,arg[0]); - ljflag = torflag = 1; if (narg == 3) { ljflag = force->inumeric(FLERR,arg[1]); torflag = force->inumeric(FLERR,arg[2]); @@ -290,10 +293,23 @@ double PairAIREBO::init_one(int i, int j) cutghost[i][j] = rcmax[ii][jj]; cutljsq[ii][jj] = cutlj*sigma[ii][jj] * cutlj*sigma[ii][jj]; - lj1[ii][jj] = 48.0 * epsilon[ii][jj] * pow(sigma[ii][jj],12.0); - lj2[ii][jj] = 24.0 * epsilon[ii][jj] * pow(sigma[ii][jj],6.0); - lj3[ii][jj] = 4.0 * epsilon[ii][jj] * pow(sigma[ii][jj],12.0); - lj4[ii][jj] = 4.0 * epsilon[ii][jj] * pow(sigma[ii][jj],6.0); + + if (morseflag) { + + // using LJ precomputed parameter arrays to store values for Morse potential + + lj1[ii][jj] = epsilonM[ii][jj] * exp(alphaM[ii][jj]*reqM[ii][jj]); + lj2[ii][jj] = exp(alphaM[ii][jj]*reqM[ii][jj]); + lj3[ii][jj] = 2*epsilonM[ii][jj]*alphaM[ii][jj]*exp(alphaM[ii][jj]*reqM[ii][jj]); + lj4[ii][jj] = alphaM[ii][jj]; + + } else { + + lj1[ii][jj] = 48.0 * epsilon[ii][jj] * pow(sigma[ii][jj],12.0); + lj2[ii][jj] = 24.0 * epsilon[ii][jj] * pow(sigma[ii][jj],6.0); + lj3[ii][jj] = 4.0 * epsilon[ii][jj] * pow(sigma[ii][jj],12.0); + lj4[ii][jj] = 4.0 * epsilon[ii][jj] * pow(sigma[ii][jj],6.0); + } cutghost[j][i] = cutghost[i][j]; cutljsq[jj][ii] = cutljsq[ii][jj]; @@ -732,11 +748,20 @@ void PairAIREBO::FLJ(int eflag, int vflag) dslw = 0.0; } - r2inv = 1.0/rijsq; - r6inv = r2inv*r2inv*r2inv; + if (morseflag) { - vdw = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]); - dvdw = -r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]) / rij; + const double exr = exp(-rij*lj4[itype][jtype]); + vdw = lj1[itype][jtype]*exr*(lj2[itype][jtype]*exr - 2); + dvdw = lj3[itype][jtype]*exr*(1-lj2[itype][jtype]*exr); + + } else { + + r2inv = 1.0/rijsq; + r6inv = r2inv*r2inv*r2inv; + + vdw = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]); + dvdw = -r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]) / rij; + } // VLJ now becomes vdw * slw, derivaties, etc. @@ -3322,6 +3347,10 @@ void PairAIREBO::read_file(char *filename) epsilon_CC,epsilon_CH,epsilon_HH; double sigma_CC,sigma_CH,sigma_HH,epsilonT_CCCC,epsilonT_CCCH,epsilonT_HCCH; + // additional parameters for Morse potential. + double epsilonM_CC,epsilonM_CH,epsilonM_HH,alphaM_CC,alphaM_CH,alphaM_HH; + double reqM_CC,reqM_CH,reqM_HH; + MPI_Comm_rank(world,&me); // read file on proc 0 @@ -3330,7 +3359,10 @@ void PairAIREBO::read_file(char *filename) FILE *fp = force->open_potential(filename); if (fp == NULL) { char str[128]; - sprintf(str,"Cannot open AIREBO potential file %s",filename); + if (morseflag) + sprintf(str,"Cannot open AIREBO-M potential file %s",filename); + else + sprintf(str,"Cannot open AIREBO potential file %s",filename); error->one(FLERR,str); } @@ -3477,6 +3509,29 @@ void PairAIREBO::read_file(char *filename) fgets(s,MAXLINE,fp); sscanf(s,"%lg",&epsilonT_HCCH); + if (morseflag) { + // lines for reading in MORSE parameters from CH.airebo_m file + fgets(s,MAXLINE,fp); + sscanf(s,"%lg",&epsilonM_CC); + fgets(s,MAXLINE,fp); + sscanf(s,"%lg",&epsilonM_CH); + fgets(s,MAXLINE,fp); + sscanf(s,"%lg",&epsilonM_HH); + fgets(s,MAXLINE,fp); + sscanf(s,"%lg",&alphaM_CC); + fgets(s,MAXLINE,fp); + sscanf(s,"%lg",&alphaM_CH); + fgets(s,MAXLINE,fp); + sscanf(s,"%lg",&alphaM_HH); + fgets(s,MAXLINE,fp); + sscanf(s,"%lg",&reqM_CC); + fgets(s,MAXLINE,fp); + sscanf(s,"%lg",&reqM_CH); + fgets(s,MAXLINE,fp); + sscanf(s,"%lg",&reqM_HH); + } + + // gC spline fgets(s,MAXLINE,fp); @@ -3805,6 +3860,25 @@ void PairAIREBO::read_file(char *filename) sigma[1][0] = sigma[0][1]; sigma[1][1] = sigma_HH; + if (morseflag) { + // Morse parameter assignments + + epsilonM[0][0] = epsilonM_CC; + epsilonM[0][1] = epsilonM_CH; + epsilonM[1][0] = epsilonM[0][1]; + epsilonM[1][1] = epsilonM_HH; + + alphaM[0][0] = alphaM_CC; + alphaM[0][1] = alphaM_CH; + alphaM[1][0] = alphaM[0][1]; + alphaM[1][1] = alphaM_HH; + + reqM[0][0] = reqM_CC; + reqM[0][1] = reqM_CH; + reqM[1][0] = reqM[0][1]; + reqM[1][1] = reqM_HH; + } + // torsional thmin = -1.0; @@ -3854,6 +3928,13 @@ void PairAIREBO::read_file(char *filename) MPI_Bcast(&sigma[0][0],4,MPI_DOUBLE,0,world); MPI_Bcast(&epsilonT[0][0],4,MPI_DOUBLE,0,world); + if (morseflag) { + // Morse parameter broadcast + MPI_Bcast(&epsilonM[0][0],4,MPI_DOUBLE,0,world); + MPI_Bcast(&alphaM[0][0],4,MPI_DOUBLE,0,world); + MPI_Bcast(&reqM[0][0],4,MPI_DOUBLE,0,world); + } + MPI_Bcast(&gCdom[0],5,MPI_DOUBLE,0,world); MPI_Bcast(&gC1[0][0],24,MPI_DOUBLE,0,world); MPI_Bcast(&gC2[0][0],24,MPI_DOUBLE,0,world); diff --git a/src/MANYBODY/pair_airebo.h b/src/MANYBODY/pair_airebo.h index b929621bd3..06bf6da554 100644 --- a/src/MANYBODY/pair_airebo.h +++ b/src/MANYBODY/pair_airebo.h @@ -42,7 +42,8 @@ class PairAIREBO : public Pair { int *map; // 0 (C), 1 (H), or -1 (NULL) for each type int me; - int ljflag,torflag; // 0/1 if LJ,torsion terms included + int ljflag,torflag; // 0/1 if LJ/Morse,torsion terms included + int morseflag; // 1 if Morse instead of LJ for non-bonded double cutlj; // user-specified LJ cutoff double cutljrebosq; // cut for when to compute @@ -68,6 +69,9 @@ class PairAIREBO : public Pair { double rcLJmin[2][2],rcLJmax[2][2],rcLJmaxsq[2][2],bLJmin[2][2],bLJmax[2][2]; double epsilon[2][2],sigma[2][2],epsilonT[2][2]; + // parameters for Morse variant + double epsilonM[2][2],alphaM[2][2],reqM[2][2]; + // spline coefficients double gCdom[5],gC1[4][6],gC2[4][6],gHdom[4],gH[3][6]; diff --git a/src/MANYBODY/pair_rebo.cpp b/src/MANYBODY/pair_rebo.cpp index f2ca30df7f..5b28eae4dd 100644 --- a/src/MANYBODY/pair_rebo.cpp +++ b/src/MANYBODY/pair_rebo.cpp @@ -12,7 +12,6 @@ ------------------------------------------------------------------------- */ #include "pair_rebo.h" -#include "force.h" #include "error.h" using namespace LAMMPS_NS; diff --git a/src/USER-OMP/pair_airebo_omp.cpp b/src/USER-OMP/pair_airebo_omp.cpp index 8215e9967d..fc6ba17f1c 100644 --- a/src/USER-OMP/pair_airebo_omp.cpp +++ b/src/USER-OMP/pair_airebo_omp.cpp @@ -2178,11 +2178,20 @@ void PairAIREBOOMP::FLJ_thr(int ifrom, int ito, int evflag, int eflag, dslw = 0.0; } - r2inv = 1.0/rijsq; - r6inv = r2inv*r2inv*r2inv; + if (morseflag) { - vdw = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]); - dvdw = -r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]) / rij; + const double exr = exp(-rij*lj4[itype][jtype]); + vdw = lj1[itype][jtype]*exr*(lj2[itype][jtype]*exr - 2); + dvdw = lj3[itype][jtype]*exr*(1-lj2[itype][jtype]*exr); + + } else { + + r2inv = 1.0/rijsq; + r6inv = r2inv*r2inv*r2inv; + + vdw = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]); + dvdw = -r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]) / rij; + } // VLJ now becomes vdw * slw, derivaties, etc. diff --git a/src/USER-REAXC/fix_qeq_reax.cpp b/src/USER-REAXC/fix_qeq_reax.cpp index 794c23da3a..bb0fe3e6d7 100644 --- a/src/USER-REAXC/fix_qeq_reax.cpp +++ b/src/USER-REAXC/fix_qeq_reax.cpp @@ -367,7 +367,8 @@ void FixQEqReax::init_shielding() int ntypes; ntypes = atom->ntypes; - memory->create(shld,ntypes+1,ntypes+1,"qeq:shileding"); + if (shld == NULL) + memory->create(shld,ntypes+1,ntypes+1,"qeq:shileding"); for( i = 1; i <= ntypes; ++i ) for( j = 1; j <= ntypes; ++j ) diff --git a/src/USER-REAXC/fix_reaxc_species.cpp b/src/USER-REAXC/fix_reaxc_species.cpp index c063869f4d..05e861d554 100644 --- a/src/USER-REAXC/fix_reaxc_species.cpp +++ b/src/USER-REAXC/fix_reaxc_species.cpp @@ -58,6 +58,8 @@ FixReaxCSpecies::FixReaxCSpecies(LAMMPS *lmp, int narg, char **arg) : size_peratom_cols = 0; peratom_freq = 1; + nvalid = -1; + MPI_Comm_rank(world,&me); MPI_Comm_size(world,&nprocs); ntypes = atom->ntypes; @@ -269,7 +271,8 @@ int FixReaxCSpecies::setmask() void FixReaxCSpecies::setup(int vflag) { ntotal = static_cast (atom->natoms); - memory->create(Name,ntypes,"reax/c/species:Name"); + if (Name == NULL) + memory->create(Name,ntypes,"reax/c/species:Name"); post_integrate(); } @@ -286,7 +289,10 @@ void FixReaxCSpecies::init() "pair_style reax/c"); reaxc->fixspecies_flag = 1; - nvalid = update->ntimestep+nfreq; + + // reset next output timestep if not yet set or timestep has been reset + if (nvalid != update->ntimestep) + nvalid = update->ntimestep+nfreq; // check if this fix has been called twice int count = 0; diff --git a/src/VORONOI/README b/src/VORONOI/README index 89e3b84671..9a83c95a8d 100644 --- a/src/VORONOI/README +++ b/src/VORONOI/README @@ -1,7 +1,7 @@ The VORONOI package adds a compute voronoi/atom command which calculates a Voronoi tesselation of the system. -It uses the Voro++ library, available at http://math.lbl.gov/voro++ to +It uses the Voro++ library, available at http://math.lbl.gov/voro++/ to compute the tesselation locally on each processor. Voro++ was developed by Chris H. Rycroft while at UC Berkeley / Lawrence Berkeley Laboratory. diff --git a/src/domain.cpp b/src/domain.cpp index 601d59cf7c..ee6d2c1bc5 100644 --- a/src/domain.cpp +++ b/src/domain.cpp @@ -503,12 +503,13 @@ void Domain::pbc() double *coord; int n3 = 3*nlocal; - if (x) coord = &x[0][0]; - int flag = 0; - for (i = 0; i < n3; i++) - if (!ISFINITE(*coord++)) flag = 1; - if (flag) error->one(FLERR,"Non-numeric atom coords - simulation unstable"); - + if (x) { + coord = &x[0][0]; + int flag = 0; + for (i = 0; i < n3; i++) + if (!ISFINITE(*coord++)) flag = 1; + if (flag) error->one(FLERR,"Non-numeric atom coords - simulation unstable"); + } // setup for PBC checks if (triclinic == 0) {