21Dec16 patch

Merge pull request #304 from timattox/USER-DPD_whitespace
USER-DPD: Whitespace cleanup to pair_dpd_fdt_energy.cpp
2016-12-21 09:53:32 -07:00 · 2016-12-21 09:51:10 -07:00 · 2016-12-21 09:50:29 -07:00 · 2016-12-20 11:17:11 -05:00 · 2016-12-19 15:11:52 -05:00 · 2016-12-16 11:36:54 -07:00
36 changed files with 269 additions and 142 deletions
--- a/doc/src/Manual.txt
+++ b/doc/src/Manual.txt
@ -1,7 +1,7 @@
 <!-- HTML_ONLY -->
 <HEAD>
 <TITLE>LAMMPS Users Manual</TITLE>
-<META NAME="docnumber" CONTENT="13 Dec 2016 version">
+<META NAME="docnumber" CONTENT="21 Dec 2016 version">
 <META NAME="author" CONTENT="http://lammps.sandia.gov - Sandia National Laboratories">
 <META NAME="copyright" CONTENT="Copyright (2003) Sandia Corporation.  This software and manual is distributed under the GNU General Public License.">
 </HEAD>
@ -21,7 +21,7 @@
 <H1></H1>

 LAMMPS Documentation :c,h3
-13 Dec 2016 version :c,h4
+21 Dec 2016 version :c,h4

 Version info: :h4

--- a/doc/src/python.txt
+++ b/doc/src/python.txt
@ -14,7 +14,7 @@ python func keyword args ... :pre

 func = name of Python function :ulb,l
 one or more keyword/args pairs must be appended :l
-keyword = {invoke} or {input} or {return} or {format} or {file} or {here} or {exists}
+keyword = {invoke} or {input} or {return} or {format} or {length} or {file} or {here} or {exists}
  {invoke} arg = none = invoke the previously defined Python function
  {input} args = N i1 i2 ... iN
    N = # of inputs to function
@ -29,6 +29,8 @@ keyword = {invoke} or {input} or {return} or {format} or {file} or {here} or {ex
    M = N+1 if there is a return value
    fstring = each character (i,f,s,p) corresponds in order to an input or return value
    'i' = integer, 'f' = floating point, 's' = string, 'p' = SELF
+  {length} arg = Nlen
+    Nlen = max length of string returned from Python function
  {file} arg = filename
    filename = file of Python code, which defines func
  {here} arg = inline
@ -165,6 +167,17 @@ equal-style variable as an argument, but only if the output of the
 Python function is flagged as a numeric value ("i" or "f") via the
 {format} keyword.

+If the {return} keyword is used and the {format} keyword specifies the
+output as a string, then the default maximum length of that string is
+63 characters (64-1 for the string terminator).  If you want to return
+a longer string, the {length} keyword can be specified with its {Nlen}
+value set to a larger number (the code allocates space for Nlen+1 to
+include the string terminator).  If the Python function generates a
+string longer than the default 63 or the specified {Nlen}, it will be
+trunctated.
+
+:line
+
 Either the {file}, {here}, or {exists} keyword must be used, but only
 one of them.  These keywords specify what Python code to load into the
 Python interpreter.  The {file} keyword gives the name of a file,
--- a/src/PYTHON/python.cpp
+++ b/src/PYTHON/python.cpp
@ -89,6 +89,7 @@ void Python::command(int narg, char **arg)
  istr = NULL;
  ostr = NULL;
  format = NULL;
+  length_longstr = 0;
  char *pyfile = NULL;
  char *herestr = NULL;
  int existflag = 0;
@ -115,6 +116,11 @@ void Python::command(int narg, char **arg)
      format = new char[n];
      strcpy(format,arg[iarg+1]);
      iarg += 2;
+    } else if (strcmp(arg[iarg],"length") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Invalid python command");
+      length_longstr = force->inumeric(FLERR,arg[iarg+1]);
+      if (length_longstr <= 0) error->all(FLERR,"Invalid python command");
+      iarg += 2;
    } else if (strcmp(arg[iarg],"file") == 0) {
      if (iarg+2 > narg) error->all(FLERR,"Invalid python command");
      delete[] pyfile;
@ -249,6 +255,7 @@ void Python::invoke_function(int ifunc, char *result)

  // function returned a value
  // assign it to result string stored by python-style variable
+  // or if user specified a length, assign it to longstr

  if (pfuncs[ifunc].noutput) {
    int otype = pfuncs[ifunc].otype;
@ -258,7 +265,9 @@ void Python::invoke_function(int ifunc, char *result)
      sprintf(result,"%.15g",PyFloat_AsDouble(pValue));
    } else if (otype == STRING) {
      char *pystr = PyString_AsString(pValue);
-      strncpy(result,pystr,VALUELENGTH-1);
+      if (pfuncs[ifunc].longstr) 
+        strncpy(pfuncs[ifunc].longstr,pystr,pfuncs[ifunc].length_longstr);
+      else strncpy(result,pystr,VALUELENGTH-1);
    }
    Py_DECREF(pValue);
  }
@ -287,6 +296,13 @@ int Python::variable_match(char *name, char *varname, int numeric)

 /* ------------------------------------------------------------------ */

+char *Python::long_string(int ifunc)
+{
+  return pfuncs[ifunc].longstr;
+}
+
+/* ------------------------------------------------------------------ */
+
 int Python::create_entry(char *name)
 {
  // ifunc = index to entry by name in pfuncs vector, can be old or new
@ -370,6 +386,7 @@ int Python::create_entry(char *name)
  // process output as value or variable

  pfuncs[ifunc].ovarname = NULL;
+  pfuncs[ifunc].longstr = NULL;
  if (!noutput) return ifunc;

  char type = format[ninput];
@ -378,6 +395,14 @@ int Python::create_entry(char *name)
  else if (type == 's') pfuncs[ifunc].otype = STRING;
  else error->all(FLERR,"Invalid python command");

+  if (length_longstr) {
+    if (pfuncs[ifunc].otype != STRING) 
+      error->all(FLERR,"Python command length keyword "
+                 "cannot be used unless output is a string");
+    pfuncs[ifunc].length_longstr = length_longstr;
+    pfuncs[ifunc].longstr = new char[length_longstr+1];
+  }
+
  if (strstr(ostr,"v_") != ostr) error->all(FLERR,"Invalid python command");
  int n = strlen(&ostr[2]) + 1;
  pfuncs[ifunc].ovarname = new char[n];
@ -398,4 +423,5 @@ void Python::deallocate(int i)
    delete [] pfuncs[i].svalue[j];
  delete [] pfuncs[i].svalue;
  delete [] pfuncs[i].ovarname;
+  delete [] pfuncs[i].longstr;
 }
