git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@14756 f3b2605a-c512-4ea7-a41b-209d697bcdaa

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 (mode != MOLECULE) error->all(FLERR,"Illegal fix pour command");
-      if (iarg+nmol+1 > narg) error->all(FLERR,"Illegal fix deposit 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");
 

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 <math.h>
@@ -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);
+  if (morseflag) {
-  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);
+    // 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;
+      if (morseflag) {
-      r6inv = r2inv*r2inv*r2inv;
 
-      vdw = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]);
+        const double exr = exp(-rij*lj4[itype][jtype]);
-      dvdw = -r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]) / rij;
+        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);

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];

src/MANYBODY/pair_rebo.cpp

@@ -12,7 +12,6 @@
 ------------------------------------------------------------------------- */
 
 #include "pair_rebo.h"
-#include "force.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;

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;
+      if (morseflag) {
-      r6inv = r2inv*r2inv*r2inv;
 
-      vdw = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]);
+        const double exr = exp(-rij*lj4[itype][jtype]);
-      dvdw = -r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]) / rij;
+        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.
 

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 )

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<int> (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;

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.

src/domain.cpp

@@ -503,12 +503,13 @@ void Domain::pbc()
 
   double *coord;
   int n3 = 3*nlocal;
-  if (x) coord = &x[0][0];
+  if (x) {
-  int flag = 0;
+    coord = &x[0][0];
-  for (i = 0; i < n3; i++)
+    int flag = 0;
-    if (!ISFINITE(*coord++)) flag = 1;
+    for (i = 0; i < n3; i++)
-  if (flag) error->one(FLERR,"Non-numeric atom coords - simulation unstable");
+      if (!ISFINITE(*coord++)) flag = 1;
-
+    if (flag) error->one(FLERR,"Non-numeric atom coords - simulation unstable");
+  }
   // setup for PBC checks
 
   if (triclinic == 0) {