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

This commit is contained in:
sjplimp 2016-02-03 21:33:41 +00:00
parent 94da81bae4
commit de6ff01ba7
19 changed files with 4076 additions and 0 deletions

35
src/USER-DPD/README Normal file
View File

@ -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.

View File

@ -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;
}

View File

@ -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
*/

View File

@ -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;
}

View File

@ -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.
*/

View File

@ -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;
}

View File

@ -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
*/

108
src/USER-DPD/fix_eos_cv.cpp Normal file
View File

@ -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");
}
}

65
src/USER-DPD/fix_eos_cv.h Normal file
View File

@ -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.
*/

View File

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

View File

@ -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
*/

View File

@ -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++];
}
}
}

115
src/USER-DPD/fix_shardlow.h Normal file
View File

@ -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.
*/

View File

@ -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;
}

View File

@ -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.
*/

View File

@ -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;
}

View File

@ -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.
*/

View File

@ -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;
}

View File

@ -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.
*/