--- a/src/PYTHON/python.h
+++ b/src/PYTHON/python.h
@ -28,9 +28,10 @@ class Python : protected Pointers {
  void invoke_function(int, char *);
  int find(char *);
  int variable_match(char *, char *, int);
+  char *long_string(int);

 private:
-  int ninput,noutput;
+  int ninput,noutput,length_longstr;
  char **istr;
  char *ostr,*format;
  void *pyMain;
@ -44,6 +45,8 @@ class Python : protected Pointers {
    char **svalue;
    int otype;
    char *ovarname;
+    char *longstr;
+    int length_longstr;
    void *pFunc;
  };

--- a/src/USER-DPD/compute_dpd_atom.cpp
+++ b/src/USER-DPD/compute_dpd_atom.cpp
@ -43,7 +43,7 @@ ComputeDpdAtom::ComputeDpdAtom(LAMMPS *lmp, int narg, char **arg) :
  size_peratom_cols = 4;

  nmax = 0;
-  
+
  if (atom->dpd_flag != 1) error->all(FLERR,"compute dpd requires atom_style with internal temperature and energies (e.g. dpd)");
 }

--- a/src/USER-DPD/fix_eos_cv.cpp
+++ b/src/USER-DPD/fix_eos_cv.cpp
@ -70,8 +70,8 @@ void FixEOScv::init()
      if (mask[i] & groupbit) {
        if(dpdTheta[i] <= 0.0)
          error->one(FLERR,"Internal temperature <= zero");
-        uCond[i] = 0.5*cvEOS*dpdTheta[i];
-        uMech[i] = 0.5*cvEOS*dpdTheta[i];
+        uCond[i] = 0.0;
+        uMech[i] = cvEOS*dpdTheta[i];
      }
  }
 }
--- a/src/USER-DPD/fix_eos_cv.h
+++ b/src/USER-DPD/fix_eos_cv.h
@ -50,6 +50,10 @@ 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:  FixEOScv requires atom_style with internal temperature and energies (e.g. dpd)
+
+Self-explanatory.
+
 E: EOS cv must be > 0.0

 The constant volume heat capacity must be larger than zero.
--- a/src/USER-DPD/fix_eos_table.cpp
+++ b/src/USER-DPD/fix_eos_table.cpp
@ -122,8 +122,8 @@ void FixEOStable::init()
        if(dpdTheta[i] <= 0.0)
          error->one(FLERR,"Internal temperature <= zero");
        energy_lookup(dpdTheta[i],tmp);
-        uCond[i] = tmp / 2.0;
-        uMech[i] = tmp / 2.0;
+        uCond[i] = 0.0;
+        uMech[i] = tmp;
      }
  }
 }
--- a/src/USER-DPD/fix_eos_table.h
+++ b/src/USER-DPD/fix_eos_table.h
@ -92,6 +92,10 @@ E: eos/table values are not increasing

 The EOS must be a monotonically, increasing function

+E:  FixEOStable requires atom_style with internal temperature and energies (e.g. dpd)
+
+Self-explanatory.
+
 E: Internal temperature < zero

 Self-explanatory.  EOS may not be valid under current simulation conditions.
--- a/src/USER-DPD/fix_eos_table_rx.cpp
+++ b/src/USER-DPD/fix_eos_table_rx.cpp
@ -34,7 +34,7 @@ using namespace FixConst;
 /* ---------------------------------------------------------------------- */

 FixEOStableRX::FixEOStableRX(LAMMPS *lmp, int narg, char **arg) :
-  Fix(lmp, narg, arg), ntables(0), tables(NULL), 
+  Fix(lmp, narg, arg), ntables(0), tables(NULL),
  tables2(NULL), dHf(NULL), eosSpecies(NULL)
 {
  if (narg != 8) error->all(FLERR,"Illegal fix eos/table/rx command");
@ -162,13 +162,15 @@ void FixEOStableRX::setup(int vflag)
  double *uCG   = atom->uCG;
  double *uCGnew = atom->uCGnew;

-  for (int i = 0; i < nlocal; i++)
-    if (mask[i] & groupbit){
-      duChem = uCG[i] - uCGnew[i];
-      uChem[i] += duChem;
-      uCG[i] = 0.0;
-      uCGnew[i] = 0.0;
-    }
+  if(!this->restart_reset){
+    for (int i = 0; i < nlocal; i++)
+      if (mask[i] & groupbit){
+        duChem = uCG[i] - uCGnew[i];
+        uChem[i] += duChem;
+        uCG[i] = 0.0;
+        uCGnew[i] = 0.0;
+      }
+  }

  // Communicate the updated momenta and velocities to all nodes
  comm->forward_comm_fix(this);
@ -200,8 +202,8 @@ void FixEOStableRX::init()
        if(dpdTheta[i] <= 0.0)
          error->one(FLERR,"Internal temperature <= zero");
        energy_lookup(i,dpdTheta[i],tmp);
-        uCond[i] = tmp / 2.0;
-        uMech[i] = tmp / 2.0;
+        uCond[i] = 0.0;
+        uMech[i] = tmp;
        uChem[i] = 0.0;
      }
  }
--- a/src/USER-DPD/fix_eos_table_rx.h
+++ b/src/USER-DPD/fix_eos_table_rx.h
@ -106,6 +106,10 @@ E:  eos/table/rx values are not increasing

 The equation-of-state must an increasing function

+E:  FixEOStableRX requires atom_style with internal temperature and energies (e.g. dpd)
+
+Self-explanatory.
+
 E:  Internal temperature <= zero.

 Self-explanatory.
