diff --git a/src/USER-DPD/README b/src/USER-DPD/README
new file mode 100644
index 0000000000000000000000000000000000000000..1626375c4b6531363055bfb58f4c10ed90467a1b
--- /dev/null
+++ b/src/USER-DPD/README
@@ -0,0 +1,35 @@
+This package implements the dissipative particle dynamics (DPD) method
+under isothermal, isoenergetic, isobaric and isenthalpic conditions.
+The DPD equations of motion are integrated through the Shardlow
+splitting algorithm.
+
+Currently, the package has the following features:
+
+* A new DPD atom style for tracking the DPD particle internal energies
+ and internal temperature
+
+* Compute commands for accessing the DPD particle internal energies
+ and internal temperature
+
+* "fix eos" commands for relating the DPD internal energy to the DPD
+ internal temperature through a coarse-grained particle
+ equation-of-state
+
+* "fix shardlow" command for integrating the stochastic ODEs
+
+* Pair styles for modeling a DPD fluid
+
+* Commands for setting the particle internal temperature
+
+See the doc pages for "atom style dpd", "compute dpd" and "compute
+dpd/atom", "fix eos/cv" and "fix eos/table", "fix shardlow", "pair
+dpd/conservative" and "pair dpd/fdt" and "pair dpd/fdt/energy"
+commands to get started. At the bottom of the doc page are many links
+to additional documentation contained in the doc/USER/dpd directory.
+
+There are example scripts for using this package in examples/USER/dpd.
+
+The primary people who created this package are James Larentzos
+(james.p.larentzos.civ at mail.mil), Timothy Mattox (Timothy.Mattox at
+engilitycorp.com) and John Brennan (john.k.brennan.civ at mail.mil).
+Contact them directly if you have questions.
diff --git a/src/USER-DPD/atom_vec_dpd.cpp b/src/USER-DPD/atom_vec_dpd.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..7f48293389e6937b9acccd334bd6fbd330b076c4
--- /dev/null
+++ b/src/USER-DPD/atom_vec_dpd.cpp
@@ -0,0 +1,871 @@
+/* ----------------------------------------------------------------------
+ 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: James Larentzos (Engility Corporation)
+------------------------------------------------------------------------- */
+
+#include <stdlib.h>
+#include "atom_vec_dpd.h"
+#include "atom.h"
+#include "comm.h"
+#include "domain.h"
+#include "modify.h"
+#include "fix.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+AtomVecDPD::AtomVecDPD(LAMMPS *lmp) : AtomVec(lmp)
+{
+ molecular = 0;
+ mass_type = 1;
+
+ comm_x_only = comm_f_only = 0; // we communicate not only x forward but also dpdTheta
+ size_forward = 6; // 3 + dpdTheta + uCond + uMech
+ size_reverse = 3; // 3
+ size_border = 9; // 6 + dpdTheta + uCond + uMech
+ size_velocity = 3;
+ size_data_atom = 6; // we read id + type + dpdTheta + x + y + z
+ size_data_vel = 4;
+ xcol_data = 4; // 1=id 2=type 3=dpdTheta 4=x
+
+ atom->dpd_flag = 1;
+}
+
+/* ----------------------------------------------------------------------
+ grow atom arrays
+ n = 0 grows arrays by a chunk
+ n > 0 allocates arrays to size n
+------------------------------------------------------------------------- */
+
+void AtomVecDPD::grow(int n)
+{
+ if (n == 0) grow_nmax();
+ else nmax = n;
+ atom->nmax = nmax;
+ if (nmax < 0)
+ error->one(FLERR,"Per-processor system is too big");
+
+ tag = memory->grow(atom->tag,nmax,"atom:tag");
+ type = memory->grow(atom->type,nmax,"atom:type");
+ mask = memory->grow(atom->mask,nmax,"atom:mask");
+ image = memory->grow(atom->image,nmax,"atom:image");
+ x = memory->grow(atom->x,nmax,3,"atom:x");
+ v = memory->grow(atom->v,nmax,3,"atom:v");
+ f = memory->grow(atom->f,nmax*comm->nthreads,3,"atom:f");
+
+ dpdTheta = memory->grow(atom->dpdTheta, nmax, "atom:dpdTheta");
+ uCond = memory->grow(atom->uCond,nmax,"atom:uCond");
+ uMech = memory->grow(atom->uMech,nmax,"atom:uMech");
+ duCond = memory->grow(atom->duCond,nmax,"atom:duCond");
+ duMech = memory->grow(atom->duMech,nmax,"atom:duMech");
+
+ if (atom->nextra_grow)
+ for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
+ modify->fix[atom->extra_grow[iextra]]->grow_arrays(nmax);
+}
+
+/* ----------------------------------------------------------------------
+ reset local array ptrs
+------------------------------------------------------------------------- */
+
+void AtomVecDPD::grow_reset()
+{
+ tag = atom->tag; type = atom->type;
+ mask = atom->mask; image = atom->image;
+ x = atom->x; v = atom->v; f = atom->f;
+ dpdTheta = atom->dpdTheta;
+ uCond = atom->uCond;
+ uMech = atom->uMech;
+ duCond = atom->duCond;
+ duMech = atom->duMech;
+}
+
+/* ----------------------------------------------------------------------
+ copy atom I info to atom J
+------------------------------------------------------------------------- */
+
+void AtomVecDPD::copy(int i, int j, int delflag)
+{
+ tag[j] = tag[i];
+ type[j] = type[i];
+ mask[j] = mask[i];
+ image[j] = image[i];
+ x[j][0] = x[i][0];
+ x[j][1] = x[i][1];
+ x[j][2] = x[i][2];
+ v[j][0] = v[i][0];
+ v[j][1] = v[i][1];
+ v[j][2] = v[i][2];
+ dpdTheta[j] = dpdTheta[i];
+ uCond[j] = uCond[i];
+ uMech[j] = uMech[i];
+
+ if (atom->nextra_grow)
+ for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
+ modify->fix[atom->extra_grow[iextra]]->copy_arrays(i,j,delflag);
+}
+
+/* ---------------------------------------------------------------------- */
+
+int AtomVecDPD::pack_comm(int n, int *list, double *buf,
+ int pbc_flag, int *pbc)
+{
+ int i,j,m;
+ double dx,dy,dz;
+
+ m = 0;
+ if (pbc_flag == 0) {
+ for (i = 0; i < n; i++) {
+ j = list[i];
+ buf[m++] = x[j][0];
+ buf[m++] = x[j][1];
+ buf[m++] = x[j][2];
+ buf[m++] = dpdTheta[j];
+ buf[m++] = uCond[j];
+ buf[m++] = uMech[j];
+ }
+ } else {
+ if (domain->triclinic == 0) {
+ dx = pbc[0]*domain->xprd;
+ dy = pbc[1]*domain->yprd;
+ dz = pbc[2]*domain->zprd;
+ } else {
+ dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz;
+ dy = pbc[1]*domain->yprd + pbc[3]*domain->yz;
+ dz = pbc[2]*domain->zprd;
+ }
+ for (i = 0; i < n; i++) {
+ j = list[i];
+ buf[m++] = x[j][0] + dx;
+ buf[m++] = x[j][1] + dy;
+ buf[m++] = x[j][2] + dz;
+ buf[m++] = dpdTheta[j];
+ buf[m++] = uCond[j];
+ buf[m++] = uMech[j];
+ }
+ }
+ return m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int AtomVecDPD::pack_comm_vel(int n, int *list, double *buf,
+ int pbc_flag, int *pbc)
+{
+ int i,j,m;
+ double dx,dy,dz,dvx,dvy,dvz;
+
+ m = 0;
+ if (pbc_flag == 0) {
+ for (i = 0; i < n; i++) {
+ j = list[i];
+ buf[m++] = x[j][0];
+ buf[m++] = x[j][1];
+ buf[m++] = x[j][2];
+ buf[m++] = v[j][0];
+ buf[m++] = v[j][1];
+ buf[m++] = v[j][2];
+ buf[m++] = dpdTheta[j];
+ buf[m++] = uCond[j];
+ buf[m++] = uMech[j];
+ }
+ } else {
+ if (domain->triclinic == 0) {
+ dx = pbc[0]*domain->xprd;
+ dy = pbc[1]*domain->yprd;
+ dz = pbc[2]*domain->zprd;
+ } else {
+ dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz;
+ dy = pbc[1]*domain->yprd + pbc[3]*domain->yz;
+ dz = pbc[2]*domain->zprd;
+ }
+ if (!deform_vremap) {
+ for (i = 0; i < n; i++) {
+ j = list[i];
+ buf[m++] = x[j][0] + dx;
+ buf[m++] = x[j][1] + dy;
+ buf[m++] = x[j][2] + dz;
+ buf[m++] = v[j][0];
+ buf[m++] = v[j][1];
+ buf[m++] = v[j][2];
+ buf[m++] = dpdTheta[j];
+ buf[m++] = uCond[j];
+ buf[m++] = uMech[j];
+ }
+ } else {
+ dvx = pbc[0]*h_rate[0] + pbc[5]*h_rate[5] + pbc[4]*h_rate[4];
+ dvy = pbc[1]*h_rate[1] + pbc[3]*h_rate[3];
+ dvz = pbc[2]*h_rate[2];
+ for (i = 0; i < n; i++) {
+ j = list[i];
+ buf[m++] = x[j][0] + dx;
+ buf[m++] = x[j][1] + dy;
+ buf[m++] = x[j][2] + dz;
+ if (mask[i] & deform_groupbit) {
+ buf[m++] = v[j][0] + dvx;
+ buf[m++] = v[j][1] + dvy;
+ buf[m++] = v[j][2] + dvz;
+ } else {
+ buf[m++] = v[j][0];
+ buf[m++] = v[j][1];
+ buf[m++] = v[j][2];
+ }
+ buf[m++] = dpdTheta[j];
+ buf[m++] = uCond[j];
+ buf[m++] = uMech[j];
+ }
+ }
+ }
+ return m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void AtomVecDPD::unpack_comm(int n, int first, double *buf)
+{
+ int i,m,last;
+
+ m = 0;
+ last = first + n;
+ for (i = first; i < last; i++) {
+ x[i][0] = buf[m++];
+ x[i][1] = buf[m++];
+ x[i][2] = buf[m++];
+ dpdTheta[i] = buf[m++];
+ uCond[i] = buf[m++];
+ uMech[i] = buf[m++];
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void AtomVecDPD::unpack_comm_vel(int n, int first, double *buf)
+{
+ int i,m,last;
+
+ m = 0;
+ last = first + n;
+ for (i = first; i < last; i++) {
+ x[i][0] = buf[m++];
+ x[i][1] = buf[m++];
+ x[i][2] = buf[m++];
+ v[i][0] = buf[m++];
+ v[i][1] = buf[m++];
+ v[i][2] = buf[m++];
+ dpdTheta[i] = buf[m++];
+ uCond[i] = buf[m++];
+ uMech[i] = buf[m++];
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+int AtomVecDPD::pack_reverse(int n, int first, double *buf)
+{
+ int i,m,last;
+
+ m = 0;
+ last = first + n;
+ for (i = first; i < last; i++) {
+ buf[m++] = f[i][0];
+ buf[m++] = f[i][1];
+ buf[m++] = f[i][2];
+ }
+ return m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void AtomVecDPD::unpack_reverse(int n, int *list, double *buf)
+{
+ int i,j,m;
+
+ m = 0;
+ for (i = 0; i < n; i++) {
+ j = list[i];
+ f[j][0] += buf[m++];
+ f[j][1] += buf[m++];
+ f[j][2] += buf[m++];
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+int AtomVecDPD::pack_border(int n, int *list, double *buf,
+ int pbc_flag, int *pbc)
+{
+ int i,j,m;
+ double dx,dy,dz;
+
+ m = 0;
+ if (pbc_flag == 0) {
+ for (i = 0; i < n; i++) {
+ j = list[i];
+ buf[m++] = x[j][0];
+ buf[m++] = x[j][1];
+ buf[m++] = x[j][2];
+ buf[m++] = ubuf(tag[j]).d;
+ buf[m++] = ubuf(type[j]).d;
+ buf[m++] = ubuf(mask[j]).d;
+ buf[m++] = dpdTheta[j];
+ buf[m++] = uCond[j];
+ buf[m++] = uMech[j];
+ }
+ } else {
+ if (domain->triclinic == 0) {
+ dx = pbc[0]*domain->xprd;
+ dy = pbc[1]*domain->yprd;
+ dz = pbc[2]*domain->zprd;
+ } else {
+ dx = pbc[0];
+ dy = pbc[1];
+ dz = pbc[2];
+ }
+ for (i = 0; i < n; i++) {
+ j = list[i];
+ buf[m++] = x[j][0] + dx;
+ buf[m++] = x[j][1] + dy;
+ buf[m++] = x[j][2] + dz;
+ buf[m++] = ubuf(tag[j]).d;
+ buf[m++] = ubuf(type[j]).d;
+ buf[m++] = ubuf(mask[j]).d;
+ buf[m++] = dpdTheta[j];
+ buf[m++] = uCond[j];
+ buf[m++] = uMech[j];
+ }
+ }
+
+ if (atom->nextra_border)
+ for (int iextra = 0; iextra < atom->nextra_border; iextra++)
+ m += modify->fix[atom->extra_border[iextra]]->pack_border(n,list,&buf[m]);
+
+ return m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int AtomVecDPD::pack_border_vel(int n, int *list, double *buf,
+ int pbc_flag, int *pbc)
+{
+ int i,j,m;
+ double dx,dy,dz,dvx,dvy,dvz;
+
+ m = 0;
+ if (pbc_flag == 0) {
+ for (i = 0; i < n; i++) {
+ j = list[i];
+ buf[m++] = x[j][0];
+ buf[m++] = x[j][1];
+ buf[m++] = x[j][2];
+ buf[m++] = ubuf(tag[j]).d;
+ buf[m++] = ubuf(type[j]).d;
+ buf[m++] = ubuf(mask[j]).d;
+ buf[m++] = v[j][0];
+ buf[m++] = v[j][1];
+ buf[m++] = v[j][2];
+ buf[m++] = dpdTheta[j];
+ buf[m++] = uCond[j];
+ buf[m++] = uMech[j];
+ }
+ } else {
+ if (domain->triclinic == 0) {
+ dx = pbc[0]*domain->xprd;
+ dy = pbc[1]*domain->yprd;
+ dz = pbc[2]*domain->zprd;
+ } else {
+ dx = pbc[0];
+ dy = pbc[1];
+ dz = pbc[2];
+ }
+ if (!deform_vremap) {
+ for (i = 0; i < n; i++) {
+ j = list[i];
+ buf[m++] = x[j][0] + dx;
+ buf[m++] = x[j][1] + dy;
+ buf[m++] = x[j][2] + dz;
+ buf[m++] = ubuf(tag[j]).d;
+ buf[m++] = ubuf(type[j]).d;
+ buf[m++] = ubuf(mask[j]).d;
+ buf[m++] = v[j][0];
+ buf[m++] = v[j][1];
+ buf[m++] = v[j][2];
+ buf[m++] = dpdTheta[j];
+ buf[m++] = uCond[j];
+ buf[m++] = uMech[j];
+ }
+ } else {
+ dvx = pbc[0]*h_rate[0] + pbc[5]*h_rate[5] + pbc[4]*h_rate[4];
+ dvy = pbc[1]*h_rate[1] + pbc[3]*h_rate[3];
+ dvz = pbc[2]*h_rate[2];
+ for (i = 0; i < n; i++) {
+ j = list[i];
+ buf[m++] = x[j][0] + dx;
+ buf[m++] = x[j][1] + dy;
+ buf[m++] = x[j][2] + dz;
+ buf[m++] = ubuf(tag[j]).d;
+ buf[m++] = ubuf(type[j]).d;
+ buf[m++] = ubuf(mask[j]).d;
+ if (mask[i] & deform_groupbit) {
+ buf[m++] = v[j][0] + dvx;
+ buf[m++] = v[j][1] + dvy;
+ buf[m++] = v[j][2] + dvz;
+ } else {
+ buf[m++] = v[j][0];
+ buf[m++] = v[j][1];
+ buf[m++] = v[j][2];
+ }
+ buf[m++] = dpdTheta[j];
+ buf[m++] = uCond[j];
+ buf[m++] = uMech[j];
+ }
+ }
+ }
+
+ if (atom->nextra_border)
+ for (int iextra = 0; iextra < atom->nextra_border; iextra++)
+ m += modify->fix[atom->extra_border[iextra]]->pack_border(n,list,&buf[m]);
+
+ return m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int AtomVecDPD::pack_comm_hybrid(int n, int *list, double *buf)
+{
+ int i,j,m;
+
+ m = 0;
+ for (i = 0; i < n; i++) {
+ j = list[i];
+ buf[m++] = dpdTheta[j];
+ buf[m++] = uCond[j];
+ buf[m++] = uMech[j];
+ }
+ return m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int AtomVecDPD::pack_border_hybrid(int n, int *list, double *buf)
+{
+ int i,j,m;
+
+ m = 0;
+ for (i = 0; i < n; i++) {
+ j = list[i];
+ buf[m++] = dpdTheta[j];
+ buf[m++] = uCond[j];
+ buf[m++] = uMech[j];
+ }
+ return m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void AtomVecDPD::unpack_border(int n, int first, double *buf)
+{
+ int i,m,last;
+
+ m = 0;
+ last = first + n;
+ for (i = first; i < last; i++) {
+ if (i == nmax) grow(0);
+ x[i][0] = buf[m++];
+ x[i][1] = buf[m++];
+ x[i][2] = buf[m++];
+ tag[i] = (tagint) ubuf(buf[m++]).i;
+ type[i] = (int) ubuf(buf[m++]).i;
+ mask[i] = (int) ubuf(buf[m++]).i;
+ dpdTheta[i] = buf[m++];
+ uCond[i] = buf[m++];
+ uMech[i] = buf[m++];
+ }
+
+ if (atom->nextra_border)
+ for (int iextra = 0; iextra < atom->nextra_border; iextra++)
+ m += modify->fix[atom->extra_border[iextra]]->
+ unpack_border(n,first,&buf[m]);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void AtomVecDPD::unpack_border_vel(int n, int first, double *buf)
+{
+ int i,m,last;
+
+ m = 0;
+ last = first + n;
+ for (i = first; i < last; i++) {
+ if (i == nmax) grow(0);
+ x[i][0] = buf[m++];
+ x[i][1] = buf[m++];
+ x[i][2] = buf[m++];
+ tag[i] = (tagint) ubuf(buf[m++]).i;
+ type[i] = (int) ubuf(buf[m++]).i;
+ mask[i] = (int) ubuf(buf[m++]).i;
+ v[i][0] = buf[m++];
+ v[i][1] = buf[m++];
+ v[i][2] = buf[m++];
+ dpdTheta[i] = buf[m++];
+ uCond[i] = buf[m++];
+ uMech[i] = buf[m++];
+ }
+
+ if (atom->nextra_border)
+ for (int iextra = 0; iextra < atom->nextra_border; iextra++)
+ m += modify->fix[atom->extra_border[iextra]]->
+ unpack_border(n,first,&buf[m]);
+}
+
+/* ---------------------------------------------------------------------- */
+
+int AtomVecDPD::unpack_comm_hybrid(int n, int first, double *buf)
+{
+ int i,m,last;
+
+ m = 0;
+ last = first + n;
+ for (i = first; i < last; i++) {
+ dpdTheta[i] = buf[m++];
+ uCond[i] = buf[m++];
+ uMech[i] = buf[m++];
+ }
+ return m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int AtomVecDPD::unpack_border_hybrid(int n, int first, double *buf)
+{
+ int i,m,last;
+
+ m = 0;
+ last = first + n;
+ for (i = first; i < last; i++) {
+ dpdTheta[i] = buf[m++];
+ uCond[i] = buf[m++];
+ uMech[i] = buf[m++];
+ }
+ return m;
+}
+
+/* ----------------------------------------------------------------------
+ pack data for atom I for sending to another proc
+ xyz must be 1st 3 values, so comm::exchange() can test on them
+------------------------------------------------------------------------- */
+
+int AtomVecDPD::pack_exchange(int i, double *buf)
+{
+ int m = 1;
+ buf[m++] = x[i][0];
+ buf[m++] = x[i][1];
+ buf[m++] = x[i][2];
+ buf[m++] = v[i][0];
+ buf[m++] = v[i][1];
+ buf[m++] = v[i][2];
+ buf[m++] = ubuf(tag[i]).d;
+ buf[m++] = ubuf(type[i]).d;
+ buf[m++] = ubuf(mask[i]).d;
+ buf[m++] = ubuf(image[i]).d;
+ buf[m++] = dpdTheta[i];
+ buf[m++] = uCond[i];
+ buf[m++] = uMech[i];
+
+ if (atom->nextra_grow)
+ for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
+ m += modify->fix[atom->extra_grow[iextra]]->pack_exchange(i,&buf[m]);
+
+ buf[0] = m;
+ return m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int AtomVecDPD::unpack_exchange(double *buf)
+{
+ int nlocal = atom->nlocal;
+ if (nlocal == nmax) grow(0);
+
+ int m = 1;
+ x[nlocal][0] = buf[m++];
+ x[nlocal][1] = buf[m++];
+ x[nlocal][2] = buf[m++];
+ v[nlocal][0] = buf[m++];
+ v[nlocal][1] = buf[m++];
+ v[nlocal][2] = buf[m++];
+ tag[nlocal] = (tagint) ubuf(buf[m++]).i;
+ type[nlocal] = (int) ubuf(buf[m++]).i;
+ mask[nlocal] = (int) ubuf(buf[m++]).i;
+ image[nlocal] = (imageint) ubuf(buf[m++]).i;
+ dpdTheta[nlocal] = buf[m++];
+ uCond[nlocal] = buf[m++];
+ uMech[nlocal] = buf[m++];
+
+ if (atom->nextra_grow)
+ for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
+ m += modify->fix[atom->extra_grow[iextra]]->
+ unpack_exchange(nlocal,&buf[m]);
+
+ atom->nlocal++;
+ return m;
+}
+
+/* ----------------------------------------------------------------------
+ size of restart data for all atoms owned by this proc
+ include extra data stored by fixes
+------------------------------------------------------------------------- */
+
+int AtomVecDPD::size_restart()
+{
+ int i;
+
+ int nlocal = atom->nlocal;
+ int n = 14 * nlocal; // 11 + dpdTheta + uCond + uMech
+
+ if (atom->nextra_restart)
+ for (int iextra = 0; iextra < atom->nextra_restart; iextra++)
+ for (i = 0; i < nlocal; i++)
+ n += modify->fix[atom->extra_restart[iextra]]->size_restart(i);
+
+ return n;
+}
+
+/* ----------------------------------------------------------------------
+ pack atom I's data for restart file including extra quantities
+ xyz must be 1st 3 values, so that read_restart can test on them
+ molecular types may be negative, but write as positive
+------------------------------------------------------------------------- */
+
+int AtomVecDPD::pack_restart(int i, double *buf)
+{
+ int m = 1;
+ buf[m++] = x[i][0];
+ buf[m++] = x[i][1];
+ buf[m++] = x[i][2];
+ buf[m++] = ubuf(tag[i]).d;
+ buf[m++] = ubuf(type[i]).d;
+ buf[m++] = ubuf(mask[i]).d;
+ buf[m++] = ubuf(image[i]).d;
+ buf[m++] = v[i][0];
+ buf[m++] = v[i][1];
+ buf[m++] = v[i][2];
+ buf[m++] = dpdTheta[i];
+ buf[m++] = uCond[i];
+ buf[m++] = uMech[i];
+
+ if (atom->nextra_restart)
+ for (int iextra = 0; iextra < atom->nextra_restart; iextra++)
+ m += modify->fix[atom->extra_restart[iextra]]->pack_restart(i,&buf[m]);
+
+ buf[0] = m;
+ return m;
+}
+
+/* ----------------------------------------------------------------------
+ unpack data for one atom from restart file including extra quantities
+------------------------------------------------------------------------- */
+
+int AtomVecDPD::unpack_restart(double *buf)
+{
+ int nlocal = atom->nlocal;
+ if (nlocal == nmax) {
+ grow(0);
+ if (atom->nextra_store)
+ memory->grow(atom->extra,nmax,atom->nextra_store,"atom:extra");
+ }
+
+ int m = 1;
+ x[nlocal][0] = buf[m++];
+ x[nlocal][1] = buf[m++];
+ x[nlocal][2] = buf[m++];
+ tag[nlocal] = (tagint) ubuf(buf[m++]).i;
+ type[nlocal] = (int) ubuf(buf[m++]).i;
+ mask[nlocal] = (int) ubuf(buf[m++]).i;
+ image[nlocal] = (imageint) ubuf(buf[m++]).i;
+ v[nlocal][0] = buf[m++];
+ v[nlocal][1] = buf[m++];
+ v[nlocal][2] = buf[m++];
+ dpdTheta[nlocal] = buf[m++];
+ uCond[nlocal] = buf[m++];
+ uMech[nlocal] = buf[m++];
+
+ double **extra = atom->extra;
+ if (atom->nextra_store) {
+ int size = static_cast<int> (buf[0]) - m;
+ for (int i = 0; i < size; i++) extra[nlocal][i] = buf[m++];
+ }
+
+ atom->nlocal++;
+ return m;
+}
+
+/* ----------------------------------------------------------------------
+ create one atom of itype at coord
+ set other values to defaults
+------------------------------------------------------------------------- */
+
+void AtomVecDPD::create_atom(int itype, double *coord)
+{
+ int nlocal = atom->nlocal;
+ if (nlocal == nmax) grow(0);
+
+ tag[nlocal] = 0;
+ type[nlocal] = itype;
+ x[nlocal][0] = coord[0];
+ x[nlocal][1] = coord[1];
+ x[nlocal][2] = coord[2];
+ mask[nlocal] = 1;
+ image[nlocal] = ((imageint) IMGMAX << IMG2BITS) |
+ ((imageint) IMGMAX << IMGBITS) | IMGMAX;
+ v[nlocal][0] = 0.0;
+ v[nlocal][1] = 0.0;
+ v[nlocal][2] = 0.0;
+ dpdTheta[nlocal] = 0.0;
+ uCond[nlocal] = 0.0;
+ uMech[nlocal] = 0.0;
+ duCond[nlocal] = 0.0;
+ duMech[nlocal] = 0.0;
+
+ atom->nlocal++;
+}
+
+/* ----------------------------------------------------------------------
+ unpack one line from Atoms section of data file
+ initialize other atom quantities
+------------------------------------------------------------------------- */
+
+void AtomVecDPD::data_atom(double *coord, tagint imagetmp, char **values)
+{
+ int nlocal = atom->nlocal;
+ if (nlocal == nmax) grow(0);
+
+ tag[nlocal] = ATOTAGINT(values[0]);
+ type[nlocal] = atoi(values[1]);
+ if (type[nlocal] <= 0 || type[nlocal] > atom->ntypes)
+ error->one(FLERR,"Invalid atom type in Atoms section of data file");
+
+ dpdTheta[nlocal] = atof(values[2]);
+ if (dpdTheta[nlocal] <= 0)
+ error->one(FLERR,"Internal temperature in Atoms section of date file must be > zero");
+
+ x[nlocal][0] = coord[0];
+ x[nlocal][1] = coord[1];
+ x[nlocal][2] = coord[2];
+
+ image[nlocal] = imagetmp;
+
+ mask[nlocal] = 1;
+ v[nlocal][0] = 0.0;
+ v[nlocal][1] = 0.0;
+ v[nlocal][2] = 0.0;
+
+ uCond[nlocal] = 0.0;
+ uMech[nlocal] = 0.0;
+
+ atom->nlocal++;
+}
+
+/* ----------------------------------------------------------------------
+ unpack hybrid quantities from one line in Atoms section of data file
+ initialize other atom quantities for this sub-style
+------------------------------------------------------------------------- */
+
+int AtomVecDPD::data_atom_hybrid(int nlocal, char **values)
+{
+ dpdTheta[nlocal] = atof(values[0]);
+
+ return 1;
+}
+
+/* ----------------------------------------------------------------------
+ pack atom info for data file including 3 image flags
+------------------------------------------------------------------------- */
+
+void AtomVecDPD::pack_data(double **buf)
+{
+ int nlocal = atom->nlocal;
+ for (int i = 0; i < nlocal; i++) {
+ buf[i][0] = ubuf(tag[i]).d;
+ buf[i][1] = ubuf(type[i]).d;
+ buf[i][2] = dpdTheta[i];
+ buf[i][3] = x[i][0];
+ buf[i][4] = x[i][1];
+ buf[i][5] = x[i][2];
+ buf[i][6] = ubuf((image[i] & IMGMASK) - IMGMAX).d;
+ buf[i][7] = ubuf((image[i] >> IMGBITS & IMGMASK) - IMGMAX).d;
+ buf[i][8] = ubuf((image[i] >> IMG2BITS) - IMGMAX).d;
+ }
+}
+
+/* ----------------------------------------------------------------------
+ pack hybrid atom info for data file
+------------------------------------------------------------------------- */
+
+int AtomVecDPD::pack_data_hybrid(int i, double *buf)
+{
+ buf[0] = dpdTheta[i];
+ return 1;
+}
+
+/* ----------------------------------------------------------------------
+ write atom info to data file including 3 image flags
+------------------------------------------------------------------------- */
+
+void AtomVecDPD::write_data(FILE *fp, int n, double **buf)
+{
+ for (int i = 0; i < n; i++)
+ fprintf(fp,TAGINT_FORMAT " %d %-1.16e %-1.16e %-1.16e %-1.16e %d %d %d\n",
+ (tagint) ubuf(buf[i][0]).i,(int) ubuf(buf[i][1]).i,
+ buf[i][2],buf[i][3],buf[i][4],buf[i][5],
+ (int) ubuf(buf[i][6]).i,(int) ubuf(buf[i][7]).i,
+ (int) ubuf(buf[i][8]).i);
+}
+
+/* ----------------------------------------------------------------------
+ write hybrid atom info to data file
+------------------------------------------------------------------------- */
+
+int AtomVecDPD::write_data_hybrid(FILE *fp, double *buf)
+{
+ fprintf(fp," %-1.16e",buf[0]);
+ return 1;
+}
+
+/* ----------------------------------------------------------------------
+ return # of bytes of allocated memory
+------------------------------------------------------------------------- */
+
+bigint AtomVecDPD::memory_usage()
+{
+ bigint bytes = 0;
+
+ if (atom->memcheck("tag")) bytes += memory->usage(tag,nmax);
+ if (atom->memcheck("type")) bytes += memory->usage(type,nmax);
+ if (atom->memcheck("mask")) bytes += memory->usage(mask,nmax);
+ if (atom->memcheck("image")) bytes += memory->usage(image,nmax);
+ if (atom->memcheck("x")) bytes += memory->usage(x,nmax,3);
+ if (atom->memcheck("v")) bytes += memory->usage(v,nmax,3);
+ if (atom->memcheck("f")) bytes += memory->usage(f,nmax*comm->nthreads,3);
+ if (atom->memcheck("dpdTheta")) bytes += memory->usage(dpdTheta,nmax);
+ if (atom->memcheck("uCond")) bytes += memory->usage(uCond,nmax);
+ if (atom->memcheck("uMech")) bytes += memory->usage(uMech,nmax);
+ if (atom->memcheck("duCond")) bytes += memory->usage(duCond,nmax);
+ if (atom->memcheck("duMech")) bytes += memory->usage(duMech,nmax);
+
+ return bytes;
+}
diff --git a/src/USER-DPD/atom_vec_dpd.h b/src/USER-DPD/atom_vec_dpd.h
new file mode 100644
index 0000000000000000000000000000000000000000..0f5c348f5cb71d06c5a1d05aba40af867aa1e6a2
--- /dev/null
+++ b/src/USER-DPD/atom_vec_dpd.h
@@ -0,0 +1,95 @@
+/* ----------------------------------------------------------------------
+ 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: James Larentzos (Engility Corporation)
+------------------------------------------------------------------------- */
+
+#ifdef ATOM_CLASS
+
+AtomStyle(dpd,AtomVecDPD)
+
+#else
+
+#ifndef LMP_ATOM_VEC_DPD_H
+#define LMP_ATOM_VEC_DPD_H
+
+#include "atom_vec.h"
+
+namespace LAMMPS_NS {
+
+class AtomVecDPD : public AtomVec {
+ public:
+ AtomVecDPD(class LAMMPS *);
+ virtual ~AtomVecDPD() {}
+ void grow(int);
+ void grow_reset();
+ void copy(int, int, int);
+ virtual int pack_comm(int, int *, double *, int, int *);
+ virtual int pack_comm_vel(int, int *, double *, int, int *);
+ int pack_comm_hybrid(int, int *, double *);
+ virtual void unpack_comm(int, int, double *);
+ virtual void unpack_comm_vel(int, int, double *);
+ int unpack_comm_hybrid(int, int, double *);
+ int pack_reverse(int, int, double *);
+ void unpack_reverse(int, int *, double *);
+ virtual int pack_border(int, int *, double *, int, int *);
+ virtual int pack_border_vel(int, int *, double *, int, int *);
+ int pack_border_hybrid(int, int *, double *);
+ virtual void unpack_border(int, int, double *);
+ virtual void unpack_border_vel(int, int, double *);
+ int unpack_border_hybrid(int, int, double *);
+ virtual int pack_exchange(int, double *);
+ virtual int unpack_exchange(double *);
+ int size_restart();
+ int pack_restart(int, double *);
+ int unpack_restart(double *);
+ void create_atom(int, double *);
+ void data_atom(double *, imageint, char **);
+ int data_atom_hybrid(int, char **);
+ void pack_data(double **);
+ int pack_data_hybrid(int, double *);
+ void write_data(FILE *, int, double **);
+ int write_data_hybrid(FILE *, double *);
+ bigint memory_usage();
+ double *uCond,*uMech,*dpdTheta;
+ double *duCond,*duMech;
+
+ protected:
+ tagint *tag;
+ int *type,*mask;
+ imageint *image;
+ double **x,**v,**f;
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Per-processor system is too big
+
+The number of owned atoms plus ghost atoms on a single
+processor must fit in 32-bit integer.
+
+E: Invalid atom type in Atoms section of data file
+
+Atom types must range from 1 to specified # of types.
+
+E: Internal temperature in Atoms section of data file must be > zero
+
+All internal temperatures must be > zero
+
+*/
diff --git a/src/USER-DPD/compute_dpd.cpp b/src/USER-DPD/compute_dpd.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a1299ca28a6bdd857645b2bca93fd209e29a763c
--- /dev/null
+++ b/src/USER-DPD/compute_dpd.cpp
@@ -0,0 +1,88 @@
+/* ----------------------------------------------------------------------
+ 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: James Larentzos (Engility Corporation)
+------------------------------------------------------------------------- */
+
+#include <mpi.h>
+#include "compute_dpd.h"
+#include "atom.h"
+#include "update.h"
+#include "force.h"
+#include "domain.h"
+#include "group.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+ComputeDpd::ComputeDpd(LAMMPS *lmp, int narg, char **arg) :
+ Compute(lmp, narg, arg)
+{
+ if (narg != 3) error->all(FLERR,"Illegal compute dpd command");
+
+ vector_flag = 1;
+ size_vector = 5;
+ extvector = 0;
+
+ vector = new double[size_vector];
+
+ if (atom->dpd_flag != 1) error->all(FLERR,"compute dpd requires atom_style with internal temperature and energies (e.g. dpd)");
+}
+
+/* ---------------------------------------------------------------------- */
+
+ComputeDpd::~ComputeDpd()
+{
+
+ delete [] vector;
+
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeDpd::compute_vector()
+{
+ invoked_vector = update->ntimestep;
+
+ double *uCond = atom->uCond;
+ double *uMech = atom->uMech;
+ double *dpdTheta = atom->dpdTheta;
+ int nlocal = atom->nlocal;
+ int *mask = atom->mask;
+ int natoms;
+
+ dpdU = new double[size_vector];
+
+ for (int i = 0; i < size_vector; i++) dpdU[i] = double(0.0);
+
+ for (int i = 0; i < nlocal; i++){
+ if (mask[i] & groupbit){
+ dpdU[0] += uCond[i];
+ dpdU[1] += uMech[i];
+ dpdU[2] += uCond[i] + uMech[i];
+ dpdU[3] += 1.0 / dpdTheta[i];
+ dpdU[4]++;
+ }
+ }
+
+ MPI_Allreduce(dpdU,vector,size_vector,MPI_DOUBLE,MPI_SUM,world);
+
+ natoms = vector[4];
+ vector[3] = natoms / vector[3];
+
+ delete [] dpdU;
+
+}
diff --git a/src/USER-DPD/compute_dpd.h b/src/USER-DPD/compute_dpd.h
new file mode 100644
index 0000000000000000000000000000000000000000..022beb68b3542545c03a9721517b7d42a1acdce9
--- /dev/null
+++ b/src/USER-DPD/compute_dpd.h
@@ -0,0 +1,60 @@
+/* ----------------------------------------------------------------------
+ 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: James Larentzos (Engility Corporation)
+------------------------------------------------------------------------- */
+
+#ifdef COMPUTE_CLASS
+
+ComputeStyle(dpd,ComputeDpd)
+
+#else
+
+#ifndef LMP_COMPUTE_DPD_H
+#define LMP_COMPUTE_DPD_H
+
+#include "compute.h"
+
+namespace LAMMPS_NS {
+
+class ComputeDpd : public Compute {
+ public:
+ ComputeDpd(class LAMMPS *, int, char **);
+ ~ComputeDpd();
+ void init() {}
+ void compute_vector();
+
+ private:
+ double *dpdU;
+
+};
+
+}
+
+#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: compute dpd requires atom_style with internal temperature and energies (e.g. dpd)
+
+Self-explanatory.
+
+*/
diff --git a/src/USER-DPD/compute_dpd_atom.cpp b/src/USER-DPD/compute_dpd_atom.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..52f2abc68f44ca18822f5f3fb943f0c631da290c
--- /dev/null
+++ b/src/USER-DPD/compute_dpd_atom.cpp
@@ -0,0 +1,107 @@
+/* ----------------------------------------------------------------------
+ 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: James Larentzos (Engility Corporation)
+------------------------------------------------------------------------- */
+
+#include "math.h"
+#include <string.h>
+#include <stdlib.h>
+#include "compute_dpd_atom.h"
+#include "atom.h"
+#include "update.h"
+#include "modify.h"
+#include "domain.h"
+#include "group.h"
+#include "memory.h"
+#include "error.h"
+#include "comm.h"
+
+#include <vector>
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+ComputeDpdAtom::ComputeDpdAtom(LAMMPS *lmp, int narg, char **arg) :
+ Compute(lmp, narg, arg)
+{
+ if (narg != 3) error->all(FLERR,"Illegal compute dpd/atom command");
+
+ peratom_flag = 1;
+ size_peratom_cols = 3;
+
+ nmax = 0;
+ dpdAtom = NULL;
+
+ if (atom->dpd_flag != 1) error->all(FLERR,"compute dpd requires atom_style with internal temperature and energies (e.g. dpd)");
+}
+
+/* ---------------------------------------------------------------------- */
+
+ComputeDpdAtom::~ComputeDpdAtom()
+{
+ memory->destroy(dpdAtom);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeDpdAtom::init()
+{
+ int count = 0;
+ for (int i = 0; i < modify->ncompute; i++)
+ if (strcmp(modify->compute[i]->style,"dpd/atom") == 0) count++;
+ if (count > 1 && comm->me == 0)
+ error->warning(FLERR,"More than one compute dpd/atom command");
+}
+
+/* ----------------------------------------------------------------------
+ gather compute vector data from other nodes
+------------------------------------------------------------------------- */
+
+void ComputeDpdAtom::compute_peratom()
+{
+
+ invoked_peratom = update->ntimestep;
+
+ double *uCond = atom->uCond;
+ double *uMech = atom->uMech;
+ double *dpdTheta = atom->dpdTheta;
+ int nlocal = atom->nlocal;
+ int *mask = atom->mask;
+ if (nlocal > nmax) {
+ memory->destroy(dpdAtom);
+ nmax = atom->nmax;
+ memory->create(dpdAtom,nmax,size_peratom_cols,"dpd/atom:dpdAtom");
+ array_atom = dpdAtom;
+ }
+
+ for (int i = 0; i < nlocal; i++){
+ if (mask[i] & groupbit){
+ dpdAtom[i][0] = uCond[i];
+ dpdAtom[i][1] = uMech[i];
+ dpdAtom[i][2] = dpdTheta[i];
+ }
+ }
+}
+
+/* ----------------------------------------------------------------------
+ memory usage of local atom-based array
+------------------------------------------------------------------------- */
+
+double ComputeDpdAtom::memory_usage()
+{
+ double bytes = size_peratom_cols * nmax * sizeof(double);
+ return bytes;
+}
diff --git a/src/USER-DPD/compute_dpd_atom.h b/src/USER-DPD/compute_dpd_atom.h
new file mode 100644
index 0000000000000000000000000000000000000000..771b01510eddcd526ba3afaebdd0d2e6c3cc0aa0
--- /dev/null
+++ b/src/USER-DPD/compute_dpd_atom.h
@@ -0,0 +1,65 @@
+/* ----------------------------------------------------------------------
+ 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: James Larentzos (Engility Corporation)
+------------------------------------------------------------------------- */
+
+#ifdef COMPUTE_CLASS
+
+ComputeStyle(dpd/atom,ComputeDpdAtom)
+
+#else
+
+#ifndef LMP_COMPUTE_DPD_ATOM_H
+#define LMP_COMPUTE_DPD_ATOM_H
+
+#include "compute.h"
+
+namespace LAMMPS_NS {
+
+class ComputeDpdAtom : public Compute {
+ public:
+ ComputeDpdAtom(class LAMMPS *, int, char **);
+ ~ComputeDpdAtom();
+ void init();
+ void compute_peratom();
+ double memory_usage();
+
+ private:
+ int nmax;
+ double **dpdAtom;
+};
+
+}
+
+#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: compute dpd requires atom_style with internal temperature and energies (e.g. dpd)
+
+Self-explanatory
+
+W: More than one compute dpd/atom command
+
+Self-explanatory
+
+*/
diff --git a/src/USER-DPD/fix_eos_cv.cpp b/src/USER-DPD/fix_eos_cv.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..7207a80a95425c621af4f455b87373963cad3d31
--- /dev/null
+++ b/src/USER-DPD/fix_eos_cv.cpp
@@ -0,0 +1,108 @@
+/* ----------------------------------------------------------------------
+ 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: James Larentzos (Engility Corporation)
+------------------------------------------------------------------------- */
+
+#include <stdlib.h>
+#include <string.h>
+#include "fix_eos_cv.h"
+#include "atom.h"
+#include "error.h"
+#include "force.h"
+
+using namespace LAMMPS_NS;
+using namespace FixConst;
+
+/* ---------------------------------------------------------------------- */
+
+FixEOScv::FixEOScv(LAMMPS *lmp, int narg, char **arg) :
+ Fix(lmp, narg, arg)
+{
+ if (narg != 4) error->all(FLERR,"Illegal fix eos/cv command");
+ cvEOS = force->numeric(FLERR,arg[3]);
+ if(cvEOS <= double(0.0)) error->all(FLERR,"EOS cv must be > 0.0");
+
+ restart_peratom = 1;
+ nevery = 1;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixEOScv::setmask()
+{
+ int mask = 0;
+ mask |= POST_INTEGRATE;
+ mask |= END_OF_STEP;
+ return mask;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixEOScv::init()
+{
+ int nlocal = atom->nlocal;
+ int *mask = atom->mask;
+ double *uCond = atom->uCond;
+ double *uMech = atom->uMech;
+ double *dpdTheta = atom->dpdTheta;
+
+ if(this->restart_reset){
+ for (int i = 0; i < nlocal; i++)
+ if (mask[i] & groupbit)
+ dpdTheta[i] = (uCond[i]+uMech[i])/cvEOS;
+ } else {
+ for (int i = 0; i < nlocal; i++)
+ if (mask[i] & groupbit) {
+ uCond[i] = double(0.5)*cvEOS*dpdTheta[i];
+ uMech[i] = double(0.5)*cvEOS*dpdTheta[i];
+ }
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixEOScv::post_integrate()
+{
+ int nlocal = atom->nlocal;
+ int *mask = atom->mask;
+ double *uCond = atom->uCond;
+ double *uMech = atom->uMech;
+ double *dpdTheta = atom->dpdTheta;
+
+ for (int i = 0; i < nlocal; i++)
+ if (mask[i] & groupbit){
+ dpdTheta[i] = (uCond[i]+uMech[i])/cvEOS;
+ if(dpdTheta[i] <= double(0.0))
+ error->one(FLERR,"Internal temperature < zero");
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixEOScv::end_of_step()
+{
+ int nlocal = atom->nlocal;
+ int *mask = atom->mask;
+ double *uCond = atom->uCond;
+ double *uMech = atom->uMech;
+ double *dpdTheta = atom->dpdTheta;
+
+ for (int i = 0; i < nlocal; i++)
+ if (mask[i] & groupbit){
+ dpdTheta[i] = (uCond[i]+uMech[i])/cvEOS;
+ if(dpdTheta[i] <= double(0.0))
+ error->one(FLERR,"Internal temperature < zero");
+ }
+}
diff --git a/src/USER-DPD/fix_eos_cv.h b/src/USER-DPD/fix_eos_cv.h
new file mode 100644
index 0000000000000000000000000000000000000000..ab59275ec8c3ff31e236abd97d6a0b165cb79da0
--- /dev/null
+++ b/src/USER-DPD/fix_eos_cv.h
@@ -0,0 +1,65 @@
+/* ----------------------------------------------------------------------
+ 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: James Larentzos (Engility Corporation)
+------------------------------------------------------------------------- */
+
+#ifdef FIX_CLASS
+
+FixStyle(eos/cv,FixEOScv)
+
+#else
+
+#ifndef LMP_FIX_EOS_CV_H
+#define LMP_FIX_EOS_CV_H
+
+#include "fix.h"
+
+namespace LAMMPS_NS {
+
+class FixEOScv : public Fix {
+ public:
+ FixEOScv(class LAMMPS *, int, char **);
+ virtual ~FixEOScv() {}
+ int setmask();
+ virtual void init();
+ virtual void post_integrate();
+ virtual void end_of_step();
+
+ protected:
+ double cvEOS;
+};
+
+}
+
+#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: EOS cv must be > 0.0
+
+The constant volume heat capacity must be larger than zero.
+
+E: Internal temperature < zero
+
+Self-explanatory. EOS may not be valid under current simulation conditions.
+
+*/
diff --git a/src/USER-DPD/fix_eos_table.cpp b/src/USER-DPD/fix_eos_table.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..85aa20ee981d10e308806f2838be31e7afc5d1d6
--- /dev/null
+++ b/src/USER-DPD/fix_eos_table.cpp
@@ -0,0 +1,433 @@
+/* ----------------------------------------------------------------------
+ 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: James Larentzos (Engility Corporation)
+------------------------------------------------------------------------- */
+
+#include <stdlib.h>
+#include <string.h>
+#include "fix_eos_table.h"
+#include "atom.h"
+#include "error.h"
+#include "force.h"
+#include "memory.h"
+
+#define MAXLINE 1024
+
+using namespace LAMMPS_NS;
+using namespace FixConst;
+
+/* ---------------------------------------------------------------------- */
+
+FixEOStable::FixEOStable(LAMMPS *lmp, int narg, char **arg) :
+ Fix(lmp, narg, arg)
+{
+ if (narg != 7) error->all(FLERR,"Illegal fix eos/table command");
+ restart_peratom = 1;
+ nevery = 1;
+
+ if (strcmp(arg[3],"linear") == 0) tabstyle = LINEAR;
+ else error->all(FLERR,"Unknown table style in fix eos/table");
+
+ tablength = force->inumeric(FLERR,arg[5]);
+ if (tablength < 2) error->all(FLERR,"Illegal number of eos/table entries");
+
+ ntables = 0;
+ tables = NULL;
+ int me;
+ MPI_Comm_rank(world,&me);
+ tables = (Table *)
+ memory->srealloc(tables,(ntables+2)*sizeof(Table),"eos:tables");
+ Table *tb = &tables[ntables];
+ Table *tb2 = &tables[ntables+1];
+ null_table(tb);
+ null_table(tb2);
+ if (me == 0) read_table(tb,tb2,arg[4],arg[6]);
+ bcast_table(tb);
+ bcast_table(tb2);
+
+ // error check on table parameters
+
+ if (tb->ninput <= 1) error->one(FLERR,"Invalid eos/table length");
+
+ tb->lo = tb->rfile[0];
+ tb->hi = tb->rfile[tb->ninput-1];
+ if (tb->lo >= tb->hi) error->all(FLERR,"eos/table values are not increasing");
+
+ if (tb2->ninput <= 1) error->one(FLERR,"Invalid eos/table length");
+
+ tb2->lo = tb2->rfile[0];
+ tb2->hi = tb2->rfile[tb2->ninput-1];
+ if (tb2->lo >= tb2->hi) error->all(FLERR,"eos/table values are not increasing");
+
+ // spline read-in and compute r,e,f vectors within table
+
+ spline_table(tb);
+ compute_table(tb);
+ spline_table(tb2);
+ compute_table(tb2);
+ ntables++;
+}
+
+/* ---------------------------------------------------------------------- */
+
+FixEOStable::~FixEOStable()
+{
+ for (int m = 0; m < 2*ntables; m++) free_table(&tables[m]);
+ memory->sfree(tables);
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixEOStable::setmask()
+{
+ int mask = 0;
+ mask |= POST_INTEGRATE;
+ mask |= END_OF_STEP;
+ return mask;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixEOStable::init()
+{
+ int nlocal = atom->nlocal;
+ int *mask = atom->mask;
+ double *uCond = atom->uCond;
+ double *uMech = atom->uMech;
+ double *dpdTheta = atom->dpdTheta;
+ double tmp;
+
+ if(this->restart_reset){
+ for (int i = 0; i < nlocal; i++)
+ if (mask[i] & groupbit)
+ temperature_lookup(uCond[i]+uMech[i],dpdTheta[i]);
+ } else {
+ for (int i = 0; i < nlocal; i++)
+ if (mask[i] & groupbit) {
+ energy_lookup(dpdTheta[i],tmp);
+ uCond[i] = tmp / double(2.0);
+ uMech[i] = tmp / double(2.0);
+ }
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixEOStable::post_integrate()
+{
+ int nlocal = atom->nlocal;
+ int *mask = atom->mask;
+ double *uCond = atom->uCond;
+ double *uMech = atom->uMech;
+ double *dpdTheta = atom->dpdTheta;
+
+ for (int i = 0; i < nlocal; i++)
+ if (mask[i] & groupbit){
+ temperature_lookup(uCond[i]+uMech[i],dpdTheta[i]);
+ if(dpdTheta[i] <= double(0.0))
+ error->one(FLERR,"Internal temperature < zero");
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixEOStable::end_of_step()
+{
+ int nlocal = atom->nlocal;
+ int *mask = atom->mask;
+ double *uCond = atom->uCond;
+ double *uMech = atom->uMech;
+ double *dpdTheta = atom->dpdTheta;
+
+ for (int i = 0; i < nlocal; i++)
+ if (mask[i] & groupbit){
+ temperature_lookup(uCond[i]+uMech[i],dpdTheta[i]);
+ if(dpdTheta[i] <= double(0.0))
+ error->one(FLERR,"Internal temperature < zero");
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixEOStable::null_table(Table *tb)
+{
+ tb->rfile = tb->efile = NULL;
+ tb->e2file = NULL;
+ tb->r = tb->e = tb->de = NULL;
+ tb->e2 = NULL;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixEOStable::free_table(Table *tb)
+{
+ memory->destroy(tb->rfile);
+ memory->destroy(tb->efile);
+ memory->destroy(tb->e2file);
+
+ memory->destroy(tb->r);
+ memory->destroy(tb->e);
+ memory->destroy(tb->de);
+ memory->destroy(tb->e2);
+}
+
+/* ----------------------------------------------------------------------
+ read table file, only called by proc 0
+------------------------------------------------------------------------- */
+
+void FixEOStable::read_table(Table *tb, Table *tb2, char *file, char *keyword)
+{
+ char line[MAXLINE];
+
+ // open file
+
+ FILE *fp = fopen(file,"r");
+ if (fp == NULL) {
+ char str[128];
+ sprintf(str,"Cannot open file %s",file);
+ error->one(FLERR,str);
+ }
+
+ // loop until section found with matching keyword
+
+ while (1) {
+ if (fgets(line,MAXLINE,fp) == NULL)
+ error->one(FLERR,"Did not find keyword in table file");
+ if (strspn(line," \t\n\r") == strlen(line)) continue; // blank line
+ if (line[0] == '#') continue; // comment
+ char *word = strtok(line," \t\n\r");
+ if (strcmp(word,keyword) == 0) break; // matching keyword
+ fgets(line,MAXLINE,fp); // no match, skip section
+ param_extract(tb,tb2,line);
+ fgets(line,MAXLINE,fp);
+ for (int i = 0; i < tb->ninput; i++) fgets(line,MAXLINE,fp);
+ }
+
+ // read args on 2nd line of section
+ // allocate table arrays for file values
+
+ fgets(line,MAXLINE,fp);
+ param_extract(tb,tb2,line);
+ memory->create(tb->rfile,tb->ninput,"eos:rfile");
+ memory->create(tb->efile,tb->ninput,"eos:efile");
+ memory->create(tb2->rfile,tb2->ninput,"eos:rfile2");
+ memory->create(tb2->efile,tb2->ninput,"eos:efile2");
+
+ // read r,e table values from file
+
+ int itmp;
+ fgets(line,MAXLINE,fp);
+ for (int i = 0; i < tb->ninput; i++) {
+ fgets(line,MAXLINE,fp);
+ sscanf(line,"%d %lg %lg",&itmp,&tb->rfile[i],&tb->efile[i]);
+ sscanf(line,"%d %lg %lg",&itmp,&tb2->efile[i],&tb2->rfile[i]);
+ }
+
+ fclose(fp);
+}
+
+/* ----------------------------------------------------------------------
+ build spline representation of e,f over entire range of read-in table
+ this function sets these values in e2file
+------------------------------------------------------------------------- */
+
+void FixEOStable::spline_table(Table *tb)
+{
+ memory->create(tb->e2file,tb->ninput,"eos:e2file");
+}
+
+/* ----------------------------------------------------------------------
+ compute r,e,f vectors from splined values
+------------------------------------------------------------------------- */
+
+void FixEOStable::compute_table(Table *tb)
+{
+ // delta = table spacing for N-1 bins
+ int tlm1 = tablength-1;
+
+ tb->delta = (tb->hi - tb->lo)/ tlm1;
+ tb->invdelta = 1.0/tb->delta;
+ tb->deltasq6 = tb->delta*tb->delta / 6.0;
+
+ // N-1 evenly spaced bins in r from min to max
+ // r,e = value at lower edge of bin
+ // de values = delta values of e,f
+ // r,e are N in length so de arrays can compute difference
+
+ memory->create(tb->r,tablength,"eos:r");
+ memory->create(tb->e,tablength,"eos:e");
+ memory->create(tb->de,tlm1,"eos:de");
+ memory->create(tb->e2,tablength,"eos:e2");
+
+ double a;
+ for (int i = 0; i < tablength; i++) {
+ a = tb->lo + i*tb->delta;
+ tb->r[i] = a;
+ tb->e[i] = splint(tb->rfile,tb->efile,tb->e2file,tb->ninput,a);
+ }
+
+ for (int i = 0; i < tlm1; i++) {
+ tb->de[i] = tb->e[i+1] - tb->e[i];
+ }
+}
+
+/* ----------------------------------------------------------------------
+ extract attributes from parameter line in table section
+ format of line: N value
+ N is required, other params are optional
+------------------------------------------------------------------------- */
+
+void FixEOStable::param_extract(Table *tb, Table *tb2, char *line)
+{
+ tb->ninput = 0;
+ tb2->ninput = 0;
+
+ char *word = strtok(line," \t\n\r\f");
+ while (word) {
+ if (strcmp(word,"N") == 0) {
+ word = strtok(NULL," \t\n\r\f");
+ tb->ninput = atoi(word);
+ tb2->ninput = atoi(word);
+ } else {
+ error->one(FLERR,"Invalid keyword in fix eos/table parameters");
+ }
+ word = strtok(NULL," \t\n\r\f");
+ }
+
+ if (tb->ninput == 0) error->one(FLERR,"fix eos/table parameters did not set N");
+ if (tb2->ninput == 0) error->one(FLERR,"fix eos/table parameters did not set N");
+}
+
+/* ----------------------------------------------------------------------
+ broadcast read-in table info from proc 0 to other procs
+ this function communicates these values in Table:
+ ninput,rfile,efile
+------------------------------------------------------------------------- */
+
+void FixEOStable::bcast_table(Table *tb)
+{
+ MPI_Bcast(&tb->ninput,1,MPI_INT,0,world);
+
+ int me;
+ MPI_Comm_rank(world,&me);
+ if (me > 0) {
+ memory->create(tb->rfile,tb->ninput,"eos:rfile");
+ memory->create(tb->efile,tb->ninput,"eos:efile");
+ }
+
+ MPI_Bcast(tb->rfile,tb->ninput,MPI_DOUBLE,0,world);
+ MPI_Bcast(tb->efile,tb->ninput,MPI_DOUBLE,0,world);
+}
+
+/* ----------------------------------------------------------------------
+ spline and splint routines modified from Numerical Recipes
+------------------------------------------------------------------------- */
+
+void FixEOStable::spline(double *x, double *y, int n,
+ double yp1, double ypn, double *y2)
+{
+ int i,k;
+ double p,qn,sig,un;
+ double *u = new double[n];
+
+ if (yp1 > 0.99e30) y2[0] = u[0] = 0.0;
+ else {
+ y2[0] = -0.5;
+ u[0] = (3.0/(x[1]-x[0])) * ((y[1]-y[0]) / (x[1]-x[0]) - yp1);
+ }
+ for (i = 1; i < n-1; i++) {
+ sig = (x[i]-x[i-1]) / (x[i+1]-x[i-1]);
+ p = sig*y2[i-1] + 2.0;
+ y2[i] = (sig-1.0) / p;
+ u[i] = (y[i+1]-y[i]) / (x[i+1]-x[i]) - (y[i]-y[i-1]) / (x[i]-x[i-1]);
+ u[i] = (6.0*u[i] / (x[i+1]-x[i-1]) - sig*u[i-1]) / p;
+ }
+ if (ypn > 0.99e30) qn = un = 0.0;
+ else {
+ qn = 0.5;
+ un = (3.0/(x[n-1]-x[n-2])) * (ypn - (y[n-1]-y[n-2]) / (x[n-1]-x[n-2]));
+ }
+ y2[n-1] = (un-qn*u[n-2]) / (qn*y2[n-2] + 1.0);
+ for (k = n-2; k >= 0; k--) y2[k] = y2[k]*y2[k+1] + u[k];
+
+ delete [] u;
+}
+
+/* ---------------------------------------------------------------------- */
+
+double FixEOStable::splint(double *xa, double *ya, double *y2a, int n, double x)
+{
+ int klo,khi,k;
+ double h,b,a,y;
+
+ klo = 0;
+ khi = n-1;
+ while (khi-klo > 1) {
+ k = (khi+klo) >> 1;
+ if (xa[k] > x) khi = k;
+ else klo = k;
+ }
+ h = xa[khi]-xa[klo];
+ a = (xa[khi]-x) / h;
+ b = (x-xa[klo]) / h;
+ y = a*ya[klo] + b*ya[khi] +
+ ((a*a*a-a)*y2a[klo] + (b*b*b-b)*y2a[khi]) * (h*h)/6.0;
+ return y;
+}
+
+/* ----------------------------------------------------------------------
+ calculate internal energy u at temperature t
+ insure t is between min/max
+------------------------------------------------------------------------- */
+
+void FixEOStable::energy_lookup(double t, double &u)
+{
+ int itable;
+ double fraction,a,b;
+
+ Table *tb = &tables[0];
+ if(t < tb->lo || t > tb->hi){
+ printf("Temperature=%lf TableMin=%lf TableMax=%lf\n",t,tb->lo,tb->hi);
+ error->one(FLERR,"Temperature is not within table cutoffs");
+ }
+
+ if (tabstyle == LINEAR) {
+ itable = static_cast<int> ((t - tb->lo) * tb->invdelta);
+ fraction = (t - tb->r[itable]) * tb->invdelta;
+ u = tb->e[itable] + fraction*tb->de[itable];
+ }
+}
+/* ----------------------------------------------------------------------
+ calculate temperature t at energy u
+ insure u is between min/max
+------------------------------------------------------------------------- */
+
+void FixEOStable::temperature_lookup(double u, double &t)
+{
+ int itable;
+ double fraction,a,b;
+
+ Table *tb = &tables[1];
+ if(u < tb->lo || u > tb->hi){
+ printf("Energy=%lf TableMin=%lf TableMax=%lf\n",u,tb->lo,tb->hi);
+ error->one(FLERR,"Energy is not within table cutoffs");
+ }
+
+ if (tabstyle == LINEAR) {
+ itable = static_cast<int> ((u - tb->lo) * tb->invdelta);
+ fraction = (u - tb->r[itable]) * tb->invdelta;
+ t = tb->e[itable] + fraction*tb->de[itable];
+ }
+}
diff --git a/src/USER-DPD/fix_eos_table.h b/src/USER-DPD/fix_eos_table.h
new file mode 100644
index 0000000000000000000000000000000000000000..b508ce3f39c6f0d5eb15c1b6d5da3ab28c52bb0b
--- /dev/null
+++ b/src/USER-DPD/fix_eos_table.h
@@ -0,0 +1,130 @@
+/* ----------------------------------------------------------------------
+ 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: James Larentzos (Engility Corporation)
+------------------------------------------------------------------------- */
+
+#ifdef FIX_CLASS
+
+FixStyle(eos/table,FixEOStable)
+
+#else
+
+#ifndef LMP_FIX_EOS_TABLE_H
+#define LMP_FIX_EOS_TABLE_H
+
+#include "fix.h"
+
+namespace LAMMPS_NS {
+
+class FixEOStable : public Fix {
+ public:
+ FixEOStable(class LAMMPS *, int, char **);
+ virtual ~FixEOStable();
+ int setmask();
+ virtual void init();
+ virtual void post_integrate();
+ virtual void end_of_step();
+ void energy_lookup(double, double &);
+ void temperature_lookup(double, double &);
+
+ protected:
+ enum{LINEAR};
+
+ int tabstyle,tablength;
+ struct Table {
+ int ninput;
+ double lo,hi;
+ double *rfile,*efile;
+ double *e2file;
+ double delta,invdelta,deltasq6;
+ double *r,*e,*de,*e2;
+ };
+ int ntables;
+ Table *tables;
+
+ void allocate();
+ void null_table(Table *);
+ void free_table(Table *);
+ void read_table(Table *, Table *, char *, char *);
+ void bcast_table(Table *);
+ void spline_table(Table *);
+ void compute_table(Table *);
+
+ void param_extract(Table *, Table *, char *);
+ void spline(double *, double *, int, double, double, double *);
+ double splint(double *, double *, double *, int, 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: Unknown table style in fix eos/table
+
+Style of table is invalid for use with fix eos/table command.
+
+E: Illegal number of eos/table entries
+
+There must be at least 2 table entries.
+
+E: Invalid eos/table length
+
+Length of read-in fix eos/table is invalid
+
+E: eos/table values are not increasing
+
+The EOS must be a monotonically, increasing function
+
+E: Internal temperature < zero
+
+Self-explanatory. EOS may not be valid under current simulation conditions.
+
+E: Cannot open file %s
+
+The specified file cannot be opened. Check that the path and name are
+correct.
+
+E: Did not find keyword in table file
+
+Keyword used in fix eos/table command was not found in table file.
+
+E: Invalid keyword in fix eos/table parameters
+
+Keyword used in list of table parameters is not recognized.
+
+E: fix eos/table parameters did not set N
+
+List of fix eos/table parameters must include N setting.
+
+E: Temperature is not within table cutoffs
+
+The internal temperature does not lie with the minimum
+and maximum temperature cutoffs of the table
+
+E: Energy is not within table cutoffs
+
+The internal energy does not lie with the minimum
+and maximum energy cutoffs of the table
+
+*/
diff --git a/src/USER-DPD/fix_shardlow.cpp b/src/USER-DPD/fix_shardlow.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..8dd6e4e2e709072372e4f78fd926fb8e0eb724b0
--- /dev/null
+++ b/src/USER-DPD/fix_shardlow.cpp
@@ -0,0 +1,568 @@
+/* ----------------------------------------------------------------------
+ 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 authors:
+ James Larentzos and Timothy I. Mattox (Engility Corporation)
+
+ Martin Lisal (Institute of Chemical Process Fundamentals
+ of the Czech Academy of Sciences and J. E. Purkinje University)
+
+ John Brennan, Joshua Moore and William Mattson (Army Research Lab)
+
+ Please cite the related publications:
+ J. P. Larentzos, J. K. Brennan, J. D. Moore, M. Lisal, W. D. Mattson,
+ "Parallel implementation of isothermal and isoenergetic Dissipative
+ Particle Dynamics using Shardlow-like splitting algorithms",
+ Computer Physics Communications, 2014, 185, pp 1987--1998.
+
+ M. Lisal, J. K. Brennan, J. Bonet Avalos, "Dissipative particle dynamics
+ at isothermal, isobaric, isoenergetic, and isoenthalpic conditions using
+ Shardlow-like splitting algorithms", Journal of Chemical Physics, 2011,
+ 135, 204105.
+------------------------------------------------------------------------- */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "fix_shardlow.h"
+#include "atom.h"
+#include "force.h"
+#include "update.h"
+#include "respa.h"
+#include "error.h"
+#include <math.h>
+#include "atom_vec.h"
+#include "comm.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "neigh_request.h"
+#include "random_mars.h"
+#include "memory.h"
+#include "domain.h"
+#include "modify.h"
+#include "pair_dpd_fdt.h"
+#include "pair_dpd_fdt_energy.h"
+#include "pair.h"
+#include "citeme.h"
+
+using namespace LAMMPS_NS;
+using namespace FixConst;
+
+#define EPSILON 1.0e-10
+
+static const char cite_fix_shardlow[] =
+ "fix shardlow command:\n\n"
+ "@Article{Larentzos14,\n"
+ " author = {J. P. Larentzos, J. K. Brennan, J. D. Moore, M. Lisal, W. D. Mattson},\n"
+ " title = {Parallel implementation of isothermal and isoenergetic Dissipative Particle Dynamics using Shardlow-like splitting algorithms},\n"
+ " journal = {Computer Physics Communications},\n"
+ " year = 2014,\n"
+ " volume = 185\n"
+ " pages = {1987--1998}\n"
+ "}\n\n"
+ "@Article{Lisal11,\n"
+ " author = {M. Lisal, J. K. Brennan, J. Bonet Avalos},\n"
+ " title = {Dissipative particle dynamics at isothermal, isobaric, isoenergetic, and isoenthalpic conditions using Shardlow-like splitting algorithms},\n"
+ " journal = {Journal of Chemical Physics},\n"
+ " year = 2011,\n"
+ " volume = 135\n"
+ " pages = {204105}\n"
+ "}\n\n";
+
+/* ---------------------------------------------------------------------- */
+
+FixShardlow::FixShardlow(LAMMPS *lmp, int narg, char **arg) :
+ Fix(lmp, narg, arg)
+{
+ if (lmp->citeme) lmp->citeme->add(cite_fix_shardlow);
+
+ if (narg != 3) error->all(FLERR,"Illegal fix shardlow command");
+
+ pairDPD = NULL;
+ pairDPDE = NULL;
+ pairDPD = (PairDPDfdt *) force->pair_match("dpd/fdt",1);
+ pairDPDE = (PairDPDfdtEnergy *) force->pair_match("dpd/fdt/energy",1);
+
+ if(pairDPDE){
+ comm_forward = 10;
+ comm_reverse = 5;
+ } else {
+ comm_forward = 6;
+ comm_reverse = 3;
+ }
+
+ if(pairDPD == NULL && pairDPDE == NULL)
+ error->all(FLERR,"Must use pair_style dpd/fdt or dpd/fdt/energy with fix shardlow");
+
+ for (int i = 0; i < modify->nfix; i++)
+ if (strcmp(modify->fix[i]->style,"nve") == 0 || strcmp(modify->fix[i]->style,"nph") == 0)
+ error->all(FLERR,"A deterministic integrator must be specified after fix shardlow in input file (e.g. fix nve or fix nph).");
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixShardlow::setmask()
+{
+ int mask = 0;
+ mask |= INITIAL_INTEGRATE;
+ mask |= PRE_NEIGHBOR;
+ return mask;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixShardlow::init_list(int id, NeighList *ptr)
+{
+ list = ptr;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixShardlow::setup(int vflag)
+{
+ for (int i = 0; i < modify->nfix; i++)
+ if (strcmp(modify->fix[i]->style,"nvt") == 0 || strcmp(modify->fix[i]->style,"npt") == 0)
+ error->all(FLERR,"Cannot use constant temperature integration routines with DPD.");
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixShardlow::setup_pre_force(int vflag)
+{
+ neighbor->build_one(list);
+}
+
+/* ----------------------------------------------------------------------
+ Perform the stochastic integration and Shardlow update
+ Allow for both per-type and per-atom mass
+
+ NOTE: only implemented for orthogonal boxes, not triclinic
+------------------------------------------------------------------------- */
+
+void FixShardlow::initial_integrate(int vflag)
+{
+ int i,j,ii,jj,inum,jnum,itype,jtype;
+ int *ilist,*jlist,*numneigh,**firstneigh;
+ double xtmp,ytmp,ztmp,delx,dely,delz;
+ double vxtmp,vytmp,vztmp,delvx,delvy,delvz;
+ double rsq,r,rinv;
+ double dot,wd,wr,randnum,factor_dpd,factor_dpd1;
+ double dpx,dpy,dpz;
+ double denom, mu_ij;
+
+ double **x = atom->x;
+ double **v = atom->v;
+ double *rmass = atom->rmass;
+ double *mass = atom->mass;
+ int *type = atom->type;
+ int nlocal = atom->nlocal;
+ int nghost = atom->nghost;
+ int nall = nlocal + nghost;
+ double *special_lj = force->special_lj;
+ int newton_pair = force->newton_pair;
+ double randPair;
+
+ double *uCond = atom->uCond;
+ double *uMech = atom->uMech;
+ double *duCond = atom->duCond;
+ double *duMech = atom->duMech;
+ double *dpdTheta = atom->dpdTheta;
+ double kappa_ij, alpha_ij, theta_ij, gamma_ij, sigma_ij, u_ij;
+ double vxi, vyi, vzi, vxj, vyj, vzj;
+ double vx0i, vy0i, vz0i, vx0j, vy0j, vz0j;
+ double dot1, dot2, dot3, dot4;
+ double mass_i, mass_j;
+ double massinv_i, massinv_j;
+ double cut, cut2;
+
+ const double dt = update->dt;
+ const double dtsqrt = sqrt(dt);
+
+ // NOTE: this logic is specific to orthogonal boxes, not triclinic
+
+ // Enforce the constraint that ghosts must be contained in the nearest sub-domains
+ double bbx = domain->subhi[0] - domain->sublo[0];
+ double bby = domain->subhi[1] - domain->sublo[1];
+ double bbz = domain->subhi[2] - domain->sublo[2];
+
+ double rcut = double(2.0)*neighbor->cutneighmax;
+
+ if (domain->triclinic)
+ error->all(FLERR,"Fix shardlow does not yet support triclinic geometries");
+
+ if(rcut >= bbx || rcut >= bby || rcut>= bbz )
+ error->all(FLERR,"Shardlow algorithm requires sub-domain length > 2*(rcut+skin). Either reduce the number of processors requested, or change the cutoff/skin\n");
+
+ // Allocate memory for the dvSSA arrays
+ dvSSA = new double*[nall];
+ for (ii = 0; ii < nall; ii++) {
+ dvSSA[ii] = new double[3];
+ }
+
+ // Zero the momenta
+ for (ii = 0; ii < nlocal; ii++) {
+ dvSSA[ii][0] = double(0.0);
+ dvSSA[ii][1] = double(0.0);
+ dvSSA[ii][2] = double(0.0);
+ if(pairDPDE){
+ duCond[ii] = double(0.0);
+ duMech[ii] = double(0.0);
+ }
+ }
+
+ // Communicate the updated momenta and velocities to all nodes
+ comm->forward_comm_fix(this);
+
+ // Define pointers to access the neighbor list
+ if(pairDPDE){
+ inum = pairDPDE->list->inum;
+ ilist = pairDPDE->list->ilist;
+ numneigh = pairDPDE->list->numneigh;
+ firstneigh = pairDPDE->list->firstneigh;
+ } else {
+ inum = pairDPD->list->inum;
+ ilist = pairDPD->list->ilist;
+ numneigh = pairDPD->list->numneigh;
+ firstneigh = pairDPD->list->firstneigh;
+ }
+
+ //Loop over all 14 directions (8 stages)
+ for (int idir = 1; idir <=8; idir++){
+
+ // Loop over neighbors of my atoms
+ for (ii = 0; ii < inum; ii++) {
+ i = ilist[ii];
+
+ xtmp = x[i][0];
+ ytmp = x[i][1];
+ ztmp = x[i][2];
+
+ itype = type[i];
+ jlist = firstneigh[i];
+ jnum = numneigh[i];
+
+ // Loop over Directional Neighbors only
+ for (jj = 0; jj < jnum; jj++) {
+ j = jlist[jj];
+ j &= NEIGHMASK;
+ if (neighbor->ssa_airnum[j] != idir) continue;
+ jtype = type[j];
+
+ delx = xtmp - x[j][0];
+ dely = ytmp - x[j][1];
+ delz = ztmp - x[j][2];
+ rsq = delx*delx + dely*dely + delz*delz;
+
+ if(pairDPDE){
+ cut2 = pairDPDE->cutsq[itype][jtype];
+ cut = pairDPDE->cut[itype][jtype];
+ } else {
+ cut2 = pairDPD->cutsq[itype][jtype];
+ cut = pairDPD->cut[itype][jtype];
+ }
+
+ // if (rsq < pairDPD->cutsq[itype][jtype]) {
+ if (rsq < cut2) {
+ r = sqrt(rsq);
+ if (r < EPSILON) continue; // r can be 0.0 in DPD systems
+ rinv = double(1.0)/r;
+
+ // Store the velocities from previous Shardlow step
+ vx0i = v[i][0] + dvSSA[i][0];
+ vy0i = v[i][1] + dvSSA[i][1];
+ vz0i = v[i][2] + dvSSA[i][2];
+
+ vx0j = v[j][0] + dvSSA[j][0];
+ vy0j = v[j][1] + dvSSA[j][1];
+ vz0j = v[j][2] + dvSSA[j][2];
+
+ // Compute the velocity difference between atom i and atom j
+ delvx = vx0i - vx0j;
+ delvy = vy0i - vy0j;
+ delvz = vz0i - vz0j;
+
+ dot = (delx*delvx + dely*delvy + delz*delvz);
+ // wr = double(1.0) - r/pairDPD->cut[itype][jtype];
+ wr = double(1.0) - r/cut;
+ wd = wr*wr;
+
+ if(pairDPDE){
+ // Compute the current temperature
+ theta_ij = double(0.5)*(double(1.0)/dpdTheta[i] + double(1.0)/dpdTheta[j]);
+ theta_ij = double(1.0)/theta_ij;
+ sigma_ij = pairDPDE->sigma[itype][jtype];
+ randnum = pairDPDE->random->gaussian();
+ } else {
+ theta_ij = pairDPD->temperature;
+ sigma_ij = pairDPD->sigma[itype][jtype];
+ randnum = pairDPD->random->gaussian();
+ }
+
+ gamma_ij = sigma_ij*sigma_ij / (2.0*force->boltz*theta_ij);
+ randPair = sigma_ij*wr*randnum*dtsqrt;
+
+ factor_dpd = -dt*gamma_ij*wd*dot*rinv;
+ factor_dpd += randPair;
+ factor_dpd *= double(0.5);
+
+ // Compute momentum change between t and t+dt
+ dpx = factor_dpd*delx*rinv;
+ dpy = factor_dpd*dely*rinv;
+ dpz = factor_dpd*delz*rinv;
+
+ if (rmass) {
+ mass_i = rmass[i];
+ mass_j = rmass[j];
+ } else {
+ mass_i = mass[itype];
+ mass_j = mass[jtype];
+ }
+ massinv_i = double(1.0) / mass_i;
+ massinv_j = double(1.0) / mass_j;
+
+ // Update the delta velocity on i
+ dvSSA[i][0] += dpx*force->ftm2v*massinv_i;
+ dvSSA[i][1] += dpy*force->ftm2v*massinv_i;
+ dvSSA[i][2] += dpz*force->ftm2v*massinv_i;
+
+ if (newton_pair || j < nlocal) {
+ // Update the delta velocity on j
+ dvSSA[j][0] -= dpx*force->ftm2v*massinv_j;
+ dvSSA[j][1] -= dpy*force->ftm2v*massinv_j;
+ dvSSA[j][2] -= dpz*force->ftm2v*massinv_j;
+ }
+
+ //ii. Compute the velocity diff
+ delvx = v[i][0] + dvSSA[i][0] - v[j][0] - dvSSA[j][0];
+ delvy = v[i][1] + dvSSA[i][1] - v[j][1] - dvSSA[j][1];
+ delvz = v[i][2] + dvSSA[i][2] - v[j][2] - dvSSA[j][2];
+
+ dot = delx*delvx + dely*delvy + delz*delvz;
+
+ //iii. Compute dpi again
+ mu_ij = massinv_i + massinv_j;
+ denom = double(1.0) + double(0.5)*mu_ij*gamma_ij*wd*dt*force->ftm2v;
+ factor_dpd = -double(0.5)*dt*gamma_ij*wd*force->ftm2v/denom;
+ factor_dpd1 = factor_dpd*(mu_ij*randPair);
+ factor_dpd1 += randPair;
+ factor_dpd1 *= double(0.5);
+
+ // Compute the momentum change between t and t+dt
+ dpx = (factor_dpd*dot*rinv/force->ftm2v + factor_dpd1)*delx*rinv;
+ dpy = (factor_dpd*dot*rinv/force->ftm2v + factor_dpd1)*dely*rinv;
+ dpz = (factor_dpd*dot*rinv/force->ftm2v + factor_dpd1)*delz*rinv;
+
+ //Update the velocity change on i
+ dvSSA[i][0] += dpx*force->ftm2v*massinv_i;
+ dvSSA[i][1] += dpy*force->ftm2v*massinv_i;
+ dvSSA[i][2] += dpz*force->ftm2v*massinv_i;
+
+ if (newton_pair || j < nlocal) {
+ //Update the velocity change on j
+ dvSSA[j][0] -= dpx*force->ftm2v*massinv_j;
+ dvSSA[j][1] -= dpy*force->ftm2v*massinv_j;
+ dvSSA[j][2] -= dpz*force->ftm2v*massinv_j;
+ }
+
+ if(pairDPDE){
+ // Compute uCond
+ randnum = pairDPDE->random->gaussian();
+ kappa_ij = pairDPDE->kappa[itype][jtype];
+ alpha_ij = sqrt(2.0*force->boltz*kappa_ij);
+ randPair = alpha_ij*wr*randnum*dtsqrt;
+
+ factor_dpd = kappa_ij*(double(1.0)/dpdTheta[i] - double(1.0)/dpdTheta[j])*wd*dt;
+ factor_dpd += randPair;
+
+ duCond[i] += factor_dpd;
+ if (newton_pair || j < nlocal) {
+ duCond[j] -= factor_dpd;
+ }
+
+ // Compute uMech
+ vxi = v[i][0] + dvSSA[i][0];
+ vyi = v[i][1] + dvSSA[i][1];
+ vzi = v[i][2] + dvSSA[i][2];
+
+ vxj = v[j][0] + dvSSA[j][0];
+ vyj = v[j][1] + dvSSA[j][1];
+ vzj = v[j][2] + dvSSA[j][2];
+
+ dot1 = vxi*vxi + vyi*vyi + vzi*vzi;
+ dot2 = vxj*vxj + vyj*vyj + vzj*vzj;
+ dot3 = vx0i*vx0i + vy0i*vy0i + vz0i*vz0i;
+ dot4 = vx0j*vx0j + vy0j*vy0j + vz0j*vz0j;
+
+ dot1 = dot1*mass_i;
+ dot2 = dot2*mass_j;
+ dot3 = dot3*mass_i;
+ dot4 = dot4*mass_j;
+
+ factor_dpd = double(0.25)*(dot1+dot2-dot3-dot4)/force->ftm2v;
+ duMech[i] -= factor_dpd;
+ if (newton_pair || j < nlocal) {
+ duMech[j] -= factor_dpd;
+ }
+ }
+ }
+ }
+ }
+
+ // Communicate the ghost delta velocities to the locally owned atoms
+ comm->reverse_comm_fix(this);
+
+ // Zero the dv
+ for (ii = 0; ii < nlocal; ii++) {
+ // Shardlow update
+ v[ii][0] += dvSSA[ii][0];
+ v[ii][1] += dvSSA[ii][1];
+ v[ii][2] += dvSSA[ii][2];
+ dvSSA[ii][0] = double(0.0);
+ dvSSA[ii][1] = double(0.0);
+ dvSSA[ii][2] = double(0.0);
+ if(pairDPDE){
+ uCond[ii] += duCond[ii];
+ uMech[ii] += duMech[ii];
+ duCond[ii] = double(0.0);
+ duMech[ii] = double(0.0);
+ }
+ }
+
+ // Communicate the updated momenta and velocities to all nodes
+ comm->forward_comm_fix(this);
+
+ } //End Loop over all directions For idir = Top, Top-Right, Right, Bottom-Right, Back
+
+ for (ii = 0; ii < nall; ii++) {
+ delete dvSSA[ii];
+ }
+ delete [] dvSSA;
+}
+
+/* ----------------------------------------------------------------------
+ * assign owned and ghost atoms their ssa active interaction region numbers
+------------------------------------------------------------------------- */
+
+void FixShardlow::setup_pre_neighbor()
+{
+ neighbor->assign_ssa_airnums();
+}
+
+void FixShardlow::pre_neighbor()
+{
+ neighbor->assign_ssa_airnums();
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixShardlow::pack_forward_comm(int n, int *list, double *buf, int pbc_flag, int *pbc)
+{
+ int ii,jj,m;
+ double **v = atom->v;
+ double *duCond = atom->duCond;
+ double *duMech = atom->duMech;
+ double *uCond = atom->uCond;
+ double *uMech = atom->uMech;
+
+ m = 0;
+ for (ii = 0; ii < n; ii++) {
+ jj = list[ii];
+ buf[m++] = dvSSA[jj][0];
+ buf[m++] = dvSSA[jj][1];
+ buf[m++] = dvSSA[jj][2];
+ buf[m++] = v[jj][0];
+ buf[m++] = v[jj][1];
+ buf[m++] = v[jj][2];
+ if(pairDPDE){
+ buf[m++] = duCond[jj];
+ buf[m++] = duMech[jj];
+ buf[m++] = uCond[jj];
+ buf[m++] = uMech[jj];
+ }
+ }
+ return m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixShardlow::unpack_forward_comm(int n, int first, double *buf)
+{
+ int ii,m,last;
+ double **v = atom->v;
+ double *duCond = atom->duCond;
+ double *duMech = atom->duMech;
+ double *uCond = atom->uCond;
+ double *uMech = atom->uMech;
+
+ m = 0;
+ last = first + n ;
+ for (ii = first; ii < last; ii++) {
+ dvSSA[ii][0] = buf[m++];
+ dvSSA[ii][1] = buf[m++];
+ dvSSA[ii][2] = buf[m++];
+ v[ii][0] = buf[m++];
+ v[ii][1] = buf[m++];
+ v[ii][2] = buf[m++];
+ if(pairDPDE){
+ duCond[ii] = buf[m++];
+ duMech[ii] = buf[m++];
+ uCond[ii] = buf[m++];
+ uMech[ii] = buf[m++];
+ }
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixShardlow::pack_reverse_comm(int n, int first, double *buf)
+{
+ int i,m,last;
+ double *duCond = atom->duCond;
+ double *duMech = atom->duMech;
+
+ m = 0;
+ last = first + n;
+ for (i = first; i < last; i++) {
+ buf[m++] = dvSSA[i][0];
+ buf[m++] = dvSSA[i][1];
+ buf[m++] = dvSSA[i][2];
+ if(pairDPDE){
+ buf[m++] = duCond[i];
+ buf[m++] = duMech[i];
+ }
+ }
+ return m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixShardlow::unpack_reverse_comm(int n, int *list, double *buf)
+{
+ int i,j,m;
+ double *duCond = atom->duCond;
+ double *duMech = atom->duMech;
+
+ m = 0;
+ for (i = 0; i < n; i++) {
+ j = list[i];
+
+ dvSSA[j][0] += buf[m++];
+ dvSSA[j][1] += buf[m++];
+ dvSSA[j][2] += buf[m++];
+ if(pairDPDE){
+ duCond[j] += buf[m++];
+ duMech[j] += buf[m++];
+ }
+ }
+}
diff --git a/src/USER-DPD/fix_shardlow.h b/src/USER-DPD/fix_shardlow.h
new file mode 100644
index 0000000000000000000000000000000000000000..6c25fd86e28b6e83490d3ac2f98685ac3bcbedf1
--- /dev/null
+++ b/src/USER-DPD/fix_shardlow.h
@@ -0,0 +1,115 @@
+/* ----------------------------------------------------------------------
+ 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 authors:
+ James Larentzos and Timothy I. Mattox (Engility Corporation)
+
+ Martin Lisal (Institute of Chemical Process Fundamentals
+ of the Czech Academy of Sciences and J. E. Purkinje University)
+
+ John Brennan, Joshua Moore and William Mattson (Army Research Lab)
+
+ Please cite the related publications:
+ J. P. Larentzos, J. K. Brennan, J. D. Moore, M. Lisal, W. D. Mattson,
+ "Parallel implementation of isothermal and isoenergetic Dissipative
+ Particle Dynamics using Shardlow-like splitting algorithms",
+ Computer Physics Communications, 2014, 185, pp 1987--1998.
+
+ M. Lisal, J. K. Brennan, J. Bonet Avalos, "Dissipative particle dynamics
+ at isothermal, isobaric, isoenergetic, and isoenthalpic conditions using
+ Shardlow-like splitting algorithms", Journal of Chemical Physics, 2011,
+ 135, 204105.
+------------------------------------------------------------------------- */
+
+#ifdef FIX_CLASS
+
+FixStyle(shardlow,FixShardlow)
+
+#else
+
+#ifndef LMP_FIX_SHARDLOW_H
+#define LMP_FIX_SHARDLOW_H
+
+#include "fix.h"
+
+namespace LAMMPS_NS {
+
+class FixShardlow : public Fix {
+ public:
+ FixShardlow(class LAMMPS *, int, char **);
+ virtual ~FixShardlow() {}
+ int setmask();
+ virtual void init_list(int,class NeighList *);
+ virtual void setup(int);
+ virtual void setup_pre_force(int);
+ virtual void initial_integrate(int);
+
+ void setup_pre_neighbor();
+ void pre_neighbor();
+
+ protected:
+ int pack_reverse_comm(int, int, double *);
+ void unpack_reverse_comm(int, int *, double *);
+ int pack_forward_comm(int , int *, double *, int, int *);
+ void unpack_forward_comm(int , int , double *);
+
+ class PairDPDfdt *pairDPD;
+ class PairDPDfdtEnergy *pairDPDE;
+ double **dvSSA;
+
+ private:
+ class NeighList *list;
+
+};
+
+}
+
+#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: Must use dpd/fdt pair_style with fix shardlow
+
+Self-explanatory.
+
+E: Must use pair_style dpd/fdt or dpd/fdt/energy with fix shardlow
+
+E: A deterministic integrator must be specified after fix shardlow in input
+file (e.g. fix nve or fix nph).
+
+Self-explanatory.
+
+E: Cannot use constant temperature integration routines with DPD
+
+Self-explanatory. Must use deterministic integrators such as nve or nph
+
+E: Fix shardlow does not yet support triclinic geometries
+
+Self-explanatory.
+
+E: Shardlow algorithm requires sub-domain length > 2*(rcut+skin). Either
+reduce the number of processors requested, or change the cutoff/skin
+
+The Shardlow splitting algorithm requires the size of the sub-domain lengths
+to be are larger than twice the cutoff+skin. Generally, the domain decomposition
+is dependant on the number of processors requested.
+
+*/
diff --git a/src/USER-DPD/pair_dpd_conservative.cpp b/src/USER-DPD/pair_dpd_conservative.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..870702311085e7e94227f5000467254112916e34
--- /dev/null
+++ b/src/USER-DPD/pair_dpd_conservative.cpp
@@ -0,0 +1,320 @@
+/* ----------------------------------------------------------------------
+ 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: James Larentzos (Engility Corporation)
+------------------------------------------------------------------------- */
+
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "pair_dpd_conservative.h"
+#include "atom.h"
+#include "atom_vec.h"
+#include "comm.h"
+#include "update.h"
+#include "fix.h"
+#include "force.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "neigh_request.h"
+#include "memory.h"
+#include "modify.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+#define EPSILON 1.0e-10
+
+/* ---------------------------------------------------------------------- */
+
+PairDPDconservative::PairDPDconservative(LAMMPS *lmp) : Pair(lmp)
+{
+}
+
+/* ---------------------------------------------------------------------- */
+
+PairDPDconservative::~PairDPDconservative()
+{
+ if (allocated) {
+ memory->destroy(setflag);
+ memory->destroy(cutsq);
+
+ memory->destroy(cut);
+ memory->destroy(a0);
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairDPDconservative::compute(int eflag, int vflag)
+{
+ int i,j,ii,jj,inum,jnum,itype,jtype;
+ double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
+ double rsq,r,rinv,wd,wr,factor_dpd;
+ int *ilist,*jlist,*numneigh,**firstneigh;
+
+ evdwl = 0.0;
+ if (eflag || vflag) ev_setup(eflag,vflag);
+ else evflag = vflag_fdotr = 0;
+
+ double **x = atom->x;
+ double **f = atom->f;
+ int *type = atom->type;
+ int nlocal = atom->nlocal;
+ double *special_lj = force->special_lj;
+ int newton_pair = force->newton_pair;
+
+ inum = list->inum;
+ ilist = list->ilist;
+ numneigh = list->numneigh;
+ firstneigh = list->firstneigh;
+
+ // loop over neighbors of my atoms
+
+ for (ii = 0; ii < inum; ii++) {
+ i = ilist[ii];
+ xtmp = x[i][0];
+ ytmp = x[i][1];
+ ztmp = x[i][2];
+ itype = type[i];
+ jlist = firstneigh[i];
+ jnum = numneigh[i];
+
+ for (jj = 0; jj < jnum; jj++) {
+ j = jlist[jj];
+ factor_dpd = special_lj[sbmask(j)];
+ j &= NEIGHMASK;
+
+ delx = xtmp - x[j][0];
+ dely = ytmp - x[j][1];
+ delz = ztmp - x[j][2];
+ rsq = delx*delx + dely*dely + delz*delz;
+ jtype = type[j];
+
+ if (rsq < cutsq[itype][jtype]) {
+ r = sqrt(rsq);
+ 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;
+
+ // conservative force = a0 * wr
+ fpair = a0[itype][jtype]*wr;
+ fpair *= factor_dpd*rinv;
+
+ f[i][0] += delx*fpair;
+ f[i][1] += dely*fpair;
+ f[i][2] += delz*fpair;
+ if (newton_pair || j < nlocal) {
+ f[j][0] -= delx*fpair;
+ f[j][1] -= dely*fpair;
+ f[j][2] -= delz*fpair;
+ }
+
+ if (eflag) {
+ // unshifted eng of conservative term:
+ // evdwl = -a0[itype][jtype]*r * (1.0-0.5*r/cut[itype][jtype]);
+ // eng shifted to 0.0 at cutoff
+ evdwl = 0.5*a0[itype][jtype]*cut[itype][jtype] * wd;
+ evdwl *= factor_dpd;
+ }
+
+ if (evflag) ev_tally(i,j,nlocal,newton_pair,
+ evdwl,0.0,fpair,delx,dely,delz);
+ }
+ }
+ }
+
+ if (vflag_fdotr) virial_fdotr_compute();
+}
+
+/* ----------------------------------------------------------------------
+ allocate all arrays
+------------------------------------------------------------------------- */
+
+void PairDPDconservative::allocate()
+{
+ allocated = 1;
+ int n = atom->ntypes;
+
+ memory->create(setflag,n+1,n+1,"pair:setflag");
+ for (int i = 1; i <= n; i++)
+ for (int j = i; j <= n; j++)
+ setflag[i][j] = 0;
+
+ memory->create(cutsq,n+1,n+1,"pair:cutsq");
+ memory->create(cut,n+1,n+1,"pair:cut");
+ memory->create(a0,n+1,n+1,"pair:a0");
+
+}
+
+/* ----------------------------------------------------------------------
+ global settings
+------------------------------------------------------------------------- */
+
+void PairDPDconservative::settings(int narg, char **arg)
+{
+ // process keywords
+ if (narg != 1) error->all(FLERR,"Illegal pair_style command");
+ cut_global = force->numeric(FLERR,arg[0]);
+
+ // reset cutoffs that have been explicitly set
+
+ if (allocated) {
+ int i,j;
+ for (i = 1; i <= atom->ntypes; i++)
+ for (j = i+1; j <= atom->ntypes; j++)
+ if (setflag[i][j]) cut[i][j] = cut_global;
+ }
+}
+
+/* ----------------------------------------------------------------------
+ set coeffs for one or more type pairs
+------------------------------------------------------------------------- */
+
+void PairDPDconservative::coeff(int narg, char **arg)
+{
+ if (narg < 3 || narg > 4) error->all(FLERR,"Incorrect args for pair coefficients");
+ if (!allocated) allocate();
+
+ int ilo,ihi,jlo,jhi;
+ force->bounds(arg[0],atom->ntypes,ilo,ihi);
+ force->bounds(arg[1],atom->ntypes,jlo,jhi);
+
+ double a0_one = force->numeric(FLERR,arg[2]);
+ double cut_one = cut_global;
+
+ if (narg == 4) cut_one = force->numeric(FLERR,arg[3]);
+
+ int count = 0;
+ for (int i = ilo; i <= ihi; i++) {
+ for (int j = MAX(jlo,i); j <= jhi; j++) {
+ a0[i][j] = a0_one;
+ cut[i][j] = cut_one;
+ setflag[i][j] = 1;
+ count++;
+ }
+ }
+
+ if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
+}
+
+/* ----------------------------------------------------------------------
+ init for one type pair i,j and corresponding j,i
+------------------------------------------------------------------------- */
+
+double PairDPDconservative::init_one(int i, int j)
+{
+ if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
+
+ cut[j][i] = cut[i][j];
+ a0[j][i] = a0[i][j];
+
+ return cut[i][j];
+}
+
+/* ----------------------------------------------------------------------
+ proc 0 writes to restart file
+------------------------------------------------------------------------- */
+
+void PairDPDconservative::write_restart(FILE *fp)
+{
+ write_restart_settings(fp);
+
+ int i,j;
+ for (i = 1; i <= atom->ntypes; i++)
+ for (j = i; j <= atom->ntypes; j++) {
+ fwrite(&setflag[i][j],sizeof(int),1,fp);
+ if (setflag[i][j]) {
+ fwrite(&a0[i][j],sizeof(double),1,fp);
+ fwrite(&cut[i][j],sizeof(double),1,fp);
+ }
+ }
+}
+
+/* ----------------------------------------------------------------------
+ proc 0 reads from restart file, bcasts
+------------------------------------------------------------------------- */
+
+void PairDPDconservative::read_restart(FILE *fp)
+{
+ read_restart_settings(fp);
+
+ allocate();
+
+ int i,j;
+ int me = comm->me;
+ for (i = 1; i <= atom->ntypes; i++)
+ for (j = i; j <= atom->ntypes; j++) {
+ if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp);
+ MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
+ if (setflag[i][j]) {
+ if (me == 0) {
+ fread(&a0[i][j],sizeof(double),1,fp);
+ fread(&cut[i][j],sizeof(double),1,fp);
+ }
+ MPI_Bcast(&a0[i][j],1,MPI_DOUBLE,0,world);
+ MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world);
+ }
+ }
+}
+
+/* ----------------------------------------------------------------------
+ proc 0 writes to restart file
+------------------------------------------------------------------------- */
+
+void PairDPDconservative::write_restart_settings(FILE *fp)
+{
+ fwrite(&cut_global,sizeof(double),1,fp);
+ fwrite(&mix_flag,sizeof(int),1,fp);
+}
+
+/* ----------------------------------------------------------------------
+ proc 0 reads from restart file, bcasts
+------------------------------------------------------------------------- */
+
+void PairDPDconservative::read_restart_settings(FILE *fp)
+{
+ if (comm->me == 0) {
+ fread(&cut_global,sizeof(double),1,fp);
+ fread(&mix_flag,sizeof(int),1,fp);
+ }
+ MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world);
+ MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
+
+}
+
+/* ---------------------------------------------------------------------- */
+
+double PairDPDconservative::single(int i, int j, int itype, int jtype, double rsq,
+ double factor_coul, double factor_dpd, double &fforce)
+{
+ double r,rinv,wr,wd,phi;
+
+ r = sqrt(rsq);
+ if (r < EPSILON) {
+ fforce = 0.0;
+ return 0.0;
+ }
+
+ rinv = 1.0/r;
+ wr = 1.0 - r/cut[itype][jtype];
+ wd = wr*wr;
+ fforce = a0[itype][jtype]*wr * factor_dpd*rinv;
+
+ phi = 0.5*a0[itype][jtype]*cut[itype][jtype] * wd;
+ return factor_dpd*phi;
+}
+
diff --git a/src/USER-DPD/pair_dpd_conservative.h b/src/USER-DPD/pair_dpd_conservative.h
new file mode 100644
index 0000000000000000000000000000000000000000..15402be6d8bb69954354d2a7ee0ec55531e0f917
--- /dev/null
+++ b/src/USER-DPD/pair_dpd_conservative.h
@@ -0,0 +1,77 @@
+/* ----------------------------------------------------------------------
+ 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: James Larentzos (Engility Corporation)
+------------------------------------------------------------------------- */
+
+#ifdef PAIR_CLASS
+
+PairStyle(dpd/conservative,PairDPDconservative)
+
+#else
+
+#ifndef LMP_PAIR_DPD_CONSERVATIVE_H
+#define LMP_PAIR_DPD_CONSERVATIVE_H
+
+#include "pair.h"
+
+namespace LAMMPS_NS {
+
+class PairDPDconservative : public Pair {
+ public:
+ PairDPDconservative(class LAMMPS *);
+ virtual ~PairDPDconservative();
+ virtual void compute(int, int);
+ virtual void settings(int, char **);
+ virtual void coeff(int, char **);
+ double init_one(int, int);
+ void write_restart(FILE *);
+ void read_restart(FILE *);
+ virtual void write_restart_settings(FILE *);
+ virtual void read_restart_settings(FILE *);
+ double single(int, int, int, int, double, double, double, double &);
+
+ double **cut;
+ double **a0;
+
+ protected:
+ double cut_global;
+
+ void allocate();
+
+};
+
+}
+
+#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: Incorrect args for pair coefficients
+
+Self-explanatory. Check the input script or data file.
+
+E: All pair coeffs are not set
+
+All pair coefficients must be set in the data file or by the
+pair_coeff command before running a simulation.
+
+*/
diff --git a/src/USER-DPD/pair_dpd_fdt.cpp b/src/USER-DPD/pair_dpd_fdt.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e03d913f9f3524155f68cfd4971c283dbfa71535
--- /dev/null
+++ b/src/USER-DPD/pair_dpd_fdt.cpp
@@ -0,0 +1,376 @@
+/* ----------------------------------------------------------------------
+ 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: James Larentzos (Engility Corporation)
+------------------------------------------------------------------------- */
+
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "atom.h"
+#include "atom_vec.h"
+#include "comm.h"
+#include "update.h"
+#include "fix.h"
+#include "force.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "neigh_request.h"
+#include "random_mars.h"
+#include "memory.h"
+#include "modify.h"
+#include "pair_dpd_fdt.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+#define EPSILON 1.0e-10
+
+/* ---------------------------------------------------------------------- */
+
+PairDPDfdt::PairDPDfdt(LAMMPS *lmp) : Pair(lmp)
+{
+ random = NULL;
+}
+
+/* ---------------------------------------------------------------------- */
+
+PairDPDfdt::~PairDPDfdt()
+{
+ if (allocated) {
+ memory->destroy(setflag);
+ memory->destroy(cutsq);
+
+ memory->destroy(cut);
+ memory->destroy(a0);
+ memory->destroy(sigma);
+ }
+
+
+ if (random) delete random;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairDPDfdt::compute(int eflag, int vflag)
+{
+ int i,j,ii,jj,inum,jnum,itype,jtype;
+ double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
+ double rsq,r,rinv,wd,wr,factor_dpd;
+ int *ilist,*jlist,*numneigh,**firstneigh;
+
+ evdwl = 0.0;
+ if (eflag || vflag) ev_setup(eflag,vflag);
+ else evflag = vflag_fdotr = 0;
+
+ double **x = atom->x;
+ double **f = atom->f;
+ int *type = atom->type;
+ int nlocal = atom->nlocal;
+ double *special_lj = force->special_lj;
+ int newton_pair = force->newton_pair;
+
+ inum = list->inum;
+ ilist = list->ilist;
+ numneigh = list->numneigh;
+ firstneigh = list->firstneigh;
+
+ // loop over neighbors of my atoms
+
+ for (ii = 0; ii < inum; ii++) {
+ i = ilist[ii];
+ xtmp = x[i][0];
+ ytmp = x[i][1];
+ ztmp = x[i][2];
+ itype = type[i];
+ jlist = firstneigh[i];
+ jnum = numneigh[i];
+
+ for (jj = 0; jj < jnum; jj++) {
+ j = jlist[jj];
+ factor_dpd = special_lj[sbmask(j)];
+ j &= NEIGHMASK;
+
+ delx = xtmp - x[j][0];
+ dely = ytmp - x[j][1];
+ delz = ztmp - x[j][2];
+ rsq = delx*delx + dely*dely + delz*delz;
+ jtype = type[j];
+
+ if (rsq < cutsq[itype][jtype]) {
+ r = sqrt(rsq);
+ 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;
+
+ // conservative force = a0 * wr
+ fpair = a0[itype][jtype]*wr;
+ fpair *= factor_dpd*rinv;
+
+ f[i][0] += delx*fpair;
+ f[i][1] += dely*fpair;
+ f[i][2] += delz*fpair;
+ if (newton_pair || j < nlocal) {
+ f[j][0] -= delx*fpair;
+ f[j][1] -= dely*fpair;
+ f[j][2] -= delz*fpair;
+ }
+
+ if (eflag) {
+ // unshifted eng of conservative term:
+ // evdwl = -a0[itype][jtype]*r * (1.0-0.5*r/cut[itype][jtype]);
+ // eng shifted to 0.0 at cutoff
+ evdwl = 0.5*a0[itype][jtype]*cut[itype][jtype] * wd;
+ evdwl *= factor_dpd;
+ }
+
+ if (evflag) ev_tally(i,j,nlocal,newton_pair,
+ evdwl,0.0,fpair,delx,dely,delz);
+ }
+ }
+ }
+
+ if (vflag_fdotr) virial_fdotr_compute();
+}
+
+/* ----------------------------------------------------------------------
+ allocate all arrays
+------------------------------------------------------------------------- */
+
+void PairDPDfdt::allocate()
+{
+ allocated = 1;
+ int n = atom->ntypes;
+
+ memory->create(setflag,n+1,n+1,"pair:setflag");
+ for (int i = 1; i <= n; i++)
+ for (int j = i; j <= n; j++)
+ setflag[i][j] = 0;
+
+ memory->create(cutsq,n+1,n+1,"pair:cutsq");
+
+ memory->create(cut,n+1,n+1,"pair:cut");
+ memory->create(a0,n+1,n+1,"pair:a0");
+ memory->create(sigma,n+1,n+1,"pair:sigma");
+}
+
+/* ----------------------------------------------------------------------
+ global settings
+------------------------------------------------------------------------- */
+
+void PairDPDfdt::settings(int narg, char **arg)
+{
+ // process keywords
+ if (narg != 3) error->all(FLERR,"Illegal pair_style command");
+
+ temperature = force->numeric(FLERR,arg[0]);
+ cut_global = force->numeric(FLERR,arg[1]);
+ seed = force->inumeric(FLERR,arg[2]);
+
+ // initialize Marsaglia RNG with processor-unique seed
+
+ if (seed <= 0) error->all(FLERR,"Illegal pair_style command");
+ delete random;
+ random = new RanMars(lmp,seed + comm->me);
+
+ // reset cutoffs that have been explicitly set
+
+ if (allocated) {
+ int i,j;
+ for (i = 1; i <= atom->ntypes; i++)
+ for (j = i+1; j <= atom->ntypes; j++)
+ if (setflag[i][j]) cut[i][j] = cut_global;
+ }
+}
+
+/* ----------------------------------------------------------------------
+ set coeffs for one or more type pairs
+------------------------------------------------------------------------- */
+
+void PairDPDfdt::coeff(int narg, char **arg)
+{
+ if (narg < 4 || narg > 5) error->all(FLERR,"Incorrect args for pair coefficients");
+ if (!allocated) allocate();
+
+ int ilo,ihi,jlo,jhi;
+ force->bounds(arg[0],atom->ntypes,ilo,ihi);
+ force->bounds(arg[1],atom->ntypes,jlo,jhi);
+
+ double a0_one = force->numeric(FLERR,arg[2]);
+ double sigma_one = force->numeric(FLERR,arg[3]);
+ double cut_one = cut_global;
+ double kappa_one;
+
+ if (narg == 5) cut_one = force->numeric(FLERR,arg[4]);
+
+ int count = 0;
+ for (int i = ilo; i <= ihi; i++) {
+ for (int j = MAX(jlo,i); j <= jhi; j++) {
+ a0[i][j] = a0_one;
+ sigma[i][j] = sigma_one;
+ cut[i][j] = cut_one;
+ setflag[i][j] = 1;
+ count++;
+ }
+ }
+
+ if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
+}
+
+/* ----------------------------------------------------------------------
+ init specific to this pair style
+------------------------------------------------------------------------- */
+
+void PairDPDfdt::init_style()
+{
+ if (comm->ghost_velocity == 0)
+ error->all(FLERR,"Pair dpd/fdt requires ghost atoms store velocity");
+
+ // if newton off, forces between atoms ij will be double computed
+ // using different random numbers
+
+ if (force->newton_pair == 0 && comm->me == 0) error->warning(FLERR,
+ "Pair dpd/fdt requires newton pair on");
+
+ int irequest = neighbor->request(this,instance_me);
+ neighbor->requests[irequest]->ssa = 0;
+ for (int i = 0; i < modify->nfix; i++)
+ if (strcmp(modify->fix[i]->style,"shardlow") == 0)
+ neighbor->requests[irequest]->ssa = 1;
+}
+
+/* ----------------------------------------------------------------------
+ init for one type pair i,j and corresponding j,i
+------------------------------------------------------------------------- */
+
+double PairDPDfdt::init_one(int i, int j)
+{
+ if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
+
+ cut[j][i] = cut[i][j];
+ a0[j][i] = a0[i][j];
+ sigma[j][i] = sigma[i][j];
+
+ return cut[i][j];
+}
+
+/* ----------------------------------------------------------------------
+ proc 0 writes to restart file
+------------------------------------------------------------------------- */
+
+void PairDPDfdt::write_restart(FILE *fp)
+{
+ write_restart_settings(fp);
+
+ int i,j;
+ for (i = 1; i <= atom->ntypes; i++)
+ for (j = i; j <= atom->ntypes; j++) {
+ fwrite(&setflag[i][j],sizeof(int),1,fp);
+ if (setflag[i][j]) {
+ fwrite(&a0[i][j],sizeof(double),1,fp);
+ fwrite(&sigma[i][j],sizeof(double),1,fp);
+ fwrite(&cut[i][j],sizeof(double),1,fp);
+ }
+ }
+}
+
+/* ----------------------------------------------------------------------
+ proc 0 reads from restart file, bcasts
+------------------------------------------------------------------------- */
+
+void PairDPDfdt::read_restart(FILE *fp)
+{
+ read_restart_settings(fp);
+
+ allocate();
+
+ int i,j;
+ int me = comm->me;
+ for (i = 1; i <= atom->ntypes; i++)
+ for (j = i; j <= atom->ntypes; j++) {
+ if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp);
+ MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
+ if (setflag[i][j]) {
+ if (me == 0) {
+ fread(&a0[i][j],sizeof(double),1,fp);
+ fread(&sigma[i][j],sizeof(double),1,fp);
+ fread(&cut[i][j],sizeof(double),1,fp);
+ }
+ MPI_Bcast(&a0[i][j],1,MPI_DOUBLE,0,world);
+ MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world);
+ MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world);
+ }
+ }
+}
+
+/* ----------------------------------------------------------------------
+ proc 0 writes to restart file
+------------------------------------------------------------------------- */
+
+void PairDPDfdt::write_restart_settings(FILE *fp)
+{
+ fwrite(&temperature,sizeof(double),1,fp);
+ fwrite(&cut_global,sizeof(double),1,fp);
+ fwrite(&seed,sizeof(int),1,fp);
+ fwrite(&mix_flag,sizeof(int),1,fp);
+}
+
+/* ----------------------------------------------------------------------
+ proc 0 reads from restart file, bcasts
+------------------------------------------------------------------------- */
+
+void PairDPDfdt::read_restart_settings(FILE *fp)
+{
+ if (comm->me == 0) {
+ fread(&temperature,sizeof(double),1,fp);
+ fread(&cut_global,sizeof(double),1,fp);
+ fread(&seed,sizeof(int),1,fp);
+ fread(&mix_flag,sizeof(int),1,fp);
+ }
+ MPI_Bcast(&temperature,1,MPI_DOUBLE,0,world);
+ MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world);
+ MPI_Bcast(&seed,1,MPI_INT,0,world);
+ MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
+
+ // initialize Marsaglia RNG with processor-unique seed
+ // same seed that pair_style command initially specified
+
+ if (random) delete random;
+ random = new RanMars(lmp,seed + comm->me);
+}
+
+/* ---------------------------------------------------------------------- */
+
+double PairDPDfdt::single(int i, int j, int itype, int jtype, double rsq,
+ double factor_coul, double factor_dpd, double &fforce)
+{
+ double r,rinv,wr,wd,phi;
+
+ r = sqrt(rsq);
+ if (r < EPSILON) {
+ fforce = 0.0;
+ return 0.0;
+ }
+
+ rinv = 1.0/r;
+ wr = 1.0 - r/cut[itype][jtype];
+ wd = wr*wr;
+ fforce = a0[itype][jtype]*wr * factor_dpd*rinv;
+
+ phi = 0.5*a0[itype][jtype]*cut[itype][jtype] * wd;
+ return factor_dpd*phi;
+}
+
diff --git a/src/USER-DPD/pair_dpd_fdt.h b/src/USER-DPD/pair_dpd_fdt.h
new file mode 100644
index 0000000000000000000000000000000000000000..8384e3b6828f16fb92f0af9b8551eb27a0f9a879
--- /dev/null
+++ b/src/USER-DPD/pair_dpd_fdt.h
@@ -0,0 +1,91 @@
+/* ----------------------------------------------------------------------
+ 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: James Larentzos (Engility Corporation)
+------------------------------------------------------------------------- */
+
+#ifdef PAIR_CLASS
+
+PairStyle(dpd/fdt,PairDPDfdt)
+
+#else
+
+#ifndef LMP_PAIR_DPD_FDT_H
+#define LMP_PAIR_DPD_FDT_H
+
+#include "pair.h"
+
+namespace LAMMPS_NS {
+
+class PairDPDfdt : public Pair {
+ public:
+ PairDPDfdt(class LAMMPS *);
+ virtual ~PairDPDfdt();
+ void compute(int, int);
+ virtual void settings(int, char **);
+ virtual void coeff(int, char **);
+ void init_style();
+ double init_one(int, int);
+ void write_restart(FILE *);
+ void read_restart(FILE *);
+ virtual void write_restart_settings(FILE *);
+ virtual void read_restart_settings(FILE *);
+ double single(int, int, int, int, double, double, double, double &);
+
+ double **cut;
+ double **a0;
+ double **sigma;
+ double temperature;
+
+ class RanMars *random;
+
+ protected:
+ double cut_global;
+ int seed;
+
+ void allocate();
+
+};
+
+}
+
+#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: Incorrect args for pair coefficients
+
+Self-explanatory. Check the input script or data file.
+
+E: Pair dpd/fdt requires ghost atoms store velocity
+
+Use the communicate vel yes command to enable this.
+
+E: Pair dpd/fdt requires newton pair on
+
+Self-explanatory.
+
+E: All pair coeffs are not set
+
+All pair coefficients must be set in the data file or by the
+pair_coeff command before running a simulation.
+
+*/
diff --git a/src/USER-DPD/pair_dpd_fdt_energy.cpp b/src/USER-DPD/pair_dpd_fdt_energy.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..8d71479f0bfa2438a04dc7e7e5203e62233d1b0f
--- /dev/null
+++ b/src/USER-DPD/pair_dpd_fdt_energy.cpp
@@ -0,0 +1,382 @@
+/* ----------------------------------------------------------------------
+ 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: James Larentzos (Engility Corporation)
+------------------------------------------------------------------------- */
+
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "atom.h"
+#include "atom_vec.h"
+#include "comm.h"
+#include "update.h"
+#include "fix.h"
+#include "force.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "neigh_request.h"
+#include "random_mars.h"
+#include "memory.h"
+#include "modify.h"
+#include "pair_dpd_fdt_energy.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+#define EPSILON 1.0e-10
+
+/* ---------------------------------------------------------------------- */
+
+PairDPDfdtEnergy::PairDPDfdtEnergy(LAMMPS *lmp) : Pair(lmp)
+{
+ random = NULL;
+}
+
+/* ---------------------------------------------------------------------- */
+
+PairDPDfdtEnergy::~PairDPDfdtEnergy()
+{
+ if (allocated) {
+ memory->destroy(setflag);
+ memory->destroy(cutsq);
+
+ memory->destroy(cut);
+ memory->destroy(a0);
+ memory->destroy(sigma);
+ memory->destroy(kappa);
+ }
+
+
+ if (random) delete random;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairDPDfdtEnergy::compute(int eflag, int vflag)
+{
+ int i,j,ii,jj,inum,jnum,itype,jtype;
+ double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
+ double rsq,r,rinv,wd,wr,factor_dpd;
+ int *ilist,*jlist,*numneigh,**firstneigh;
+
+ evdwl = 0.0;
+ if (eflag || vflag) ev_setup(eflag,vflag);
+ else evflag = vflag_fdotr = 0;
+
+ double **x = atom->x;
+ double **f = atom->f;
+ int *type = atom->type;
+ int nlocal = atom->nlocal;
+ double *special_lj = force->special_lj;
+ int newton_pair = force->newton_pair;
+
+ inum = list->inum;
+ ilist = list->ilist;
+ numneigh = list->numneigh;
+ firstneigh = list->firstneigh;
+
+ // loop over neighbors of my atoms
+
+ for (ii = 0; ii < inum; ii++) {
+ i = ilist[ii];
+ xtmp = x[i][0];
+ ytmp = x[i][1];
+ ztmp = x[i][2];
+ itype = type[i];
+ jlist = firstneigh[i];
+ jnum = numneigh[i];
+
+ for (jj = 0; jj < jnum; jj++) {
+ j = jlist[jj];
+ factor_dpd = special_lj[sbmask(j)];
+ j &= NEIGHMASK;
+
+ delx = xtmp - x[j][0];
+ dely = ytmp - x[j][1];
+ delz = ztmp - x[j][2];
+ rsq = delx*delx + dely*dely + delz*delz;
+ jtype = type[j];
+
+ if (rsq < cutsq[itype][jtype]) {
+ r = sqrt(rsq);
+ 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;
+
+ // conservative force = a0 * wr
+ fpair = a0[itype][jtype]*wr;
+ fpair *= factor_dpd*rinv;
+
+ f[i][0] += delx*fpair;
+ f[i][1] += dely*fpair;
+ f[i][2] += delz*fpair;
+ if (newton_pair || j < nlocal) {
+ f[j][0] -= delx*fpair;
+ f[j][1] -= dely*fpair;
+ f[j][2] -= delz*fpair;
+ }
+
+ if (eflag) {
+ // unshifted eng of conservative term:
+ // evdwl = -a0[itype][jtype]*r * (1.0-0.5*r/cut[itype][jtype]);
+ // eng shifted to 0.0 at cutoff
+ evdwl = 0.5*a0[itype][jtype]*cut[itype][jtype] * wd;
+ evdwl *= factor_dpd;
+ }
+
+ if (evflag) ev_tally(i,j,nlocal,newton_pair,
+ evdwl,0.0,fpair,delx,dely,delz);
+ }
+ }
+ }
+
+ if (vflag_fdotr) virial_fdotr_compute();
+}
+
+/* ----------------------------------------------------------------------
+ allocate all arrays
+------------------------------------------------------------------------- */
+
+void PairDPDfdtEnergy::allocate()
+{
+ allocated = 1;
+ int n = atom->ntypes;
+
+ memory->create(setflag,n+1,n+1,"pair:setflag");
+ for (int i = 1; i <= n; i++)
+ for (int j = i; j <= n; j++)
+ setflag[i][j] = 0;
+
+ memory->create(cutsq,n+1,n+1,"pair:cutsq");
+
+ memory->create(cut,n+1,n+1,"pair:cut");
+ memory->create(a0,n+1,n+1,"pair:a0");
+ memory->create(sigma,n+1,n+1,"pair:sigma");
+ memory->create(kappa,n+1,n+1,"pair:kappa");
+}
+
+/* ----------------------------------------------------------------------
+ global settings
+------------------------------------------------------------------------- */
+
+void PairDPDfdtEnergy::settings(int narg, char **arg)
+{
+ // process keywords
+ if (narg != 2) error->all(FLERR,"Illegal pair_style command");
+
+ cut_global = force->numeric(FLERR,arg[0]);
+ seed = force->inumeric(FLERR,arg[1]);
+ if (atom->dpd_flag != 1)
+ error->all(FLERR,"pair_style dpd/fdt/energy requires atom_style with internal temperature and energies (e.g. dpd)");
+
+ // initialize Marsaglia RNG with processor-unique seed
+
+ if (seed <= 0) error->all(FLERR,"Illegal pair_style command");
+ delete random;
+ random = new RanMars(lmp,seed + comm->me);
+
+ // reset cutoffs that have been explicitly set
+
+ if (allocated) {
+ int i,j;
+ for (i = 1; i <= atom->ntypes; i++)
+ for (j = i+1; j <= atom->ntypes; j++)
+ if (setflag[i][j]) cut[i][j] = cut_global;
+ }
+}
+
+/* ----------------------------------------------------------------------
+ set coeffs for one or more type pairs
+------------------------------------------------------------------------- */
+
+void PairDPDfdtEnergy::coeff(int narg, char **arg)
+{
+ if (narg < 5 || narg > 6) error->all(FLERR,"Incorrect args for pair coefficients");
+ if (!allocated) allocate();
+
+ int ilo,ihi,jlo,jhi;
+ force->bounds(arg[0],atom->ntypes,ilo,ihi);
+ force->bounds(arg[1],atom->ntypes,jlo,jhi);
+
+ double a0_one = force->numeric(FLERR,arg[2]);
+ double sigma_one = force->numeric(FLERR,arg[3]);
+ double cut_one = cut_global;
+ double kappa_one;
+
+ kappa_one = force->numeric(FLERR,arg[4]);
+ if (narg == 6) cut_one = force->numeric(FLERR,arg[5]);
+
+ int count = 0;
+ for (int i = ilo; i <= ihi; i++) {
+ for (int j = MAX(jlo,i); j <= jhi; j++) {
+ a0[i][j] = a0_one;
+ sigma[i][j] = sigma_one;
+ kappa[i][j] = kappa_one;
+ cut[i][j] = cut_one;
+ setflag[i][j] = 1;
+ count++;
+ }
+ }
+
+ if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
+}
+
+/* ----------------------------------------------------------------------
+ init specific to this pair style
+------------------------------------------------------------------------- */
+
+void PairDPDfdtEnergy::init_style()
+{
+ if (comm->ghost_velocity == 0)
+ error->all(FLERR,"Pair dpd/fdt/energy requires ghost atoms store velocity");
+
+ // if newton off, forces between atoms ij will be double computed
+ // using different random numbers
+
+ if (force->newton_pair == 0 && comm->me == 0) error->warning(FLERR,
+ "Pair dpd/fdt/energy requires newton pair on");
+
+ int irequest = neighbor->request(this,instance_me);
+ neighbor->requests[irequest]->ssa = 0;
+ for (int i = 0; i < modify->nfix; i++)
+ if (strcmp(modify->fix[i]->style,"shardlow") == 0)
+ neighbor->requests[irequest]->ssa = 1;
+}
+
+/* ----------------------------------------------------------------------
+ init for one type pair i,j and corresponding j,i
+------------------------------------------------------------------------- */
+
+double PairDPDfdtEnergy::init_one(int i, int j)
+{
+ if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
+
+ cut[j][i] = cut[i][j];
+ a0[j][i] = a0[i][j];
+ sigma[j][i] = sigma[i][j];
+ kappa[j][i] = kappa[i][j];
+
+ return cut[i][j];
+}
+
+/* ----------------------------------------------------------------------
+ proc 0 writes to restart file
+------------------------------------------------------------------------- */
+
+void PairDPDfdtEnergy::write_restart(FILE *fp)
+{
+ write_restart_settings(fp);
+
+ int i,j;
+ for (i = 1; i <= atom->ntypes; i++)
+ for (j = i; j <= atom->ntypes; j++) {
+ fwrite(&setflag[i][j],sizeof(int),1,fp);
+ if (setflag[i][j]) {
+ fwrite(&a0[i][j],sizeof(double),1,fp);
+ fwrite(&sigma[i][j],sizeof(double),1,fp);
+ fwrite(&kappa[i][j],sizeof(double),1,fp);
+ fwrite(&cut[i][j],sizeof(double),1,fp);
+ }
+ }
+}
+
+/* ----------------------------------------------------------------------
+ proc 0 reads from restart file, bcasts
+------------------------------------------------------------------------- */
+
+void PairDPDfdtEnergy::read_restart(FILE *fp)
+{
+ read_restart_settings(fp);
+
+ allocate();
+
+ int i,j;
+ int me = comm->me;
+ for (i = 1; i <= atom->ntypes; i++)
+ for (j = i; j <= atom->ntypes; j++) {
+ if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp);
+ MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
+ if (setflag[i][j]) {
+ if (me == 0) {
+ fread(&a0[i][j],sizeof(double),1,fp);
+ fread(&sigma[i][j],sizeof(double),1,fp);
+ fread(&kappa[i][j],sizeof(double),1,fp);
+ fread(&cut[i][j],sizeof(double),1,fp);
+ }
+ MPI_Bcast(&a0[i][j],1,MPI_DOUBLE,0,world);
+ MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world);
+ MPI_Bcast(&kappa[i][j],1,MPI_DOUBLE,0,world);
+ MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world);
+ }
+ }
+}
+
+/* ----------------------------------------------------------------------
+ proc 0 writes to restart file
+------------------------------------------------------------------------- */
+
+void PairDPDfdtEnergy::write_restart_settings(FILE *fp)
+{
+ fwrite(&cut_global,sizeof(double),1,fp);
+ fwrite(&seed,sizeof(int),1,fp);
+ fwrite(&mix_flag,sizeof(int),1,fp);
+}
+
+/* ----------------------------------------------------------------------
+ proc 0 reads from restart file, bcasts
+------------------------------------------------------------------------- */
+
+void PairDPDfdtEnergy::read_restart_settings(FILE *fp)
+{
+ if (comm->me == 0) {
+ fread(&cut_global,sizeof(double),1,fp);
+ fread(&seed,sizeof(int),1,fp);
+ fread(&mix_flag,sizeof(int),1,fp);
+ }
+ MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world);
+ MPI_Bcast(&seed,1,MPI_INT,0,world);
+ MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
+
+ // initialize Marsaglia RNG with processor-unique seed
+ // same seed that pair_style command initially specified
+
+ if (random) delete random;
+ random = new RanMars(lmp,seed + comm->me);
+}
+
+/* ---------------------------------------------------------------------- */
+
+double PairDPDfdtEnergy::single(int i, int j, int itype, int jtype, double rsq,
+ double factor_coul, double factor_dpd, double &fforce)
+{
+ double r,rinv,wr,wd,phi;
+
+ r = sqrt(rsq);
+ if (r < EPSILON) {
+ fforce = 0.0;
+ return 0.0;
+ }
+
+ rinv = 1.0/r;
+ wr = 1.0 - r/cut[itype][jtype];
+ wd = wr*wr;
+ fforce = a0[itype][jtype]*wr * factor_dpd*rinv;
+
+ phi = 0.5*a0[itype][jtype]*cut[itype][jtype] * wd;
+ return factor_dpd*phi;
+}
+
diff --git a/src/USER-DPD/pair_dpd_fdt_energy.h b/src/USER-DPD/pair_dpd_fdt_energy.h
new file mode 100644
index 0000000000000000000000000000000000000000..4a68383398ca5e9886c6c096feb4bfb9fcd46f41
--- /dev/null
+++ b/src/USER-DPD/pair_dpd_fdt_energy.h
@@ -0,0 +1,90 @@
+/* ----------------------------------------------------------------------
+ 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: James Larentzos (Engility Corporation)
+------------------------------------------------------------------------- */
+
+#ifdef PAIR_CLASS
+
+PairStyle(dpd/fdt/energy,PairDPDfdtEnergy)
+
+#else
+
+#ifndef LMP_PAIR_DPD_FDT_ENERGY_H
+#define LMP_PAIR_DPD_FDT_ENERGY_H
+
+#include "pair.h"
+
+namespace LAMMPS_NS {
+
+class PairDPDfdtEnergy : public Pair {
+ public:
+ PairDPDfdtEnergy(class LAMMPS *);
+ virtual ~PairDPDfdtEnergy();
+ virtual void compute(int, int);
+ virtual void settings(int, char **);
+ virtual void coeff(int, char **);
+ void init_style();
+ double init_one(int, int);
+ void write_restart(FILE *);
+ void read_restart(FILE *);
+ virtual void write_restart_settings(FILE *);
+ virtual void read_restart_settings(FILE *);
+ double single(int, int, int, int, double, double, double, double &);
+
+ double **cut;
+ double **a0;
+ double **sigma,**kappa;
+
+ class RanMars *random;
+
+ protected:
+ double cut_global;
+ int seed;
+
+ void allocate();
+
+};
+
+}
+
+#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: Incorrect args for pair coefficients
+
+Self-explanatory. Check the input script or data file.
+
+E: Pair dpd/fdt/energy requires ghost atoms store velocity
+
+Use the communicate vel yes command to enable this.
+
+E: Pair dpd/fdt/energy requires newton pair on
+
+Self-explanatory.
+
+E: All pair coeffs are not set
+
+All pair coefficients must be set in the data file or by the
+pair_coeff command before running a simulation.
+
+*/