--- a/src/USER-DPD/fix_rx.cpp
+++ b/src/USER-DPD/fix_rx.cpp
@ -60,16 +60,15 @@ double getElapsedTime( const TimerType &t0, const TimerType &t1) { return t1-t0;
 /* ---------------------------------------------------------------------- */

 FixRX::FixRX(LAMMPS *lmp, int narg, char **arg) :
-  Fix(lmp, narg, arg), mol2param(NULL), nreactions(0), 
-  params(NULL), Arr(NULL), nArr(NULL), Ea(NULL), tempExp(NULL), 
-  stoich(NULL), stoichReactants(NULL), stoichProducts(NULL), kR(NULL), 
-  pairDPDE(NULL), dpdThetaLocal(NULL), sumWeights(NULL), sparseKinetics_nu(NULL), 
-  sparseKinetics_nuk(NULL), sparseKinetics_inu(NULL), sparseKinetics_isIntegralReaction(NULL), 
-  kineticsFile(NULL), id_fix_species(NULL), 
+  Fix(lmp, narg, arg), mol2param(NULL), nreactions(0),
+  params(NULL), Arr(NULL), nArr(NULL), Ea(NULL), tempExp(NULL),
+  stoich(NULL), stoichReactants(NULL), stoichProducts(NULL), kR(NULL),
+  pairDPDE(NULL), dpdThetaLocal(NULL), sumWeights(NULL), sparseKinetics_nu(NULL),
+  sparseKinetics_nuk(NULL), sparseKinetics_inu(NULL), sparseKinetics_isIntegralReaction(NULL),
+  kineticsFile(NULL), id_fix_species(NULL),
  id_fix_species_old(NULL), fix_species(NULL), fix_species_old(NULL)
 {
  if (narg < 7 || narg > 12) error->all(FLERR,"Illegal fix rx command");
-  restart_peratom = 1;
  nevery = 1;

  nreactions = maxparam = 0;
@ -366,6 +365,7 @@ void FixRX::post_constructor()

  modify->add_fix(nspecies+5,newarg);
  fix_species = (FixPropertyAtom *) modify->fix[modify->nfix-1];
+  restartFlag = modify->fix[modify->nfix-1]->restart_reset;

  modify->add_fix(nspecies+5,newarg2);
  fix_species_old = (FixPropertyAtom *) modify->fix[modify->nfix-1];
@ -647,34 +647,38 @@ void FixRX::setup_pre_force(int vflag)
  double tmp;
  int ii;

-  if(localTempFlag){
-    int count = nlocal + (newton_pair ? nghost : 0);
-    dpdThetaLocal = new double[count];
-    memset(dpdThetaLocal, 0, sizeof(double)*count);
-    computeLocalTemperature();
+  if(restartFlag){
+    restartFlag = 0;
+  } else {
+    if(localTempFlag){
+      int count = nlocal + (newton_pair ? nghost : 0);
+      dpdThetaLocal = new double[count];
+      memset(dpdThetaLocal, 0, sizeof(double)*count);
+      computeLocalTemperature();
+    }
+
+    for (int id = 0; id < nlocal; id++)
+      for (int ispecies=0; ispecies<nspecies; ispecies++){
+        tmp = atom->dvector[ispecies][id];
+        atom->dvector[ispecies+nspecies][id] = tmp;
+      }
+    for (int i = 0; i < nlocal; i++)
+      if (mask[i] & groupbit){
+
+        // Set the reaction rate constants to zero:  no reactions occur at step 0
+        for(int irxn=0;irxn<nreactions;irxn++)
+          kR[irxn] = 0.0;
+
+        if (odeIntegrationFlag == ODE_LAMMPS_RK4)
+          rk4(i,NULL);
+        else if (odeIntegrationFlag == ODE_LAMMPS_RKF45)
+          rkf45(i,NULL);
+      }
+
+    // Communicate the updated momenta and velocities to all nodes
+    comm->forward_comm_fix(this);
+    if(localTempFlag) delete [] dpdThetaLocal;
  }
-
-  for (int id = 0; id < nlocal; id++)
-    for (int ispecies=0; ispecies<nspecies; ispecies++){
-      tmp = atom->dvector[ispecies][id];
-      atom->dvector[ispecies+nspecies][id] = tmp;
-    }
-  for (int i = 0; i < nlocal; i++)
-    if (mask[i] & groupbit){
-
-      // Set the reaction rate constants to zero:  no reactions occur at step 0
-      for(int irxn=0;irxn<nreactions;irxn++)
-        kR[irxn] = 0.0;
-
-      if (odeIntegrationFlag == ODE_LAMMPS_RK4)
-        rk4(i,NULL);
-      else if (odeIntegrationFlag == ODE_LAMMPS_RKF45)
-        rkf45(i,NULL);
-    }
-
-  // Communicate the updated momenta and velocities to all nodes
-  comm->forward_comm_fix(this);
-  if(localTempFlag) delete [] dpdThetaLocal;
 }

 /* ---------------------------------------------------------------------- */
--- a/src/USER-DPD/fix_rx.h
+++ b/src/USER-DPD/fix_rx.h
@ -132,6 +132,7 @@ class FixRX : public Fix {
  char *kineticsFile;
  char *id_fix_species,*id_fix_species_old;
  class FixPropertyAtom *fix_species,*fix_species_old;
+  int restartFlag;

 };

--- a/src/USER-DPD/npair_half_bin_newton_ssa.cpp
+++ b/src/USER-DPD/npair_half_bin_newton_ssa.cpp
@ -12,7 +12,7 @@
 ------------------------------------------------------------------------- */

 /* ----------------------------------------------------------------------
-   Contributing authors: 
+   Contributing authors:
   James Larentzos and Timothy I. Mattox (Engility Corporation)
 ------------------------------------------------------------------------- */

@ -147,7 +147,7 @@ void NPairHalfBinNewtonSSA::build(NeighList *list)
    // loop over all local atoms in other bins in "half" stencil

    for (k = 0; k < nstencil_half; k++) {
-      for (j = binhead_ssa[ibin+stencil[k]]; j >= 0; 
+      for (j = binhead_ssa[ibin+stencil[k]]; j >= 0;
           j = bins_ssa[j]) {

        jtype = type[j];
@ -183,7 +183,7 @@ void NPairHalfBinNewtonSSA::build(NeighList *list)
    // Note2: only non-pure locals can have ghosts as neighbors

    if (ssaAIR[i] == 1) for (k = 0; k < nstencil_full; k++) {
-      for (j = gbinhead_ssa[ibin+stencil[k]]; j >= 0; 
+      for (j = gbinhead_ssa[ibin+stencil[k]]; j >= 0;
           j = bins_ssa[j]) {

        jtype = type[j];
--- a/src/USER-DPD/npair_halffull_newton_ssa.cpp
+++ b/src/USER-DPD/npair_halffull_newton_ssa.cpp
@ -12,7 +12,7 @@
 ------------------------------------------------------------------------- */

 /* ----------------------------------------------------------------------
-   Contributing authors: 
+   Contributing authors:
   James Larentzos and Timothy I. Mattox (Engility Corporation)
 ------------------------------------------------------------------------- */

--- a/src/USER-DPD/npair_halffull_newton_ssa.h
+++ b/src/USER-DPD/npair_halffull_newton_ssa.h
@ -15,7 +15,7 @@

 NPairStyle(halffull/newton/ssa,
           NPairHalffullNewtonSSA,
-           NP_HALFFULL | NP_NSQ | NP_BIN | NP_MULTI | NP_NEWTON | 
+           NP_HALFFULL | NP_NSQ | NP_BIN | NP_MULTI | NP_NEWTON |
           NP_ORTHO | NP_TRI | NP_SSA)

 #else
--- a/src/USER-DPD/nstencil_half_bin_2d_newton_ssa.cpp
+++ b/src/USER-DPD/nstencil_half_bin_2d_newton_ssa.cpp
@ -12,7 +12,7 @@
 ------------------------------------------------------------------------- */

 /* ----------------------------------------------------------------------
-   Contributing authors: 
+   Contributing authors:
   James Larentzos and Timothy I. Mattox (Engility Corporation)
 ------------------------------------------------------------------------- */

@ -24,7 +24,7 @@ using namespace LAMMPS_NS;

 /* ---------------------------------------------------------------------- */

-NStencilHalfBin2dNewtonSSA::NStencilHalfBin2dNewtonSSA(LAMMPS *lmp) : 
+NStencilHalfBin2dNewtonSSA::NStencilHalfBin2dNewtonSSA(LAMMPS *lmp) :
  NStencilSSA(lmp) {}

 /* ----------------------------------------------------------------------
--- a/src/USER-DPD/nstencil_half_bin_3d_newton_ssa.cpp
+++ b/src/USER-DPD/nstencil_half_bin_3d_newton_ssa.cpp
@ -12,7 +12,7 @@
 ------------------------------------------------------------------------- */

 /* ----------------------------------------------------------------------
-   Contributing authors: 
+   Contributing authors:
   James Larentzos and Timothy I. Mattox (Engility Corporation)
 ------------------------------------------------------------------------- */

@ -24,7 +24,7 @@ using namespace LAMMPS_NS;

 /* ---------------------------------------------------------------------- */

-NStencilHalfBin3dNewtonSSA::NStencilHalfBin3dNewtonSSA(LAMMPS *lmp) : 
+NStencilHalfBin3dNewtonSSA::NStencilHalfBin3dNewtonSSA(LAMMPS *lmp) :
  NStencilSSA(lmp) {}

 /* ----------------------------------------------------------------------
@ -62,7 +62,7 @@ void NStencilHalfBin3dNewtonSSA::create()

  // For k==0, make sure to skip already included bins

-  k = 0; 
+  k = 0;
  for (j = -sy; j <= 0; j++)
    for (i = -sx; i <= sx; i++) {
      if (j == 0 && i > 0) continue;
--- a/src/USER-DPD/pair_dpd_fdt.cpp
+++ b/src/USER-DPD/pair_dpd_fdt.cpp
@ -182,7 +182,7 @@ void PairDPDfdt::compute(int eflag, int vflag)
          wr = 1.0 - r/cut[itype][jtype];
          wd = wr*wr;
          randnum = random->gaussian();
-	  gamma_ij = sigma[itype][jtype]*sigma[itype][jtype]
+          gamma_ij = sigma[itype][jtype]*sigma[itype][jtype]
                     / (2.0*force->boltz*temperature);

          // conservative force = a0 * wd
--- a/src/USER-DPD/pair_dpd_fdt_energy.cpp
+++ b/src/USER-DPD/pair_dpd_fdt_energy.cpp
@ -206,7 +206,7 @@ void PairDPDfdtEnergy::compute(int eflag, int vflag)
          if (r < EPSILON) continue;     // r can be 0.0 in DPD systems
          rinv = 1.0/r;
          wr = 1.0 - r/cut[itype][jtype];
-	  wd = wr*wr;
+          wd = wr*wr;

          delvx = vxtmp - v[j][0];
          delvy = vytmp - v[j][1];
@ -214,11 +214,11 @@ void PairDPDfdtEnergy::compute(int eflag, int vflag)
          dot = delx*delvx + dely*delvy + delz*delvz;
          randnum = random->gaussian();

-	  // Compute the current temperature
-	  theta_ij = 0.5*(1.0/dpdTheta[i] + 1.0/dpdTheta[j]);
-	  theta_ij = 1.0/theta_ij;
-	
-	  gamma_ij = sigma[itype][jtype]*sigma[itype][jtype]
+          // Compute the current temperature
+          theta_ij = 0.5*(1.0/dpdTheta[i] + 1.0/dpdTheta[j]);
+          theta_ij = 1.0/theta_ij;
+
+          gamma_ij = sigma[itype][jtype]*sigma[itype][jtype]
                     / (2.0*force->boltz*theta_ij);

          // conservative force = a0 * wr
@ -239,44 +239,44 @@ void PairDPDfdtEnergy::compute(int eflag, int vflag)
            f[j][2] -= delz*fpair;
          }

-	  if (rmass) {
-	    mass_i = rmass[i];
-	    mass_j = rmass[j];
-	  } else {
-	    mass_i = mass[itype];
-	    mass_j = mass[jtype];
-	  }
-	  massinv_i = 1.0 / mass_i;
-	  massinv_j = 1.0 / mass_j;
+          if (rmass) {
+            mass_i = rmass[i];
+            mass_j = rmass[j];
+          } else {
+            mass_i = mass[itype];
+            mass_j = mass[jtype];
+          }
+          massinv_i = 1.0 / mass_i;
+          massinv_j = 1.0 / mass_j;

-	  // Compute the mechanical and conductive energy, uMech and uCond
-	  mu_ij = massinv_i + massinv_j;
-	  mu_ij *= force->ftm2v;
+          // Compute the mechanical and conductive energy, uMech and uCond
+          mu_ij = massinv_i + massinv_j;
+          mu_ij *= force->ftm2v;

-	  uTmp = gamma_ij*wd*rinv*rinv*dot*dot
+          uTmp = gamma_ij*wd*rinv*rinv*dot*dot
                 - 0.5*sigma[itype][jtype]*sigma[itype][jtype]*mu_ij*wd;
-	  uTmp -= sigma[itype][jtype]*wr*rinv*dot*randnum*dtinvsqrt;
-	  uTmp *= 0.5;
+          uTmp -= sigma[itype][jtype]*wr*rinv*dot*randnum*dtinvsqrt;
+          uTmp *= 0.5;

-	  duMech[i] += uTmp;
-	  if (newton_pair || j < nlocal) {
-	    duMech[j] += uTmp;
-	  }
-	
-	  // Compute uCond
-	  randnum = random->gaussian();
-	  kappa_ij = kappa[itype][jtype];
-	  alpha_ij = sqrt(2.0*force->boltz*kappa_ij);
-	  randPair = alpha_ij*wr*randnum*dtinvsqrt;
+          duMech[i] += uTmp;
+          if (newton_pair || j < nlocal) {
+            duMech[j] += uTmp;
+          }
+
+          // Compute uCond
+          randnum = random->gaussian();
+          kappa_ij = kappa[itype][jtype];
+          alpha_ij = sqrt(2.0*force->boltz*kappa_ij);
+          randPair = alpha_ij*wr*randnum*dtinvsqrt;
+
+          uTmp = kappa_ij*(1.0/dpdTheta[i] - 1.0/dpdTheta[j])*wd;
+          uTmp += randPair;
+
+          duCond[i] += uTmp;
+          if (newton_pair || j < nlocal) {
+            duCond[j] -= uTmp;
+          }

-	  uTmp = kappa_ij*(1.0/dpdTheta[i] - 1.0/dpdTheta[j])*wd;
-	  uTmp += randPair;
-	
-	  duCond[i] += uTmp;
-	  if (newton_pair || j < nlocal) {
-	    duCond[j] -= uTmp;
-	  }
-	
          if (eflag) {
            // unshifted eng of conservative term:
            // evdwl = -a0[itype][jtype]*r * (1.0-0.5*r/cut[itype][jtype]);
--- a/src/USER-DPD/pair_exp6_rx.cpp
+++ b/src/USER-DPD/pair_exp6_rx.cpp
@ -310,13 +310,13 @@ void PairExp6rx::compute(int eflag, int vflag)

            uin1 = buck1*(6.0*rin1exp - alphaOld12_ij*rm6ij*rin6inv) - urc - durc*(rin1-rCut);

-            win1 = -buck1*buck2*(rin1*rin1exp*rminv - rm6ij*rin6inv) - rin1*durc;
+            win1 = buck1*buck2*(rin1*rin1exp*rminv - rm6ij*rin6inv) + rin1*durc;

-            aRep = -1.0*win1*powint(rin1,nRep)/nRep;
+            aRep = win1*powint(rin1,nRep)/nRep;

            uin1rep = aRep/powint(rin1,nRep);

-            forceExp6 = -double(nRep)*aRep/powint(r,nRep);
+            forceExp6 = double(nRep)*aRep/powint(r,nRep);
            fpairOldEXP6_12 = factor_lj*forceExp6*r2inv;

            evdwlOldEXP6_12 = uin1 - uin1rep + aRep/powint(r,nRep);
@ -350,13 +350,13 @@ void PairExp6rx::compute(int eflag, int vflag)

            uin1 = buck1*(6.0*rin1exp - alphaOld21_ij*rm6ij*rin6inv) - urc - durc*(rin1-rCut);

-            win1 = -buck1*buck2*(rin1*rin1exp*rminv - rm6ij*rin6inv) - rin1*durc;
+            win1 = buck1*buck2*(rin1*rin1exp*rminv - rm6ij*rin6inv) + rin1*durc;

-            aRep = -1.0*win1*powint(rin1,nRep)/nRep;
+            aRep = win1*powint(rin1,nRep)/nRep;

            uin1rep = aRep/powint(rin1,nRep);

-            forceExp6 = -double(nRep)*aRep/powint(r,nRep);
+            forceExp6 = double(nRep)*aRep/powint(r,nRep);
            fpairOldEXP6_21 = factor_lj*forceExp6*r2inv;

            evdwlOldEXP6_21 = uin1 - uin1rep + aRep/powint(r,nRep);
@ -405,9 +405,9 @@ void PairExp6rx::compute(int eflag, int vflag)

            uin1 = buck1*(6.0*rin1exp - alpha12_ij*rm6ij*rin6inv) - urc - durc*(rin1-rCut);

-            win1 = -buck1*buck2*(rin1*rin1exp*rminv - rm6ij*rin6inv) - rin1*durc;
+            win1 = buck1*buck2*(rin1*rin1exp*rminv - rm6ij*rin6inv) + rin1*durc;

-            aRep = -1.0*win1*powint(rin1,nRep)/nRep;
+            aRep = win1*powint(rin1,nRep)/nRep;

            uin1rep = aRep/powint(rin1,nRep);

@ -437,9 +437,9 @@ void PairExp6rx::compute(int eflag, int vflag)

            uin1 = buck1*(6.0*rin1exp - alpha21_ij*rm6ij*rin6inv) - urc - durc*(rin1-rCut);

-            win1 = -buck1*buck2*(rin1*rin1exp*rminv - rm6ij*rin6inv) - rin1*durc;
+            win1 = buck1*buck2*(rin1*rin1exp*rminv - rm6ij*rin6inv) + rin1*durc;

-            aRep = -1.0*win1*powint(rin1,nRep)/nRep;
+            aRep = win1*powint(rin1,nRep)/nRep;

            uin1rep = aRep/powint(rin1,nRep);

@ -1102,6 +1102,32 @@ void PairExp6rx::getParamsEXP6(int id,double &epsilon1,double &alpha1,double &rm
      rm2_old *= pow(nTotalOFA_old,fuchslinR);
    }
  }
+
+  // Check that no fractions are less than zero
+  if(fraction1 < 0.0){
+    if(fraction1 < -1.0e-10){
+      error->all(FLERR,"Computed fraction less than -1.0e-10");
+    }
+    fraction1 = 0.0;
+  }
+  if(fraction2 < 0.0){
+    if(fraction2 < -1.0e-10){
+      error->all(FLERR,"Computed fraction less than -1.0e-10");
+    }
+    fraction2 = 0.0;
+  }
+  if(fraction1_old < 0.0){
+    if(fraction1_old < -1.0e-10){
+      error->all(FLERR,"Computed fraction less than -1.0e-10");
+    }
+    fraction1_old = 0.0;
+  }
+  if(fraction2_old < 0.0){
+    if(fraction2_old < -1.0e-10){
+      error->all(FLERR,"Computed fraction less than -1.0e-10");
+    }
+    fraction2_old = 0.0;
+  }
 }

 /* ---------------------------------------------------------------------- */
--- a/src/USER-DPD/pair_multi_lucy.cpp
+++ b/src/USER-DPD/pair_multi_lucy.cpp
@ -182,7 +182,7 @@ void PairMultiLucy::compute(int eflag, int vflag)
          f[j][2] -= delz*fpair;
        }
        if (evflag) ev_tally(i,j,nlocal,newton_pair,
-        		     0.0,0.0,fpair,delx,dely,delz);
+                             0.0,0.0,fpair,delx,dely,delz);
      }
    }

@ -201,7 +201,7 @@ void PairMultiLucy::compute(int eflag, int vflag)
    evdwl *=(pi*cutsq[itype][itype]*cutsq[itype][itype])/84.0;

    if (evflag) ev_tally(0,0,nlocal,newton_pair,
-        		 evdwl,0.0,0.0,0.0,0.0,0.0);
+                         evdwl,0.0,0.0,0.0,0.0,0.0);
  }

 if (vflag_fdotr) virial_fdotr_compute();
--- a/src/USER-DPD/pair_multi_lucy.h
+++ b/src/USER-DPD/pair_multi_lucy.h
@ -18,7 +18,7 @@ PairStyle(multi/lucy,PairMultiLucy)
 #else

 #ifndef LMP_PAIR_MULTI_LUCY_H
-#define LMP_PAIR_MUTLI_LUCY_H
+#define LMP_PAIR_MULTI_LUCY_H

 #include "pair.h"

--- a/src/USER-DPD/pair_multi_lucy_rx.h
+++ b/src/USER-DPD/pair_multi_lucy_rx.h
@ -18,7 +18,7 @@ PairStyle(multi/lucy/rx,PairMultiLucyRX)
 #else

 #ifndef LMP_PAIR_MULTI_LUCY_RX_H
-#define LMP_PAIR_MUTLI_LUCY_RX_H
+#define LMP_PAIR_MULTI_LUCY_RX_H

 #include "pair.h"

--- a/src/USER-INTEL/fix_intel.cpp
+++ b/src/USER-INTEL/fix_intel.cpp
@ -325,8 +325,6 @@ void FixIntel::init()
    error->all(FLERR,
 	       "Currently, cannot use more than one intel style with hybrid.");

-  neighbor->fix_intel = (void *)this;
-
  check_neighbor_intel();
  if (_precision_mode == PREC_MODE_SINGLE)
    _single_buffers->zero_ev();
--- a/src/atom_vec_body.cpp
+++ b/src/atom_vec_body.cpp
@ -847,6 +847,7 @@ void AtomVecBody::unpack_border(int n, int first, double *buf)
      inertia[2] = buf[m++];
      bonus[j].ninteger = (int) ubuf(buf[m++]).i;
      bonus[j].ndouble = (int) ubuf(buf[m++]).i;
+      // corresponding put() calls are in clear_bonus()
      bonus[j].ivalue = icp->get(bonus[j].ninteger,bonus[j].iindex);
      bonus[j].dvalue = dcp->get(bonus[j].ndouble,bonus[j].dindex);
      m += bptr->unpack_border_body(&bonus[j],&buf[m]);
@ -897,6 +898,7 @@ void AtomVecBody::unpack_border_vel(int n, int first, double *buf)
      inertia[2] = buf[m++];
      bonus[j].ninteger = (int) ubuf(buf[m++]).i;
      bonus[j].ndouble = (int) ubuf(buf[m++]).i;
+      // corresponding put() calls are in clear_bonus()
      bonus[j].ivalue = icp->get(bonus[j].ninteger,bonus[j].iindex);
      bonus[j].dvalue = dcp->get(bonus[j].ndouble,bonus[j].dindex);
      m += bptr->unpack_border_body(&bonus[j],&buf[m]);
@ -946,6 +948,7 @@ int AtomVecBody::unpack_border_hybrid(int n, int first, double *buf)
      inertia[2] = buf[m++];
      bonus[j].ninteger = (int) ubuf(buf[m++]).i;
      bonus[j].ndouble = (int) ubuf(buf[m++]).i;
+      // corresponding put() calls are in clear_bonus()
      bonus[j].ivalue = icp->get(bonus[j].ninteger,bonus[j].iindex);
      bonus[j].dvalue = dcp->get(bonus[j].ndouble,bonus[j].dindex);
      m += bptr->unpack_border_body(&bonus[j],&buf[m]);
@ -1050,6 +1053,7 @@ int AtomVecBody::unpack_exchange(double *buf)
    inertia[2] = buf[m++];
    bonus[nlocal_bonus].ninteger = (int) ubuf(buf[m++]).i;
    bonus[nlocal_bonus].ndouble = (int) ubuf(buf[m++]).i;
+    // corresponding put() calls are in copy()
    bonus[nlocal_bonus].ivalue = icp->get(bonus[nlocal_bonus].ninteger,
 					  bonus[nlocal_bonus].iindex);
    bonus[nlocal_bonus].dvalue = dcp->get(bonus[nlocal_bonus].ndouble,
--- a/src/fix_shear_history.cpp
+++ b/src/fix_shear_history.cpp
@ -526,7 +526,7 @@ double FixShearHistory::memory_usage()
 {
  int nmax = atom->nmax;
  double bytes = nmax * sizeof(int);
-  bytes += nmax * sizeof(int *);
+  bytes += nmax * sizeof(tagint *);
  bytes += nmax * sizeof(double *);

  int nmypage = comm->nthreads;
--- a/src/memory.cpp
+++ b/src/memory.cpp
@ -26,6 +26,10 @@
 #endif
 #endif

+#if defined(LMP_USER_INTEL) && !defined(LAMMPS_MEMALIGN)
+#define LAMMPS_MEMALIGN 64
+#endif
+
 using namespace LAMMPS_NS;

 /* ---------------------------------------------------------------------- */
--- a/src/my_page.h
+++ b/src/my_page.h
@ -57,10 +57,10 @@ namespace LAMMPS_NS {

 template<class T>
 class MyPage {
+ public:
  int ndatum;      // total # of stored datums
  int nchunk;      // total # of stored chunks

- public:
  MyPage() {
    ndatum = nchunk = 0;
    pages = NULL;
--- a/src/my_pool_chunk.h
+++ b/src/my_pool_chunk.h
@ -33,7 +33,7 @@ methods:
                     minchunk <= N <= maxchunk required
   put(index) = return indexed chunk to pool (same index returned by get)
   int size() = return total size of allocated pages in bytes
-public varaibles:
+public variables:
   ndatum = total # of stored datums
   nchunk = total # of stored chunks
   size = total size of all allocated pages in daums
--- a/src/neigh_list.cpp
+++ b/src/neigh_list.cpp
@ -192,6 +192,7 @@ void NeighList::setup_pages(int pgsize_caller, int oneatom_caller)
     gran calls grow() in granhistory
     respaouter calls grow() in respainner, respamiddle
   triggered by neighbor list build
+   not called if a copy list
 ------------------------------------------------------------------------- */

 void NeighList::grow(int nlocal, int nall)
--- a/src/neighbor.cpp
+++ b/src/neighbor.cpp
@ -1369,6 +1369,8 @@ int Neighbor::choose_pair(NeighRequest *rq)
               "with ghost neighbors");

  // flags for settings the request + system requires of NPair class
+  // some are set to 0/1, others are set to mask bit
+  // comparisons below in loop over classes reflect that
  //   copyflag = no/yes copy request
  //   skipflag = no/yes skip request
  //   halfflag = half request (gran and respa are also half lists)
@ -1402,7 +1404,7 @@ int Neighbor::choose_pair(NeighRequest *rq)
  // NOTE: exactly one of these request flags is set (see neigh_request.h)
  //       this requires gran/respaouter also set halfflag
  //       can simplify this logic, if follow NOTE in neigh_request.h
-  //       all why do size/off2on and size/off2on/oneside set NP_HALF
+  //       why do size/off2on and size/off2on/oneside set NP_HALF
  //         either should set both half & full, or half should be in file name
  //         to be consistent with how other NP classes use "half"

@ -1433,9 +1435,8 @@ int Neighbor::choose_pair(NeighRequest *rq)
  else if (rq->newton == 1) newtflag = 1;
  else if (rq->newton == 2) newtflag = 0;

-  // use flags to match exactly one of NPair class masks, bit by bit
-  // copyflag match returns with no further checks
-  // exactly one of halfflag,fullflag,halffullflag is set and thus must match
+  // use flags to match exactly one of NPair class masks
+  // sequence of checks is bit by bit in NeighConst

  int mask;

@ -1447,13 +1448,21 @@ int Neighbor::choose_pair(NeighRequest *rq)
  for (int i = 0; i < npclass; i++) {
    mask = pairmasks[i];

-    if (copyflag && (mask & NP_COPY)) {
+    // if copyflag set, return or continue with no further checks
+
+    if (copyflag) {
+      if (!(mask & NP_COPY)) continue;
      if (kokkos_device_flag != (mask & NP_KOKKOS_DEVICE)) continue;
      if (kokkos_host_flag != (mask & NP_KOKKOS_HOST)) continue;
      return i+1;
    }
+
+    // skipflag must match along with other flags, so do not return
+
    if (skipflag != (mask & NP_SKIP)) continue;

+    // exactly one of halfflag,fullflag,halffullflag is set and must match
+
    if (halfflag) {
      if (!(mask & NP_HALF)) continue;
    } else if (fullflag) {
@ -1470,19 +1479,38 @@ int Neighbor::choose_pair(NeighRequest *rq)
    if (ssaflag != (mask & NP_SSA)) continue;
    if (ompflag != (mask & NP_OMP)) continue;
    if (intelflag != (mask & NP_INTEL)) continue;
+
+    // style is one of NSQ,BIN,MULTI and must match
+
+    if (style == NSQ) {
+      if (!(mask & NP_NSQ)) continue;
+    } else if (style == BIN) {
+      if (!(mask & NP_BIN)) continue;
+    } else if (style == MULTI) {
+      if (!(mask & NP_MULTI)) continue;
+    }
+
+    // newtflag is on or off and must match
+
+    if (newtflag) {
+      if (!(mask & NP_NEWTON)) continue;
+    } else if (!newtflag) {
+      if (!(mask & NP_NEWTOFF)) continue;
+    }
+
+    // triclinic flag is on or off and must match
+
+    if (triclinic) {
+      if (!(mask & NP_TRI)) continue;
+    } else if (!triclinic) {
+      if (!(mask & NP_ORTHO)) continue;
+    }
+
+    // Kokkos flags
+
    if (kokkos_device_flag != (mask & NP_KOKKOS_DEVICE)) continue;
    if (kokkos_host_flag != (mask & NP_KOKKOS_HOST)) continue;

-    if (style == NSQ && !(mask & NP_NSQ)) continue;
-    if (style == BIN && !(mask & NP_BIN)) continue;
-    if (style == MULTI && !(mask & NP_MULTI)) continue;
-
-    if (newtflag && !(mask & NP_NEWTON)) continue;
-    if (!newtflag && !(mask & NP_NEWTOFF)) continue;
-
-    if (!triclinic && !(mask & NP_ORTHO)) continue;
-    if (triclinic && !(mask & NP_TRI)) continue;
-
    return i+1;
  }

@ -1727,7 +1755,7 @@ void Neighbor::build(int topoflag)

  for (i = 0; i < npair_perpetual; i++) {
    m = plist[i];
-    lists[m]->grow(nlocal,nall);
+    if (!lists[m]->copy) lists[m]->grow(nlocal,nall);
    neigh_pair[m]->build_setup();
    neigh_pair[m]->build(lists[m]);
  }
--- a/src/npair_skip.h
+++ b/src/npair_skip.h
@ -15,7 +15,8 @@

 NPairStyle(skip,
           NPairSkip,
-           NP_SKIP | NP_HALF | NP_FULL | NP_NSQ | NP_BIN | NP_MULTI | 
+           NP_SKIP | NP_HALF | NP_FULL | NP_HALFFULL |
+           NP_NSQ | NP_BIN | NP_MULTI | 
           NP_NEWTON | NP_NEWTOFF | NP_ORTHO | NP_TRI)

 #else
--- a/src/python_wrapper.h
+++ b/src/python_wrapper.h
@ -38,7 +38,7 @@ class Python {
  void invoke_function(int, char *) {}
  int find(char *) {return -1;}
  int variable_match(char *, char *, int) {return -1;}
-
+  char *long_string(int) {return NULL;}
 };

 }
--- a/src/variable.cpp
+++ b/src/variable.cpp
@ -875,6 +875,10 @@ char *Variable::retrieve(char *name)
      error->all(FLERR,"Python variable does not match Python function");
    python->invoke_function(ifunc,data[ivar][1]);
    str = data[ivar][1];
+    // if Python func returns a string longer than VALUELENGTH
+    // then the Python class stores the result, query it via long_string()
+    char *strlong = python->long_string(ifunc);
+    if (strlong) str = strlong;
  } else if (style[ivar] == INTERNAL) {
    sprintf(data[ivar][0],"%.15g",dvalue[ivar]);
    str = data[ivar][0];
--- a/src/version.h
+++ b/src/version.h
@ -1 +1 @@
-#define LAMMPS_VERSION "13 Dec 2016"
+#define LAMMPS_VERSION "21 Dec 2016"
Author	SHA1	Message	Date
Steve Plimpton	62dea1bb63	21Dec16 patch	2016-12-21 09:53:32 -07:00
sjplimp	800ff43413	Merge pull request #304 from timattox/USER-DPD_whitespace USER-DPD: Whitespace cleanup to pair_dpd_fdt_energy.cpp	2016-12-21 09:51:10 -07:00
Steve Plimpton	9161bd98bf	fixed bug with pair hybrid/overaly and manybody potentials finding the right skip neighbor method	2016-12-21 09:50:29 -07:00
Tim Mattox	c4e02a5d2b	USER-DPD: more whitespace fixes	2016-12-20 11:17:11 -05:00
Tim Mattox	2ba424e1a3	USER-DPD: Whitespace cleanup to pair_dpd_fdt_energy.cpp	2016-12-19 15:11:52 -05:00
Steve Plimpton	81a2db8a0c	17Dec16 patch	2016-12-16 11:36:54 -07:00
Steve Plimpton	0a176841e7	extra python_wrapper change needed for last patch	2016-12-16 11:35:42 -07:00
Steve Plimpton	3027ac9250	patch 16Dec16	2016-12-16 10:30:57 -07:00
sjplimp	fc54ab5cea	Merge pull request #301 from akohlmey/corrections-and-bugfixes Collected corrections and bugfixes	2016-12-16 10:25:29 -07:00
Steve Plimpton	e364b80724	added length keyword to python command	2016-12-16 10:24:25 -07:00
Axel Kohlmeyer	830c9e8661	Merge branch 'USER-DPD_internal_energy' of https://github.com/timattox/lammps_USER-DPD into corrections-and-bugfixes This closes #303	2016-12-16 11:22:25 -05:00
Axel Kohlmeyer	4907b29ad2	Merge branch 'USER-DPD_bugfixes' of https://github.com/timattox/lammps_USER-DPD into corrections-and-bugfixes This closes #302	2016-12-16 11:21:15 -05:00
Tim Mattox	eff7238ff2	USER-DPD: fix_eos*: partition all internal energy into the uMech term only This makes our results more closely match a vetted serial implementation. NOTE: This does make the output different from any previous versions. Patch by Jim Larentzos. Applied by Tim Mattox.	2016-12-16 10:25:12 -05:00
Tim Mattox	126fb22e93	USER-DPD: Fix #define typo in pair_multi_lucy.h and pair_multi_lucy_rx.h	2016-12-16 10:08:30 -05:00
Tim Mattox	0a90492c44	USER-DPD: Update the header files to properly document all error statements Patch by Jim Larentzos. Applied by Tim Mattox.	2016-12-15 17:39:15 -05:00
Tim Mattox	fed629c23e	USER-DPD: Bugfix for fix_rx and fix_eos_table_rx to handle restart files. Patch by Jim Larentzos. Applied by Tim Mattox.	2016-12-15 17:10:13 -05:00
Tim Mattox	925481c3f4	USER-DPD: Fix hard-wall force interaction bug, and ensure fraction is >= 0 pair_exp6_rx.cpp patch by Jim Larentzos. Applied by Tim Mattox.	2016-12-15 16:46:25 -05:00
Axel Kohlmeyer	da2ad5b6e0	update FixIntel code for new neighbor list code	2016-12-14 15:51:12 -05:00
Axel Kohlmeyer	bfcab72268	restore change to make -DLAMMPS_MEMALIGN=64 default when USER-INTEL package is installed (which requires it)	2016-12-14 15:24:55 -05:00