diff --git a/src/KOKKOS/verlet_kokkos.cpp b/src/KOKKOS/verlet_kokkos.cpp
index d063294023e1830ee0f54a5b63d470b027a3e714..f3c4e3e07d3f60d14d267a66008559bae22393c8 100644
--- a/src/KOKKOS/verlet_kokkos.cpp
+++ b/src/KOKKOS/verlet_kokkos.cpp
@@ -68,9 +68,10 @@ void VerletKokkos::setup()
{
if (comm->me == 0 && screen) {
fprintf(screen,"Setting up Verlet run ...\n");
- fprintf(screen," Unit style : %s\n", update->unit_style);
- fprintf(screen," Current step: " BIGINT_FORMAT "\n", update->ntimestep);
- fprintf(screen," Time step : %g\n", update->dt);
+ fprintf(screen," Unit style : %s\n", update->unit_style);
+ fprintf(screen," Current step : " BIGINT_FORMAT "\n", update->ntimestep);
+ fprintf(screen," Time step : %g\n", update->dt);
+ timer->print_timeout(screen);
}
update->setupflag = 1;
@@ -294,12 +295,16 @@ void VerletKokkos::run(int n)
f_merge_copy = DAT::t_f_array("VerletKokkos::f_merge_copy",atomKK->k_f.dimension_0());
static double time = 0.0;
- static int count = 0;
atomKK->sync(Device,ALL_MASK);
Kokkos::Impl::Timer ktimer;
+ timer->init_timeout();
for (int i = 0; i < n; i++) {
+ if (timer->check_timeout(i)) {
+ update->nsteps = i;
+ break;
+ }
ntimestep = ++update->ntimestep;
ev_set(ntimestep);
@@ -383,7 +388,6 @@ void VerletKokkos::run(int n)
unsigned int datamask_read_device = 0;
unsigned int datamask_modify_device = 0;
unsigned int datamask_read_host = 0;
- unsigned int datamask_modify_host = 0;
if ( pair_compute_flag ) {
if (force->pair->execution_space==Host) {
diff --git a/src/MANYBODY/fix_qeq_comb.cpp b/src/MANYBODY/fix_qeq_comb.cpp
index c2d9cd9281f44e4a59e73e2768b6bb70aeb65e7d..3c5954677b7c493cac3e966f8639e4d7dcfe4a9a 100644
--- a/src/MANYBODY/fix_qeq_comb.cpp
+++ b/src/MANYBODY/fix_qeq_comb.cpp
@@ -45,6 +45,8 @@ FixQEQComb::FixQEQComb(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
peratom_flag = 1;
size_peratom_cols = 0;
peratom_freq = 1;
+ respa_level_support = 1;
+ ilevel_respa = 0;
nevery = force->inumeric(FLERR,arg[3]);
precision = force->numeric(FLERR,arg[4]);
@@ -125,8 +127,10 @@ void FixQEQComb::init()
if (comb == NULL && comb3 == NULL)
error->all(FLERR,"Must use pair_style comb or comb3 with fix qeq/comb");
- if (strstr(update->integrate_style,"respa"))
- nlevels_respa = ((Respa *) update->integrate)->nlevels;
+ if (strstr(update->integrate_style,"respa")) {
+ ilevel_respa = ((Respa *) update->integrate)->nlevels-1;
+ if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa);
+ }
ngroup = group->count(igroup);
if (ngroup == 0) error->all(FLERR,"Fix qeq/comb group has no atoms");
@@ -140,9 +144,9 @@ void FixQEQComb::setup(int vflag)
if (strstr(update->integrate_style,"verlet"))
post_force(vflag);
else {
- ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
- post_force_respa(vflag,nlevels_respa-1,0);
- ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
+ ((Respa *) update->integrate)->copy_flevel_f(ilevel_respa);
+ post_force_respa(vflag,ilevel_respa,0);
+ ((Respa *) update->integrate)->copy_f_flevel(ilevel_respa);
}
firstflag = 0;
}
@@ -275,7 +279,7 @@ void FixQEQComb::post_force(int vflag)
void FixQEQComb::post_force_respa(int vflag, int ilevel, int iloop)
{
- if (ilevel == nlevels_respa-1) post_force(vflag);
+ if (ilevel == ilevel_respa) post_force(vflag);
}
/* ----------------------------------------------------------------------
diff --git a/src/MANYBODY/fix_qeq_comb.h b/src/MANYBODY/fix_qeq_comb.h
index d627ef669b1a79721e4c528bc242d0eca3dd5f32..384ea4bbccacd267300b1e8d36620f722affd0a5 100644
--- a/src/MANYBODY/fix_qeq_comb.h
+++ b/src/MANYBODY/fix_qeq_comb.h
@@ -43,7 +43,7 @@ class FixQEQComb : public Fix {
protected:
int me,firstflag;
double precision;
- int nlevels_respa;
+ int ilevel_respa;
bigint ngroup;
FILE *fp;
diff --git a/src/MISC/fix_efield.cpp b/src/MISC/fix_efield.cpp
index 2151a05092e8e8432f06670dc21e8e502fe4ba3c..85af46ed3b953d530ed9c6a3c33d57a9e1cd25c5 100644
--- a/src/MISC/fix_efield.cpp
+++ b/src/MISC/fix_efield.cpp
@@ -52,6 +52,8 @@ FixEfield::FixEfield(LAMMPS *lmp, int narg, char **arg) :
global_freq = 1;
extvector = 1;
extscalar = 1;
+ respa_level_support = 1;
+ ilevel_respa = 0;
qe2f = force->qe2f;
xstr = ystr = zstr = NULL;
@@ -216,8 +218,10 @@ void FixEfield::init()
update->whichflag == 2 && estyle == NONE)
error->all(FLERR,"Must use variable energy with fix efield");
- if (strstr(update->integrate_style,"respa"))
- nlevels_respa = ((Respa *) update->integrate)->nlevels;
+ if (strstr(update->integrate_style,"respa")) {
+ ilevel_respa = ((Respa *) update->integrate)->nlevels-1;
+ if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa);
+ }
}
/* ---------------------------------------------------------------------- */
@@ -227,9 +231,9 @@ void FixEfield::setup(int vflag)
if (strstr(update->integrate_style,"verlet"))
post_force(vflag);
else {
- ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
- post_force_respa(vflag,nlevels_respa-1,0);
- ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
+ ((Respa *) update->integrate)->copy_flevel_f(ilevel_respa);
+ post_force_respa(vflag,ilevel_respa,0);
+ ((Respa *) update->integrate)->copy_f_flevel(ilevel_respa);
}
}
@@ -393,7 +397,7 @@ void FixEfield::post_force(int vflag)
void FixEfield::post_force_respa(int vflag, int ilevel, int iloop)
{
- if (ilevel == nlevels_respa-1) post_force(vflag);
+ if (ilevel == ilevel_respa) post_force(vflag);
}
/* ---------------------------------------------------------------------- */
diff --git a/src/MISC/fix_efield.h b/src/MISC/fix_efield.h
index 93272ba9bfcb3dc9ecde6c76142c4cdf3ff472a0..4a360b6464c431ae07aa68550462bfed1e96979e 100644
--- a/src/MISC/fix_efield.h
+++ b/src/MISC/fix_efield.h
@@ -45,7 +45,7 @@ class FixEfield : public Fix {
char *xstr,*ystr,*zstr,*estr;
char *idregion;
int xvar,yvar,zvar,evar,xstyle,ystyle,zstyle,estyle;
- int nlevels_respa;
+ int ilevel_respa;
double qe2f;
int qflag,muflag;
diff --git a/src/MISC/fix_orient_bcc.cpp b/src/MISC/fix_orient_bcc.cpp
index 88db82bd82b016d1a0fbe79d1ef85269221570c9..24c587645d50ed3db8825518dd1bd674394c744e 100644
--- a/src/MISC/fix_orient_bcc.cpp
+++ b/src/MISC/fix_orient_bcc.cpp
@@ -72,6 +72,8 @@ FixOrientBCC::FixOrientBCC(LAMMPS *lmp, int narg, char **arg) :
peratom_flag = 1;
size_peratom_cols = 2;
peratom_freq = 1;
+ respa_level_support = 1;
+ ilevel_respa = 0;
nstats = force->inumeric(FLERR,arg[3]);
direction_of_motion = force->inumeric(FLERR,arg[4]);
@@ -210,8 +212,10 @@ int FixOrientBCC::setmask()
void FixOrientBCC::init()
{
- if (strstr(update->integrate_style,"respa"))
- nlevels_respa = ((Respa *) update->integrate)->nlevels;
+ if (strstr(update->integrate_style,"respa")) {
+ ilevel_respa = ((Respa *) update->integrate)->nlevels-1;
+ if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa);
+ }
// need a full neighbor list, built whenever re-neighboring occurs
@@ -236,9 +240,9 @@ void FixOrientBCC::setup(int vflag)
if (strstr(update->integrate_style,"verlet"))
post_force(vflag);
else {
- ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
- post_force_respa(vflag,nlevels_respa-1,0);
- ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
+ ((Respa *) update->integrate)->copy_flevel_f(ilevel_respa);
+ post_force_respa(vflag,ilevel_respa,0);
+ ((Respa *) update->integrate)->copy_f_flevel(ilevel_respa);
}
}
@@ -467,7 +471,7 @@ void FixOrientBCC::post_force(int vflag)
void FixOrientBCC::post_force_respa(int vflag, int ilevel, int iloop)
{
- if (ilevel == nlevels_respa-1) post_force(vflag);
+ if (ilevel == ilevel_respa) post_force(vflag);
}
/* ---------------------------------------------------------------------- */
diff --git a/src/MISC/fix_orient_bcc.h b/src/MISC/fix_orient_bcc.h
index aecbaca793bd6ee20bbbed7d1a5f3f9981148b02..2035e85b155c68f3ed03cd763ec583a3c94eb05d 100644
--- a/src/MISC/fix_orient_bcc.h
+++ b/src/MISC/fix_orient_bcc.h
@@ -58,7 +58,7 @@ class FixOrientBCC : public Fix {
private:
int me;
- int nlevels_respa;
+ int ilevel_respa;
int direction_of_motion; // 1 = center shrinks, 0 = center grows
int nstats; // stats output every this many steps
diff --git a/src/MISC/fix_orient_fcc.cpp b/src/MISC/fix_orient_fcc.cpp
index df41981f3d966a910ae560f9c260283f6bf31505..c0c836dca7c7513dc53e64216be38413f204a2c2 100644
--- a/src/MISC/fix_orient_fcc.cpp
+++ b/src/MISC/fix_orient_fcc.cpp
@@ -69,6 +69,8 @@ FixOrientFCC::FixOrientFCC(LAMMPS *lmp, int narg, char **arg) :
peratom_flag = 1;
size_peratom_cols = 2;
peratom_freq = 1;
+ respa_level_support = 1;
+ ilevel_respa = 0;
nstats = force->inumeric(FLERR,arg[3]);
direction_of_motion = force->inumeric(FLERR,arg[4]);
@@ -207,8 +209,10 @@ int FixOrientFCC::setmask()
void FixOrientFCC::init()
{
- if (strstr(update->integrate_style,"respa"))
- nlevels_respa = ((Respa *) update->integrate)->nlevels;
+ if (strstr(update->integrate_style,"respa")) {
+ ilevel_respa = ((Respa *) update->integrate)->nlevels-1;
+ if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa);
+ }
// need a full neighbor list, built whenever re-neighboring occurs
@@ -233,9 +237,9 @@ void FixOrientFCC::setup(int vflag)
if (strstr(update->integrate_style,"verlet"))
post_force(vflag);
else {
- ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
- post_force_respa(vflag,nlevels_respa-1,0);
- ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
+ ((Respa *) update->integrate)->copy_flevel_f(ilevel_respa);
+ post_force_respa(vflag,ilevel_respa,0);
+ ((Respa *) update->integrate)->copy_f_flevel(ilevel_respa);
}
}
@@ -464,7 +468,7 @@ void FixOrientFCC::post_force(int vflag)
void FixOrientFCC::post_force_respa(int vflag, int ilevel, int iloop)
{
- if (ilevel == nlevels_respa-1) post_force(vflag);
+ if (ilevel == ilevel_respa) post_force(vflag);
}
/* ---------------------------------------------------------------------- */
diff --git a/src/MISC/fix_orient_fcc.h b/src/MISC/fix_orient_fcc.h
index 3d5d1207aa0239586cefcc50c84b07036b9bc9b2..8bdc098289cb68b7a8987da6e78d82d16aeb89db 100644
--- a/src/MISC/fix_orient_fcc.h
+++ b/src/MISC/fix_orient_fcc.h
@@ -58,7 +58,7 @@ class FixOrientFCC : public Fix {
private:
int me;
- int nlevels_respa;
+ int ilevel_respa;
int direction_of_motion; // 1 = center shrinks, 0 = center grows
int nstats; // stats output every this many steps
diff --git a/src/RIGID/fix_ehex.cpp b/src/RIGID/fix_ehex.cpp
index 3d39bba93d07200fd74719562843f68e68bdfe4e..c9688873747d40f425c06a0ca5cdb87c2ec5c975 100644
--- a/src/RIGID/fix_ehex.cpp
+++ b/src/RIGID/fix_ehex.cpp
@@ -244,13 +244,11 @@ void FixEHEX::end_of_step() {
------------------------------------------------------------------------- */
void FixEHEX::rescale() {
- double heat,ke,massone;
double Kr, Ke, escale;
double vsub[3],vcm[3], sfr[3];
- double dvcmdt[3];
double mi;
double dt;
- double F, mr, Fo2Kr, epsr_ik, sfvr, eta_ik;
+ double F, mr, epsr_ik, sfvr, eta_ik;
dt = update->dt;
@@ -266,8 +264,6 @@ void FixEHEX::rescale() {
mr = masstotal;
- Fo2Kr = F / (2.*Kr);
-
// energy scaling factor
escale = 1. + (F*dt)/Kr;
@@ -406,7 +402,7 @@ void FixEHEX::update_scalingmask() {
inside the region.
------------------------------------------------------------------------- */
-bool FixEHEX::check_cluster(int *shake_atom, int n, Region * region) {
+bool FixEHEX::check_cluster(tagint *shake_atom, int n, Region * region) {
// IMPORTANT NOTE: If any site of the cluster belongs to a group
// which should not be rescaled than all of the sites
diff --git a/src/RIGID/fix_ehex.h b/src/RIGID/fix_ehex.h
index 1e7a0fe482fa2572680945eaffa252cdd635cc6f..3ae974fc0dc9bbae1d50d7b57a80e7e0bd49c187 100644
--- a/src/RIGID/fix_ehex.h
+++ b/src/RIGID/fix_ehex.h
@@ -43,7 +43,7 @@ class FixEHEX : public Fix {
void com_properties(double *, double *, double *, double*, double *, double*);
bool rescale_atom(int i, Region*region);
virtual void grow_arrays(int nmax);
- bool check_cluster(int *shake_atom, int n, Region * region);
+ bool check_cluster(tagint *shake_atom, int n, Region * region);
private:
int iregion;
diff --git a/src/RIGID/fix_rattle.cpp b/src/RIGID/fix_rattle.cpp
index 3a5fef69cb2f9916fcc37bbc130c805f6fdfdb05..359be6b885e2ae3085d6b1f68bde7ce349fffe9f 100644
--- a/src/RIGID/fix_rattle.cpp
+++ b/src/RIGID/fix_rattle.cpp
@@ -240,7 +240,6 @@ void FixRattle::final_integrate_respa(int ilevel, int iloop)
void FixRattle::vrattle3angle(int m)
{
tagint i0,i1,i2;
- int nlist,list[3];
double c[3], l[3], a[3][3], r01[3], imass[3],
r02[3], r12[3], vp01[3], vp02[3], vp12[3];
@@ -323,7 +322,6 @@ void FixRattle::vrattle3angle(int m)
void FixRattle::vrattle2(int m)
{
tagint i0, i1;
- int nlist,list[2];
double imass[2], r01[3], vp01[3];
// local atom IDs and constraint distances
@@ -372,7 +370,6 @@ void FixRattle::vrattle2(int m)
void FixRattle::vrattle3(int m)
{
tagint i0,i1,i2;
- int nlist,list[3];
double imass[3], r01[3], r02[3], vp01[3], vp02[3],
a[2][2],c[2],l[2];
@@ -442,7 +439,6 @@ void FixRattle::vrattle3(int m)
void FixRattle::vrattle4(int m)
{
tagint i0,i1,i2,i3;
- int nlist,list[4];
double imass[4], c[3], l[3], a[3][3],
r01[3], r02[3], r03[3], vp01[3], vp02[3], vp03[3];
@@ -592,7 +588,6 @@ void FixRattle::solve3x3exactly(const double a[][3],
void FixRattle::reset_dt()
{
FixShake::reset_dt();
- dtfv = 0.5 * update->dt * force->ftm2v;
}
/* ----------------------------------------------------------------------
@@ -601,8 +596,8 @@ void FixRattle::reset_dt()
void FixRattle::update_v_half_nocons()
{
+ const double dtfv = 0.5 * update->dt * force->ftm2v;
double dtfvinvm;
- dtfv = 0.5 * update->dt * force->ftm2v;
if (rmass) {
for (int i = 0; i < nlocal; i++) {
@@ -632,10 +627,6 @@ void FixRattle::update_v_half_nocons()
void FixRattle::update_v_half_nocons_respa(int ilevel)
{
- // select timestep for current level
-
- dtfv = 0.5 * step_respa[ilevel] * force->ftm2v;
-
// carry out unconstrained velocity update
update_v_half_nocons();
diff --git a/src/RIGID/fix_rattle.h b/src/RIGID/fix_rattle.h
index 7b94d2cf0d63451004ff7723a8463be86822c87c..e45dfb1a4aba6639c976f6101a514558833d7cd1 100644
--- a/src/RIGID/fix_rattle.h
+++ b/src/RIGID/fix_rattle.h
@@ -28,7 +28,6 @@ namespace LAMMPS_NS {
class FixRattle : public FixShake {
public:
double **vp; // array for unconstrained velocities
- double dtfv; // timestep for velocity update
int comm_mode; // mode for communication pack/unpack
double derr_max; // distance error
double verr_max; // velocity error
diff --git a/src/USER-DPD/atom_vec_dpd.cpp b/src/USER-DPD/atom_vec_dpd.cpp
index 92028b60e621e1bfc51a3b252720ed7c86b904ea..586ddca02b396eef10c247f0012ed8bd28b539c1 100644
--- a/src/USER-DPD/atom_vec_dpd.cpp
+++ b/src/USER-DPD/atom_vec_dpd.cpp
@@ -76,9 +76,8 @@ void AtomVecDPD::grow(int n)
uChem = memory->grow(atom->uChem,nmax,"atom:uChem");
uCG = memory->grow(atom->uCG,nmax,"atom:uCG");
uCGnew = memory->grow(atom->uCGnew,nmax,"atom:uCGnew");
- duCond = memory->grow(atom->duCond,nmax,"atom:duCond");
- duMech = memory->grow(atom->duMech,nmax,"atom:duMech");
duChem = memory->grow(atom->duChem,nmax,"atom:duChem");
+ ssaAIR = memory->grow(atom->ssaAIR,nmax,"atom:ssaAIR");
if (atom->nextra_grow)
for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
@@ -101,9 +100,8 @@ void AtomVecDPD::grow_reset()
uChem = atom->uChem;
uCG = atom->uCG;
uCGnew = atom->uCGnew;
- duCond = atom->duCond;
- duMech = atom->duMech;
duChem = atom->duChem;
+ ssaAIR = atom->ssaAIR;
}
/* ----------------------------------------------------------------------
@@ -128,6 +126,7 @@ void AtomVecDPD::copy(int i, int j, int delflag)
uChem[j] = uChem[i];
uCG[j] = uCG[i];
uCGnew[j] = uCGnew[i];
+ ssaAIR[j] = ssaAIR[i];
if (atom->nextra_grow)
for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
@@ -220,10 +219,10 @@ int AtomVecDPD::pack_comm_vel(int n, int *list, double *buf,
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];
- buf[m++] = uChem[j];
+ buf[m++] = dpdTheta[j];
+ buf[m++] = uCond[j];
+ buf[m++] = uMech[j];
+ buf[m++] = uChem[j];
}
} else {
dvx = pbc[0]*h_rate[0] + pbc[5]*h_rate[5] + pbc[4]*h_rate[4];
@@ -243,10 +242,10 @@ int AtomVecDPD::pack_comm_vel(int n, int *list, double *buf,
buf[m++] = v[j][1];
buf[m++] = v[j][2];
}
- buf[m++] = dpdTheta[j];
- buf[m++] = uCond[j];
- buf[m++] = uMech[j];
- buf[m++] = uChem[j];
+ buf[m++] = dpdTheta[j];
+ buf[m++] = uCond[j];
+ buf[m++] = uMech[j];
+ buf[m++] = uChem[j];
}
}
}
@@ -428,18 +427,18 @@ int AtomVecDPD::pack_border_vel(int n, int *list, double *buf,
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++] = 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];
- buf[m++] = uChem[j];
- buf[m++] = uCG[j];
- buf[m++] = uCGnew[j];
+ buf[m++] = dpdTheta[j];
+ buf[m++] = uCond[j];
+ buf[m++] = uMech[j];
+ buf[m++] = uChem[j];
+ buf[m++] = uCG[j];
+ buf[m++] = uCGnew[j];
}
} else {
dvx = pbc[0]*h_rate[0] + pbc[5]*h_rate[5] + pbc[4]*h_rate[4];
@@ -450,9 +449,9 @@ int AtomVecDPD::pack_border_vel(int n, int *list, double *buf,
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++] = 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;
@@ -462,12 +461,12 @@ int AtomVecDPD::pack_border_vel(int n, int *list, double *buf,
buf[m++] = v[j][1];
buf[m++] = v[j][2];
}
- buf[m++] = dpdTheta[j];
- buf[m++] = uCond[j];
- buf[m++] = uMech[j];
- buf[m++] = uChem[j];
- buf[m++] = uCG[j];
- buf[m++] = uCGnew[j];
+ buf[m++] = dpdTheta[j];
+ buf[m++] = uCond[j];
+ buf[m++] = uMech[j];
+ buf[m++] = uChem[j];
+ buf[m++] = uCG[j];
+ buf[m++] = uCGnew[j];
}
}
}
@@ -517,6 +516,41 @@ int AtomVecDPD::pack_border_hybrid(int n, int *list, double *buf)
return m;
}
+/* ----------------------------------------------------------------------
+ convert atom coords into the ssa active interaction region number
+------------------------------------------------------------------------- */
+
+int AtomVecDPD::coord2ssaAIR(double *x)
+{
+ int ix, iy, iz;
+
+ ix = iy = iz = 0;
+ if (x[2] < domain->sublo[2]) iz = -1;
+ if (x[2] > domain->subhi[2]) iz = 1;
+ if (x[1] < domain->sublo[1]) iy = -1;
+ if (x[1] > domain->subhi[1]) iy = 1;
+ if (x[0] < domain->sublo[0]) ix = -1;
+ if (x[0] > domain->subhi[0]) ix = 1;
+
+ if(iz < 0){
+ return 0;
+ } else if(iz == 0){
+ if( iy<0 ) return 0; // bottom left/middle/right
+ if( (iy==0) && (ix<0) ) return 0; // left atoms
+ if( (iy==0) && (ix==0) ) return 1; // Locally owned atoms
+ if( (iy==0) && (ix>0) ) return 3; // Right atoms
+ if( (iy>0) && (ix==0) ) return 2; // Top-middle atoms
+ if( (iy>0) && (ix!=0) ) return 4; // Top-right and top-left atoms
+ } else { // iz > 0
+ if((ix==0) && (iy==0)) return 5; // Back atoms
+ if((ix==0) && (iy!=0)) return 6; // Top-back and bottom-back atoms
+ if((ix!=0) && (iy==0)) return 7; // Left-back and right-back atoms
+ if((ix!=0) && (iy!=0)) return 8; // Back corner atoms
+ }
+
+ return 0;
+}
+
/* ---------------------------------------------------------------------- */
void AtomVecDPD::unpack_border(int n, int first, double *buf)
@@ -539,6 +573,7 @@ void AtomVecDPD::unpack_border(int n, int first, double *buf)
uChem[i] = buf[m++];
uCG[i] = buf[m++];
uCGnew[i] = buf[m++];
+ ssaAIR[i] = coord2ssaAIR(x[i]);
}
if (atom->nextra_border)
@@ -572,6 +607,7 @@ void AtomVecDPD::unpack_border_vel(int n, int first, double *buf)
uChem[i] = buf[m++];
uCG[i] = buf[m++];
uCGnew[i] = buf[m++];
+ ssaAIR[i] = coord2ssaAIR(x[i]);
}
if (atom->nextra_border)
@@ -612,6 +648,7 @@ int AtomVecDPD::unpack_border_hybrid(int n, int first, double *buf)
uChem[i] = buf[m++];
uCG[i] = buf[m++];
uCGnew[i] = buf[m++];
+ ssaAIR[i] = coord2ssaAIR(x[i]);
}
return m;
}
@@ -673,6 +710,7 @@ int AtomVecDPD::unpack_exchange(double *buf)
uChem[nlocal] = buf[m++];
uCG[nlocal] = buf[m++];
uCGnew[nlocal] = buf[m++];
+ ssaAIR[nlocal] = 1; /* coord2ssaAIR(x[nlocal]) */
if (atom->nextra_grow)
for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
@@ -763,6 +801,7 @@ int AtomVecDPD::unpack_restart(double *buf)
uCond[nlocal] = buf[m++];
uMech[nlocal] = buf[m++];
uChem[nlocal] = buf[m++];
+ ssaAIR[nlocal] = 1; /* coord2ssaAIR(x[nlocal]) */
double **extra = atom->extra;
if (atom->nextra_store) {
@@ -802,9 +841,8 @@ void AtomVecDPD::create_atom(int itype, double *coord)
uChem[nlocal] = 0.0;
uCG[nlocal] = 0.0;
uCGnew[nlocal] = 0.0;
- duCond[nlocal] = 0.0;
- duMech[nlocal] = 0.0;
duChem[nlocal] = 0.0;
+ ssaAIR[nlocal] = 1; /* coord2ssaAIR(x[nlocal]) */
atom->nlocal++;
}
@@ -845,6 +883,7 @@ void AtomVecDPD::data_atom(double *coord, tagint imagetmp, char **values)
uChem[nlocal] = 0.0;
uCG[nlocal] = 0.0;
uCGnew[nlocal] = 0.0;
+ ssaAIR[nlocal] = 1; /* coord2ssaAIR(x[nlocal]) */
atom->nlocal++;
}
@@ -857,6 +896,7 @@ void AtomVecDPD::data_atom(double *coord, tagint imagetmp, char **values)
int AtomVecDPD::data_atom_hybrid(int nlocal, char **values)
{
dpdTheta[nlocal] = atof(values[0]);
+ ssaAIR[nlocal] = 1; /* coord2ssaAIR(x[nlocal]) */
return 1;
}
@@ -900,9 +940,9 @@ 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],
+ 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);
+ (int) ubuf(buf[i][8]).i);
}
/* ----------------------------------------------------------------------
@@ -937,9 +977,8 @@ bigint AtomVecDPD::memory_usage()
if (atom->memcheck("uChem")) bytes += memory->usage(uChem,nmax);
if (atom->memcheck("uCG")) bytes += memory->usage(uCG,nmax);
if (atom->memcheck("uCGnew")) bytes += memory->usage(uCGnew,nmax);
- if (atom->memcheck("duCond")) bytes += memory->usage(duCond,nmax);
- if (atom->memcheck("duMech")) bytes += memory->usage(duMech,nmax);
if (atom->memcheck("duChem")) bytes += memory->usage(duChem,nmax);
+ if (atom->memcheck("ssaAIR")) bytes += memory->usage(ssaAIR,nmax);
return bytes;
}
diff --git a/src/USER-DPD/atom_vec_dpd.h b/src/USER-DPD/atom_vec_dpd.h
index 71ef48d053180f2bd8afff21792fc025367c6c7c..077d18485f4fb7d86ef90d7ea796ddbb645a4b6b 100644
--- a/src/USER-DPD/atom_vec_dpd.h
+++ b/src/USER-DPD/atom_vec_dpd.h
@@ -1,4 +1,4 @@
-/* ----------------------------------------------------------------------
+/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
@@ -59,13 +59,16 @@ class AtomVecDPD : public AtomVec {
int write_data_hybrid(FILE *, double *);
bigint memory_usage();
double *uCond,*uMech,*uChem,*uCG,*uCGnew,*rho,*dpdTheta;
- double *duCond,*duMech,*duChem;
+ double *duChem;
+ int *ssaAIR; // Shardlow Splitting Algorithm Active Interaction Region number
protected:
tagint *tag;
int *type,*mask;
imageint *image;
double **x,**v,**f;
+
+ int coord2ssaAIR(double *); // map atom coord to an AIR number
};
}
diff --git a/src/USER-DPD/compute_dpd.cpp b/src/USER-DPD/compute_dpd.cpp
index 59716362c773e74b7947d888e6591031f39b5224..08aa19bcdfcffed106ff2d1836473840a4546525 100644
--- a/src/USER-DPD/compute_dpd.cpp
+++ b/src/USER-DPD/compute_dpd.cpp
@@ -36,9 +36,9 @@ ComputeDpd::ComputeDpd(LAMMPS *lmp, int narg, char **arg) :
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)");
}
@@ -67,7 +67,7 @@ void ComputeDpd::compute_vector()
dpdU = new double[size_vector];
- for (int i = 0; i < size_vector; i++) dpdU[i] = double(0.0);
+ for (int i = 0; i < size_vector; i++) dpdU[i] = 0.0;
for (int i = 0; i < nlocal; i++){
if (mask[i] & groupbit){
@@ -78,7 +78,7 @@ void ComputeDpd::compute_vector()
dpdU[4]++;
}
}
-
+
MPI_Allreduce(dpdU,vector,size_vector,MPI_DOUBLE,MPI_SUM,world);
natoms = vector[4];
diff --git a/src/USER-DPD/compute_dpd.h b/src/USER-DPD/compute_dpd.h
index c2e8133cd48db4aeb59317aa687efa7502009916..695048eef38ebb3a750fa089afe08474523550ab 100644
--- a/src/USER-DPD/compute_dpd.h
+++ b/src/USER-DPD/compute_dpd.h
@@ -1,4 +1,4 @@
-/* ----------------------------------------------------------------------
+/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
@@ -51,6 +51,6 @@ 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.
+Self-explanatory.
*/
diff --git a/src/USER-DPD/compute_dpd_atom.h b/src/USER-DPD/compute_dpd_atom.h
index e591e4c48e0ae233d6113cf3d8fbc53f50bc4c16..c77490be764010af3c3df9e8e2c04cce12b42de0 100644
--- a/src/USER-DPD/compute_dpd_atom.h
+++ b/src/USER-DPD/compute_dpd_atom.h
@@ -1,4 +1,4 @@
-/* ----------------------------------------------------------------------
+/* -*- c++ -*- ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
diff --git a/src/USER-DPD/fix_eos_cv.cpp b/src/USER-DPD/fix_eos_cv.cpp
index 4a8d213362e564ed76d0bdf56b955504f17ddd4c..df764f708295d7ad8af729879e54dd1fcf6ae848 100644
--- a/src/USER-DPD/fix_eos_cv.cpp
+++ b/src/USER-DPD/fix_eos_cv.cpp
@@ -32,7 +32,7 @@ FixEOScv::FixEOScv(LAMMPS *lmp, int narg, char **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");
+ if(cvEOS <= 0.0) error->all(FLERR,"EOS cv must be > 0.0");
restart_peratom = 1;
nevery = 1;
@@ -61,14 +61,14 @@ void FixEOScv::init()
if(this->restart_reset){
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
- dpdTheta[i] = (uCond[i]+uMech[i])/cvEOS;
+ dpdTheta[i] = (uCond[i]+uMech[i])/cvEOS;
} else {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
- if(dpdTheta[i] <= double(0.0))
- error->one(FLERR,"Internal temperature <= zero");
- uCond[i] = double(0.5)*cvEOS*dpdTheta[i];
- uMech[i] = double(0.5)*cvEOS*dpdTheta[i];
+ if(dpdTheta[i] <= 0.0)
+ error->one(FLERR,"Internal temperature <= zero");
+ uCond[i] = 0.5*cvEOS*dpdTheta[i];
+ uMech[i] = 0.5*cvEOS*dpdTheta[i];
}
}
}
@@ -86,8 +86,8 @@ void FixEOScv::post_integrate()
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");
+ if(dpdTheta[i] <= 0.0)
+ error->one(FLERR,"Internal temperature <= zero");
}
}
@@ -104,7 +104,7 @@ void FixEOScv::end_of_step()
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");
+ if(dpdTheta[i] <= 0.0)
+ error->one(FLERR,"Internal temperature <= zero");
}
}
diff --git a/src/USER-DPD/fix_eos_cv.h b/src/USER-DPD/fix_eos_cv.h
index 043b8b82d1361742a43af8a2e7c303756ce90c62..f53c3177af93a0f9e469e7d5fb7a0d0fe17ee737 100644
--- a/src/USER-DPD/fix_eos_cv.h
+++ b/src/USER-DPD/fix_eos_cv.h
@@ -1,4 +1,4 @@
-/* ----------------------------------------------------------------------
+/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
diff --git a/src/USER-DPD/fix_eos_table.cpp b/src/USER-DPD/fix_eos_table.cpp
index b9579cfaf7e842f1096e7fe643129fe3a9808876..98276d0f018630ff2d82607c00e3352d415e1c01 100644
--- a/src/USER-DPD/fix_eos_table.cpp
+++ b/src/USER-DPD/fix_eos_table.cpp
@@ -112,15 +112,15 @@ void FixEOStable::init()
if(this->restart_reset){
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
- temperature_lookup(uCond[i]+uMech[i],dpdTheta[i]);
+ temperature_lookup(uCond[i]+uMech[i],dpdTheta[i]);
} else {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
- if(dpdTheta[i] <= double(0.0))
- error->one(FLERR,"Internal temperature <= zero");
- energy_lookup(dpdTheta[i],tmp);
- uCond[i] = tmp / double(2.0);
- uMech[i] = tmp / double(2.0);
+ if(dpdTheta[i] <= 0.0)
+ error->one(FLERR,"Internal temperature <= zero");
+ energy_lookup(dpdTheta[i],tmp);
+ uCond[i] = tmp / 2.0;
+ uMech[i] = tmp / 2.0;
}
}
}
@@ -138,8 +138,8 @@ void FixEOStable::post_integrate()
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");
+ if(dpdTheta[i] <= 0.0)
+ error->one(FLERR,"Internal temperature <= zero");
}
}
@@ -156,8 +156,8 @@ void FixEOStable::end_of_step()
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");
+ if(dpdTheta[i] <= 0.0)
+ error->one(FLERR,"Internal temperature <= zero");
}
}
@@ -195,13 +195,13 @@ void FixEOStable::read_table(Table *tb, Table *tb2, char *file, char *keyword)
// open file
- FILE *fp = fopen(file,"r");
+ FILE *fp = force->open_potential(file);
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) {
@@ -426,7 +426,7 @@ void FixEOStable::temperature_lookup(double u, double &t)
double fraction;
Table *tb = &tables[1];
- if(u < tb->lo || u > tb->hi){
+ 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");
}
diff --git a/src/USER-DPD/fix_eos_table.h b/src/USER-DPD/fix_eos_table.h
index 7c33cf6059d31051780791b09198d22269866aa1..1623fa21f8b4d39b5247ee18dcec633699b642ec 100644
--- a/src/USER-DPD/fix_eos_table.h
+++ b/src/USER-DPD/fix_eos_table.h
@@ -1,4 +1,4 @@
-/* ----------------------------------------------------------------------
+/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
@@ -99,7 +99,7 @@ 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.
+correct.
E: Did not find keyword in table file
@@ -115,12 +115,12 @@ 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
+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
+The internal energy does not lie with the minimum
and maximum energy cutoffs of the table
*/
diff --git a/src/USER-DPD/fix_eos_table_rx.cpp b/src/USER-DPD/fix_eos_table_rx.cpp
index f7acada9d562271aff90bc398680ce2a2ae88aa7..5f382436e87a4a1abf4fb020b8a12e59ed96ba4c 100644
--- a/src/USER-DPD/fix_eos_table_rx.cpp
+++ b/src/USER-DPD/fix_eos_table_rx.cpp
@@ -57,6 +57,8 @@ FixEOStableRX::FixEOStableRX(LAMMPS *lmp, int narg, char **arg) :
ntables = 0;
tables = NULL;
tables2 = NULL;
+ eosSpecies = NULL;
+
int me;
MPI_Comm_rank(world,&me);
for (int ii=0;ii<nspecies;ii++){
@@ -70,7 +72,7 @@ FixEOStableRX::FixEOStableRX(LAMMPS *lmp, int narg, char **arg) :
null_table(tb);
null_table(tb2);
-
+
ntables++;
}
@@ -79,7 +81,7 @@ FixEOStableRX::FixEOStableRX(LAMMPS *lmp, int narg, char **arg) :
Table *tb2 = &tables2[ntables];
if (me == 0) read_table(tb,tb2,arg[4],arg[6]);
-
+
for (int ii=0;ii<nspecies;ii++){
Table *tb = &tables[ntables];
Table *tb2 = &tables2[ntables];
@@ -88,21 +90,21 @@ FixEOStableRX::FixEOStableRX(LAMMPS *lmp, int narg, char **arg) :
bcast_table(tb2);
// error check on table parameters
-
+
if (tb->ninput <= 1) error->one(FLERR,"Invalid eos/table/rx length");
-
+
tb->lo = tb->rfile[0];
tb->hi = tb->rfile[tb->ninput-1];
if (tb->lo >= tb->hi) error->all(FLERR,"eos/table/rx values are not increasing");
-
+
if (tb2->ninput <= 1) error->one(FLERR,"Invalid eos/table/rx length");
-
+
tb2->lo = tb2->rfile[0];
tb2->hi = tb2->rfile[tb2->ninput-1];
if (tb2->lo >= tb2->hi) error->all(FLERR,"eos/table/rx values are not increasing");
-
+
// spline read-in and compute r,e,f vectors within table
-
+
spline_table(tb);
compute_table(tb);
spline_table(tb2);
@@ -129,6 +131,7 @@ FixEOStableRX::~FixEOStableRX()
memory->sfree(tables2);
delete [] dHf;
+ delete [] eosSpecies;
}
/* ---------------------------------------------------------------------- */
@@ -159,8 +162,8 @@ void FixEOStableRX::setup(int vflag)
if (mask[i] & groupbit){
duChem = uCG[i] - uCGnew[i];
uChem[i] += duChem;
- uCG[i] = double(0.0);
- uCGnew[i] = double(0.0);
+ uCG[i] = 0.0;
+ uCGnew[i] = 0.0;
}
// Communicate the updated momenta and velocities to all nodes
@@ -182,20 +185,20 @@ void FixEOStableRX::init()
double *uChem = atom->uChem;
double *dpdTheta = atom->dpdTheta;
double tmp;
-
+
if(this->restart_reset){
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
- temperature_lookup(i,uCond[i]+uMech[i]+uChem[i],dpdTheta[i]);
+ temperature_lookup(i,uCond[i]+uMech[i]+uChem[i],dpdTheta[i]);
} else {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
- if(dpdTheta[i] <= double(0.0))
- error->one(FLERR,"Internal temperature <= zero");
- energy_lookup(i,dpdTheta[i],tmp);
- uCond[i] = tmp / double(2.0);
- uMech[i] = tmp / double(2.0);
- uChem[i] = double(0.0);
+ if(dpdTheta[i] <= 0.0)
+ error->one(FLERR,"Internal temperature <= zero");
+ energy_lookup(i,dpdTheta[i],tmp);
+ uCond[i] = tmp / 2.0;
+ uMech[i] = tmp / 2.0;
+ uChem[i] = 0.0;
}
}
}
@@ -215,8 +218,8 @@ void FixEOStableRX::post_integrate()
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit){
temperature_lookup(i,uCond[i]+uMech[i]+uChem[i],dpdTheta[i]);
- if(dpdTheta[i] <= double(0.0))
- error->one(FLERR,"Internal temperature <= zero");
+ if(dpdTheta[i] <= 0.0)
+ error->one(FLERR,"Internal temperature <= zero");
}
}
@@ -241,8 +244,8 @@ void FixEOStableRX::end_of_step()
if (mask[i] & groupbit){
duChem = uCG[i] - uCGnew[i];
uChem[i] += duChem;
- uCG[i] = double(0.0);
- uCGnew[i] = double(0.0);
+ uCG[i] = 0.0;
+ uCGnew[i] = 0.0;
}
// Communicate the updated momenta and velocities to all nodes
@@ -251,8 +254,8 @@ void FixEOStableRX::end_of_step()
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit){
temperature_lookup(i,uCond[i]+uMech[i]+uChem[i],dpdTheta[i]);
- if(dpdTheta[i] <= double(0.0))
- error->one(FLERR,"Internal temperature <= zero");
+ if(dpdTheta[i] <= 0.0)
+ error->one(FLERR,"Internal temperature <= zero");
}
}
@@ -332,9 +335,9 @@ void FixEOStableRX::read_file(char *file)
for (ispecies = 0; ispecies < nspecies; ispecies++)
if (strcmp(words[0],&atom->dname[ispecies][0]) == 0) break;
- if (ispecies == nspecies) continue;
- dHf[ispecies] = atof(words[1]);
+ if (ispecies < nspecies)
+ dHf[ispecies] = atof(words[1]);
}
delete [] words;
@@ -380,7 +383,7 @@ void FixEOStableRX::read_table(Table *tb, Table *tb2, char *file, char *keyword)
sprintf(str,"Cannot open file %s",file);
error->one(FLERR,str);
}
-
+
// loop until section found with matching keyword
while (1) {
@@ -406,13 +409,13 @@ void FixEOStableRX::read_table(Table *tb, Table *tb2, char *file, char *keyword)
memory->create(tb->efile,tb->ninput,"eos:efile");
memory->create(tb2->rfile,tb2->ninput,"eos:rfile");
memory->create(tb2->efile,tb2->ninput,"eos:efile");
-
+
for (int ispecies=1;ispecies<nspecies;ispecies++){
Table *tbl = &tables[ispecies];
Table *tbl2 = &tables2[ispecies];
tbl->ninput = tb->ninput;
tbl2->ninput = tb2->ninput;
-
+
memory->create(tbl->rfile,tbl->ninput,"eos:rfile");
memory->create(tbl->efile,tbl->ninput,"eos:efile");
memory->create(tbl2->rfile,tbl2->ninput,"eos:rfile");
@@ -421,7 +424,6 @@ void FixEOStableRX::read_table(Table *tb, Table *tb2, char *file, char *keyword)
// read r,e table values from file
- int itmp;
double rtmp, tmpE;
int nwords;
char * word;
@@ -431,7 +433,7 @@ void FixEOStableRX::read_table(Table *tb, Table *tb2, char *file, char *keyword)
fgets(line,MAXLINE,fp);
for (int i = 0; i < ninputs; i++) {
fgets(line,MAXLINE,fp);
-
+
nwords = atom->count_words(line);
if(nwords != nspecies+2){
printf("nwords=%d nspecies=%d\n",nwords,nspecies);
@@ -439,7 +441,6 @@ void FixEOStableRX::read_table(Table *tb, Table *tb2, char *file, char *keyword)
}
nwords = 0;
word = strtok(line," \t\n\r\f");
- itmp = atoi(word);
word = strtok(NULL," \t\n\r\f");
rtmp = atof(word);
@@ -523,10 +524,11 @@ void FixEOStableRX::param_extract(Table *tb, char *line)
int ispecies;
ncolumn = 0;
- eosSpecies = new int[nspecies];
+ if (!eosSpecies)
+ eosSpecies = new int[nspecies];
for (ispecies = 0; ispecies < nspecies; ispecies++)
eosSpecies[ispecies] = -1;
-
+
tb->ninput = 0;
char *word = strtok(line," \t\n\r\f");
@@ -540,9 +542,9 @@ void FixEOStableRX::param_extract(Table *tb, char *line)
while (word) {
for (ispecies = 0; ispecies < nspecies; ispecies++)
if (strcmp(word,&atom->dname[ispecies][0]) == 0){
- eosSpecies[ncolumn] = ispecies;
- ncolumn++;
- break;
+ eosSpecies[ncolumn] = ispecies;
+ ncolumn++;
+ break;
}
if (ispecies == nspecies){
printf("name=%s not found in species list\n",word);
@@ -552,7 +554,7 @@ void FixEOStableRX::param_extract(Table *tb, char *line)
}
for (int icolumn = 0; icolumn < ncolumn; icolumn++)
- if(eosSpecies[icolumn]==-1)
+ if(eosSpecies[icolumn]==-1)
error->one(FLERR,"EOS data is missing from fix eos/table/rx tabe");
if(ncolumn != nspecies){
printf("ncolumns=%d nspecies=%d\n",ncolumn,nspecies);
@@ -648,8 +650,8 @@ void FixEOStableRX::energy_lookup(int id, double thetai, double &ui)
int itable;
double fraction, uTmp, nTotal;
- ui = double(0.0);
- nTotal = double(0.0);
+ ui = 0.0;
+ nTotal = 0.0;
for(int ispecies=0;ispecies<nspecies;ispecies++){
Table *tb = &tables[ispecies];
thetai = MAX(thetai,tb->lo);
@@ -659,7 +661,7 @@ void FixEOStableRX::energy_lookup(int id, double thetai, double &ui)
itable = static_cast<int> ((thetai - tb->lo) * tb->invdelta);
fraction = (thetai - tb->r[itable]) * tb->invdelta;
uTmp = tb->e[itable] + fraction*tb->de[itable];
-
+
uTmp += dHf[ispecies];
// mol fraction form:
ui += atom->dvector[ispecies][id]*uTmp;
@@ -691,11 +693,11 @@ void FixEOStableRX::temperature_lookup(int id, double ui, double &thetai)
f1 = u1 - ui;
// Compute guess of t2
- t2 = (double(1.0) + double(0.001))*t1;
+ t2 = (1.0 + 0.001)*t1;
// Compute u2 at t2
energy_lookup(id,t2,u2);
-
+
// Compute f1
f2 = u2 - ui;
diff --git a/src/USER-DPD/fix_eos_table_rx.h b/src/USER-DPD/fix_eos_table_rx.h
index bb2e507c44509cea2cc01c8f76ba793d4831a4fa..b513e292befc66ce2f13ed135982444b2eb7146d 100644
--- a/src/USER-DPD/fix_eos_table_rx.h
+++ b/src/USER-DPD/fix_eos_table_rx.h
@@ -1,4 +1,4 @@
-/* ----------------------------------------------------------------------
+/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
@@ -66,7 +66,7 @@ class FixEOStableRX : public Fix {
int nspecies;
void read_file(char *);
-
+
double *dHf;
int pack_reverse_comm(int, int, double *);
@@ -116,7 +116,7 @@ Self-explanatory.
E: Incorrect format in eos table/rx file
-Self-explanatory.
+Self-explanatory.
E: Cannot open file %s
diff --git a/src/USER-DPD/fix_rx.cpp b/src/USER-DPD/fix_rx.cpp
index 10b297753ea88e8dcdadc768a45c6ad82cef67ff..ce1cb46183d14a6cb1c2c42c54abf69d64e4a563 100644
--- a/src/USER-DPD/fix_rx.cpp
+++ b/src/USER-DPD/fix_rx.cpp
@@ -50,6 +50,8 @@ FixRX::FixRX(LAMMPS *lmp, int narg, char **arg) :
params = NULL;
mol2param = NULL;
pairDPDE = NULL;
+ id_fix_species = NULL;
+ id_fix_species_old = NULL;
// Keep track of the argument list.
int iarg = 3;
@@ -110,11 +112,19 @@ FixRX::~FixRX()
delete [] stoichReactants[ii];
delete [] stoichProducts[ii];
}
+ delete [] Arr;
+ delete [] nArr;
+ delete [] Ea;
+ delete [] tempExp;
delete [] stoich;
delete [] stoichReactants;
delete [] stoichProducts;
delete [] kR;
-}
+ delete [] id_fix_species;
+ delete [] id_fix_species_old;
+
+
+}
/* ---------------------------------------------------------------------- */
@@ -125,10 +135,10 @@ void FixRX::post_constructor()
bool match;
for (int i = 0; i < modify->nfix; i++)
- if (strncmp(modify->fix[i]->style,"property/atom",13) == 0)
- error->all(FLERR,"fix rx cannot be combined with fix property/atom");
+ if (strncmp(modify->fix[i]->style,"property/atom",13) == 0)
+ error->all(FLERR,"fix rx cannot be combined with fix property/atom");
- char **tmpspecies = new char*[maxspecies];
+ char **tmpspecies = new char*[maxspecies];
for(int jj=0; jj < maxspecies; jj++)
tmpspecies[jj] = NULL;
@@ -137,7 +147,7 @@ void FixRX::post_constructor()
FILE *fp;
fp = NULL;
if (comm->me == 0) {
- fp = fopen(kineticsFile,"r");
+ fp = force->open_potential(kineticsFile);
if (fp == NULL) {
char str[128];
sprintf(str,"Cannot open rx file %s",kineticsFile);
@@ -147,7 +157,7 @@ void FixRX::post_constructor()
// Assign species names to tmpspecies array and determine the number of unique species
- int n,nwords,ispecies;
+ int n,nwords;
char line[MAXLINE],*ptr;
int eof = 0;
char * word;
@@ -179,17 +189,17 @@ void FixRX::post_constructor()
word = strtok(NULL, " \t\n\r\f");
match=false;
for(int jj=0;jj<nUniqueSpecies;jj++){
- if(strcmp(word,tmpspecies[jj])==0){
- match=true;
- break;
- }
+ if(strcmp(word,tmpspecies[jj])==0){
+ match=true;
+ break;
+ }
}
- if(!match){
- if(nUniqueSpecies+1>=maxspecies)
- error->all(FLERR,"Exceeded the maximum number of species permitted in fix rx.");
- tmpspecies[nUniqueSpecies] = new char[strlen(word)];
- strcpy(tmpspecies[nUniqueSpecies],word);
- nUniqueSpecies++;
+ if(!match){
+ if(nUniqueSpecies+1>=maxspecies)
+ error->all(FLERR,"Exceeded the maximum number of species permitted in fix rx.");
+ tmpspecies[nUniqueSpecies] = new char[strlen(word)+1];
+ strcpy(tmpspecies[nUniqueSpecies],word);
+ nUniqueSpecies++;
}
word = strtok(NULL, " \t\n\r\f");
if(strcmp(word,"+") != 0 && strcmp(word,"=") != 0) break;
@@ -202,9 +212,6 @@ void FixRX::post_constructor()
// new id = fix-ID + FIX_STORE_ATTRIBUTE
// new fix group = group for this fix
- id_fix_species = NULL;
- id_fix_species_old = NULL;
-
n = strlen(id) + strlen("_SPECIES") + 1;
id_fix_species = new char[n];
n = strlen(id) + strlen("_SPECIES_OLD") + 1;
@@ -231,8 +238,8 @@ void FixRX::post_constructor()
strncat(str1,tmpspecies[ii],strlen(tmpspecies[ii]));
strncat(str2,tmpspecies[ii],strlen(tmpspecies[ii]));
strncat(str2,"Old",3);
- newarg[ii+3] = new char[strlen(str1)];
- newarg2[ii+3] = new char[strlen(str2)];
+ newarg[ii+3] = new char[strlen(str1)+1];
+ newarg2[ii+3] = new char[strlen(str2)+1];
strcpy(newarg[ii+3],str1);
strcpy(newarg2[ii+3],str2);
}
@@ -249,11 +256,15 @@ void FixRX::post_constructor()
if(nspecies==0) error->all(FLERR,"There are no rx species specified.");
- for(int jj=0;jj<maxspecies;jj++) delete tmpspecies[jj];
- delete tmpspecies;
+ for(int jj=0;jj<nspecies;jj++) {
+ delete[] tmpspecies[jj];
+ delete[] newarg[jj+3];
+ delete[] newarg2[jj+3];
+ }
- delete newarg;
- delete newarg2;
+ delete[] newarg;
+ delete[] newarg2;
+ delete[] tmpspecies;
read_file( kineticsFile );
@@ -282,7 +293,7 @@ void FixRX::init()
bool eos_flag = false;
for (int i = 0; i < modify->nfix; i++)
if (strcmp(modify->fix[i]->style,"eos/table/rx") == 0) eos_flag = true;
- if(!eos_flag) error->all(FLERR,"fix rx requires fix eos/table/rx to be specified");
+ if(!eos_flag) error->all(FLERR,"fix rx requires fix eos/table/rx to be specified");
}
/* ---------------------------------------------------------------------- */
@@ -292,20 +303,19 @@ void FixRX::setup_pre_force(int vflag)
int nlocal = atom->nlocal;
int nghost = atom->nghost;
int *mask = atom->mask;
- double *dpdTheta = atom->dpdTheta;
int newton_pair = force->newton_pair;
- double tmp, theta;
+ double tmp;
int ii;
if(localTempFlag){
if (newton_pair) {
dpdThetaLocal = new double[nlocal+nghost];
for (ii = 0; ii < nlocal+nghost; ii++)
- dpdThetaLocal[ii] = double(0.0);
+ dpdThetaLocal[ii] = 0.0;
} else {
dpdThetaLocal = new double[nlocal];
for (ii = 0; ii < nlocal; ii++)
- dpdThetaLocal[ii] = double(0.0);
+ dpdThetaLocal[ii] = 0.0;
}
computeLocalTemperature();
}
@@ -317,12 +327,10 @@ void FixRX::setup_pre_force(int vflag)
}
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit){
- if(localTempFlag) theta = dpdThetaLocal[i];
- else theta=dpdTheta[i];
// Set the reaction rate constants to zero: no reactions occur at step 0
for(int irxn=0;irxn<nreactions;irxn++)
- kR[irxn] = double(0.0);
+ kR[irxn] = 0.0;
if(odeIntegrationFlag==ODE_LAMMPS_RK4) rk4(i);
}
@@ -347,14 +355,13 @@ void FixRX::pre_force(int vflag)
if (newton_pair) {
dpdThetaLocal = new double[nlocal+nghost];
for (ii = 0; ii < nlocal+nghost; ii++)
- dpdThetaLocal[ii] = double(0.0);
+ dpdThetaLocal[ii] = 0.0;
} else {
dpdThetaLocal = new double[nlocal];
for (ii = 0; ii < nlocal; ii++)
- dpdThetaLocal[ii] = double(0.0);
+ dpdThetaLocal[ii] = 0.0;
}
computeLocalTemperature();
- double *dpdTheta = this->dpdThetaLocal;
}
for (int i = 0; i < nlocal; i++)
@@ -364,7 +371,7 @@ void FixRX::pre_force(int vflag)
//Compute the reaction rate constants
for(int irxn=0;irxn<nreactions;irxn++)
- kR[irxn] = Arr[irxn]*pow(theta,nArr[irxn])*exp(-Ea[irxn]/force->boltz/theta);
+ kR[irxn] = Arr[irxn]*pow(theta,nArr[irxn])*exp(-Ea[irxn]/force->boltz/theta);
if(odeIntegrationFlag==ODE_LAMMPS_RK4) rk4(i);
}
@@ -384,7 +391,7 @@ void FixRX::read_file(char *file)
FILE *fp;
fp = NULL;
if (comm->me == 0) {
- fp = fopen(file,"r");
+ fp = force->open_potential(file);
if (fp == NULL) {
char str[128];
sprintf(str,"Cannot open rx file %s",file);
@@ -421,7 +428,7 @@ void FixRX::read_file(char *file)
}
// open file on proc 0
- if (comm->me == 0) fp = fopen(file,"r");
+ if (comm->me == 0) fp = force->open_potential(file);
// read each reaction from kinetics file
eof=0;
@@ -444,14 +451,14 @@ void FixRX::read_file(char *file)
kR = new double[nreactions];
for (int ii=0;ii<nreactions;ii++){
for (int jj=0;jj<nspecies;jj++){
- stoich[ii][jj] = double(0.0);
- stoichReactants[ii][jj] = double(0.0);
- stoichProducts[ii][jj] = double(0.0);
+ stoich[ii][jj] = 0.0;
+ stoichReactants[ii][jj] = 0.0;
+ stoichProducts[ii][jj] = 0.0;
}
}
nreactions=0;
- sign = double(-1.0);
+ sign = -1.0;
while (1) {
if (comm->me == 0) {
ptr = fgets(line,MAXLINE,fp);
@@ -479,31 +486,37 @@ void FixRX::read_file(char *file)
tmpStoich = atof(word);
word = strtok(NULL, " \t\n\r\f");
for (ispecies = 0; ispecies < nspecies; ispecies++){
- if (strcmp(word,&atom->dname[ispecies][0]) == 0){
- stoich[nreactions][ispecies] += sign*tmpStoich;
- if(sign<double(0.0)) stoichReactants[nreactions][ispecies] += tmpStoich;
- else stoichProducts[nreactions][ispecies] += tmpStoich;
- break;
- }
+ if (strcmp(word,&atom->dname[ispecies][0]) == 0){
+ stoich[nreactions][ispecies] += sign*tmpStoich;
+ if(sign<0.0)
+ stoichReactants[nreactions][ispecies] += tmpStoich;
+ else stoichProducts[nreactions][ispecies] += tmpStoich;
+ break;
+ }
}
if(ispecies==nspecies){
- printf("%s mol fraction is not found in data file\n",word);
- printf("nspecies=%d ispecies=%d\n",nspecies,ispecies);
- error->all(FLERR,"Illegal fix rx command");
+ if (comm->me) {
+ fprintf(stderr,"%s mol fraction is not found in data file\n",word);
+ fprintf(stderr,"nspecies=%d ispecies=%d\n",nspecies,ispecies);
+ }
+ error->all(FLERR,"Illegal fix rx command");
}
word = strtok(NULL, " \t\n\r\f");
- if(strcmp(word,"=") == 0) sign = double(1.0);
- if(strcmp(word,"+") != 0 && strcmp(word,"=") != 0){
- if(word==NULL) error->all(FLERR,"Missing parameters in reaction kinetic equation");
- Arr[nreactions] = atof(word);
- word = strtok(NULL, " \t\n\r\f");
- if(word==NULL) error->all(FLERR,"Missing parameters in reaction kinetic equation");
- nArr[nreactions] = atof(word);
- word = strtok(NULL, " \t\n\r\f");
- if(word==NULL) error->all(FLERR,"Missing parameters in reaction kinetic equation");
- Ea[nreactions] = atof(word);
- sign = double(-1.0);
- break;
+ if(word==NULL)
+ error->all(FLERR,"Missing parameters in reaction kinetic equation");
+ if(strcmp(word,"=") == 0) sign = 1.0;
+ if(strcmp(word,"+") != 0 && strcmp(word,"=") != 0) {
+ Arr[nreactions] = atof(word);
+ word = strtok(NULL, " \t\n\r\f");
+ if(word==NULL)
+ error->all(FLERR,"Missing parameters in reaction kinetic equation");
+ nArr[nreactions] = atof(word);
+ word = strtok(NULL, " \t\n\r\f");
+ if(word==NULL)
+ error->all(FLERR,"Missing parameters in reaction kinetic equation");
+ Ea[nreactions] = atof(word);
+ sign = -1.0;
+ break;
}
word = strtok(NULL, " \t\n\r\f");
}
@@ -527,8 +540,8 @@ void FixRX::setupParams()
n = -1;
for (j = 0; j < nreactions; j++) {
if (i == params[j].ispecies) {
- if (n >= 0) error->all(FLERR,"Potential file has duplicate entry");
- n = j;
+ if (n >= 0) error->all(FLERR,"Potential file has duplicate entry");
+ n = j;
}
}
mol2param[i] = n;
@@ -568,13 +581,13 @@ void FixRX::rk4(int id)
// k2
for (int ispecies = 0; ispecies < nspecies; ispecies++)
- yp[ispecies] = y[ispecies] + double(0.5)*h*k1[ispecies];
+ yp[ispecies] = y[ispecies] + 0.5*h*k1[ispecies];
rhs(0.0,yp,k2,dummyArray);
// k3
for (int ispecies = 0; ispecies < nspecies; ispecies++)
- yp[ispecies] = y[ispecies] + double(0.5)*h*k2[ispecies];
+ yp[ispecies] = y[ispecies] + 0.5*h*k2[ispecies];
rhs(0.0,yp,k3,dummyArray);
@@ -585,17 +598,17 @@ void FixRX::rk4(int id)
rhs(0.0,yp,k4,dummyArray);
for (int ispecies = 0; ispecies < nspecies; ispecies++)
- y[ispecies] += h*(k1[ispecies]/6.0 + k2[ispecies]/3.0 +
- k3[ispecies]/3.0 + k4[ispecies]/6.0);
+ y[ispecies] += h*(k1[ispecies]/6.0 + k2[ispecies]/3.0 +
+ k3[ispecies]/3.0 + k4[ispecies]/6.0);
} // end for (int step...
// Store the solution back in atom->dvector.
for (int ispecies = 0; ispecies < nspecies; ispecies++){
- if(y[ispecies] < double(-1.0e-10))
+ if(y[ispecies] < -1.0e-10)
error->one(FLERR,"Computed concentration in RK4 solver is < -1.0e-10");
- else if(y[ispecies] < double(0.0))
- y[ispecies] = double(0.0);
+ else if(y[ispecies] < 0.0)
+ y[ispecies] = 0.0;
atom->dvector[ispecies][id] = y[ispecies];
}
delete [] k1;
@@ -612,7 +625,7 @@ int FixRX::rhs(double t, const double *y, double *dydt, void *params)
int nspecies = atom->nspecies_dpd;
for(int ispecies=0; ispecies<nspecies; ispecies++)
- dydt[ispecies] = double(0.0);
+ dydt[ispecies] = 0.0;
// Construct the reaction rate laws
for(int jrxn=0; jrxn<nreactions; jrxn++){
@@ -624,7 +637,7 @@ int FixRX::rhs(double t, const double *y, double *dydt, void *params)
}
rxnRateLaw[jrxn] = rxnRateLawForward;
}
-
+
// Construct the reaction rates for each species
for(int ispecies=0; ispecies<nspecies; ispecies++)
for(int jrxn=0; jrxn<nreactions; jrxn++)
@@ -649,18 +662,18 @@ void FixRX::computeLocalTemperature()
int newton_pair = force->newton_pair;
// local temperature variables
- double wij, fr, fr4;
+ double wij = 0.0;
double *dpdTheta = atom->dpdTheta;
// Initialize the local density and local temperature arrays
if (newton_pair) {
sumWeights = new double[nlocal+nghost];
for (ii = 0; ii < nlocal+nghost; ii++)
- sumWeights[ii] = double(0.0);
+ sumWeights[ii] = 0.0;
} else {
sumWeights = new double[nlocal];
for (ii = 0; ii < nlocal; ii++)
- dpdThetaLocal[ii] = double(0.0);
+ dpdThetaLocal[ii] = 0.0;
}
inum = pairDPDE->list->inum;
@@ -690,24 +703,22 @@ void FixRX::computeLocalTemperature()
if (rsq < pairDPDE->cutsq[itype][jtype]) {
double rcut = sqrt(pairDPDE->cutsq[itype][jtype]);
- double rij = sqrt(rsq);
- double tmpFactor = double(1.0)-rij/rcut;
- double tmpFactor4 = tmpFactor*tmpFactor*tmpFactor*tmpFactor;
- double ratio = rij/rcut;
-
- // Lucy's Weight Function
- if(wtFlag==LUCY){
- wij = (double(1.0)+double(3.0)*ratio) * (double(1.0)-ratio)*(double(1.0)-ratio)*(double(1.0)-ratio);
- dpdThetaLocal[i] += wij/dpdTheta[j];
- if (newton_pair || j < nlocal)
- dpdThetaLocal[j] += wij/dpdTheta[i];
- }
-
- sumWeights[i] += wij;
+ double rij = sqrt(rsq);
+ double ratio = rij/rcut;
+
+ // Lucy's Weight Function
+ if(wtFlag==LUCY){
+ wij = (1.0+3.0*ratio) * (1.0-ratio)*(1.0-ratio)*(1.0-ratio);
+ dpdThetaLocal[i] += wij/dpdTheta[j];
+ if (newton_pair || j < nlocal)
+ dpdThetaLocal[j] += wij/dpdTheta[i];
+ }
+
+ sumWeights[i] += wij;
if (newton_pair || j < nlocal) {
- sumWeights[j] += wij;
- }
-
+ sumWeights[j] += wij;
+ }
+
}
}
}
@@ -718,17 +729,17 @@ void FixRX::computeLocalTemperature()
// Lucy Weight Function
if(wtFlag==LUCY){
- wij = double(1.0);
+ wij = 1.0;
dpdThetaLocal[i] += wij / dpdTheta[i];
}
sumWeights[i] += wij;
-
+
// Normalized local temperature
dpdThetaLocal[i] = dpdThetaLocal[i] / sumWeights[i];
-
+
if(localTempFlag == HARMONIC)
- dpdThetaLocal[i] = double(1.0) / dpdThetaLocal[i];
-
+ dpdThetaLocal[i] = 1.0 / dpdThetaLocal[i];
+
}
delete [] sumWeights;
diff --git a/src/USER-DPD/fix_rx.h b/src/USER-DPD/fix_rx.h
index 912b75ed5fdefe0719a58020069f59bd8462d3ee..97b9bd55cc5355e8359e92dfe8b7586c0608f009 100644
--- a/src/USER-DPD/fix_rx.h
+++ b/src/USER-DPD/fix_rx.h
@@ -1,4 +1,4 @@
-/* ----------------------------------------------------------------------
+/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
@@ -129,6 +129,6 @@ Self-explanatory
E: Computed concentration in RK4 solver is < -1.0e-10.
-Self-explanatory: Adjust settings for the RK4 solver.
+Self-explanatory: Adjust settings for the RK4 solver.
*/
diff --git a/src/USER-DPD/fix_shardlow.cpp b/src/USER-DPD/fix_shardlow.cpp
index 919689e6189b7892f0a58282ffdfb16188873a54..7248df2e67b40215e2e57914c894593c4925b0d2 100644
--- a/src/USER-DPD/fix_shardlow.cpp
+++ b/src/USER-DPD/fix_shardlow.cpp
@@ -12,11 +12,11 @@
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
- Contributing authors:
+ Contributing authors:
James Larentzos (U.S. Army Research Laboratory)
and Timothy I. Mattox (Engility Corporation)
- Martin Lisal (Institute of Chemical Process Fundamentals
+ 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)
@@ -24,7 +24,7 @@
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",
+ 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
@@ -96,15 +96,17 @@ FixShardlow::FixShardlow(LAMMPS *lmp, int narg, char **arg) :
pairDPDE = (PairDPDfdtEnergy *) force->pair_match("dpd/fdt/energy",1);
if(pairDPDE){
- comm_forward = 10;
+ comm_forward = 3;
comm_reverse = 5;
} else {
- comm_forward = 6;
+ comm_forward = 3;
comm_reverse = 3;
}
if(pairDPD == NULL && pairDPDE == NULL)
error->all(FLERR,"Must use pair_style dpd/fdt or dpd/fdt/energy with fix shardlow");
+ if(!(atom->dpd_flag))
+ error->all(FLERR,"Must use atom_style dpd with fix shardlow");
}
/* ---------------------------------------------------------------------- */
@@ -113,7 +115,6 @@ int FixShardlow::setmask()
{
int mask = 0;
mask |= INITIAL_INTEGRATE;
- mask |= PRE_NEIGHBOR;
return mask;
}
@@ -144,13 +145,6 @@ void FixShardlow::setup(int vflag)
}
}
-/* ---------------------------------------------------------------------- */
-
-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
@@ -176,14 +170,12 @@ void FixShardlow::initial_integrate(int vflag)
int *type = atom->type;
int nlocal = atom->nlocal;
int nghost = atom->nghost;
- int nall = nlocal + nghost;
int newton_pair = force->newton_pair;
double randPair;
+ int *ssaAIR = atom->ssaAIR;
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;
double vxi, vyi, vzi, vxj, vyj, vzj;
@@ -203,7 +195,7 @@ void FixShardlow::initial_integrate(int vflag)
double bby = domain->subhi[1] - domain->sublo[1];
double bbz = domain->subhi[2] - domain->sublo[2];
- double rcut = double(2.0)*neighbor->cutneighmax;
+ double rcut = 2.0*neighbor->cutneighmax;
if (domain->triclinic)
error->all(FLERR,"Fix shardlow does not yet support triclinic geometries");
@@ -211,25 +203,8 @@ void FixShardlow::initial_integrate(int vflag)
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);
+ // Allocate memory for v_t0 to hold the initial velocities for the ghosts
+ v_t0 = (double (*)[3]) memory->smalloc(sizeof(double)*3*nghost, "FixShardlow:v_t0");
// Define pointers to access the neighbor list
if(pairDPDE){
@@ -244,9 +219,20 @@ void FixShardlow::initial_integrate(int vflag)
firstneigh = pairDPD->list->firstneigh;
}
- //Loop over all 14 directions (8 stages)
+ //Loop over all 14 directions (8 stages)
for (int idir = 1; idir <=8; idir++){
-
+
+ if (idir > 1) {
+ // Communicate the updated velocities to all nodes
+ comm->forward_comm_fix(this);
+
+ if(pairDPDE){
+ // Zero out the ghosts' uCond & uMech to be used as delta accumulators
+ memset(&(uCond[nlocal]), 0, sizeof(double)*nghost);
+ memset(&(uMech[nlocal]), 0, sizeof(double)*nghost);
+ }
+ }
+
// Loop over neighbors of my atoms
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
@@ -254,221 +240,190 @@ void FixShardlow::initial_integrate(int vflag)
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
-
+
+ // load velocity for i from memory
+ vxi = v[i][0];
+ vyi = v[i][1];
+ vzi = v[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);
+ j = jlist[jj];
+ j &= NEIGHMASK;
+ if (ssaAIR[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 = 1.0/r;
+
+ // Keep a copy of the velocities from previous Shardlow step
+ vx0i = vxi;
+ vy0i = vyi;
+ vz0i = vzi;
+
+ vx0j = vxj = v[j][0];
+ vy0j = vyj = v[j][1];
+ vz0j = vzj = v[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 = 1.0 - r/pairDPD->cut[itype][jtype];
+ wr = 1.0 - r/cut;
+ wd = wr*wr;
+
+ if(pairDPDE){
+ // Compute the current temperature
+ theta_ij = 0.5*(1.0/dpdTheta[i] + 1.0/dpdTheta[j]);
+ theta_ij = 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 *= 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 = 1.0 / mass_i;
+ massinv_j = 1.0 / mass_j;
+
+ // Update the velocity on i
+ vxi += dpx*force->ftm2v*massinv_i;
+ vyi += dpy*force->ftm2v*massinv_i;
+ vzi += dpz*force->ftm2v*massinv_i;
+
+ if (newton_pair || j < nlocal) {
+ // Update the velocity on j
+ vxj -= dpx*force->ftm2v*massinv_j;
+ vyj -= dpy*force->ftm2v*massinv_j;
+ vzj -= dpz*force->ftm2v*massinv_j;
+ }
+
+ //ii. Compute the velocity diff
+ delvx = vxi - vxj;
+ delvy = vyi - vyj;
+ delvz = vzi - vzj;
+
+ dot = delx*delvx + dely*delvy + delz*delvz;
+
+ //iii. Compute dpi again
+ mu_ij = massinv_i + massinv_j;
+ denom = 1.0 + 0.5*mu_ij*gamma_ij*wd*dt*force->ftm2v;
+ factor_dpd = -0.5*dt*gamma_ij*wd*force->ftm2v/denom;
+ factor_dpd1 = factor_dpd*(mu_ij*randPair);
+ factor_dpd1 += randPair;
+ factor_dpd1 *= 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 on i
+ vxi += dpx*force->ftm2v*massinv_i;
+ vyi += dpy*force->ftm2v*massinv_i;
+ vzi += dpz*force->ftm2v*massinv_i;
+
+ if (newton_pair || j < nlocal) {
+ // Update the velocity on j
+ vxj -= dpx*force->ftm2v*massinv_j;
+ vyj -= dpy*force->ftm2v*massinv_j;
+ vzj -= dpz*force->ftm2v*massinv_j;
+ // Store updated velocity for j
+ v[j][0] = vxj;
+ v[j][1] = vyj;
+ v[j][2] = vzj;
+ }
+
+ 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*(1.0/dpdTheta[i] - 1.0/dpdTheta[j])*wd*dt;
+ factor_dpd += randPair;
+
+ uCond[i] += factor_dpd;
+ if (newton_pair || j < nlocal) {
+ uCond[j] -= factor_dpd;
+ }
+
+ // Compute uMech
+ 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 = 0.25*(dot1+dot2-dot3-dot4)/force->ftm2v;
+ uMech[i] -= factor_dpd;
+ if (newton_pair || j < nlocal) {
+ uMech[j] -= factor_dpd;
+ }
+ }
+ }
}
+ // store updated velocity for i
+ v[i][0] = vxi;
+ v[i][1] = vyi;
+ v[i][2] = vzi;
}
-
- // 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
-------------------------------------------------------------------------- */
+ // Communicate the ghost deltas to the atom owners
+ if (idir > 1) comm->reverse_comm_fix(this);
-void FixShardlow::setup_pre_neighbor()
-{
- neighbor->assign_ssa_airnums();
-}
+ } //End Loop over all directions For idir = Top, Top-Right, Right, Bottom-Right, Back
-void FixShardlow::pre_neighbor()
-{
- neighbor->assign_ssa_airnums();
+ memory->sfree(v_t0);
+ v_t0 = NULL;
}
/* ---------------------------------------------------------------------- */
@@ -477,26 +432,13 @@ int FixShardlow::pack_forward_comm(int n, int *list, double *buf, int pbc_flag,
{
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;
}
@@ -506,27 +448,15 @@ int FixShardlow::pack_forward_comm(int n, int *list, double *buf, int pbc_flag,
void FixShardlow::unpack_forward_comm(int n, int first, double *buf)
{
int ii,m,last;
+ int nlocal = atom->nlocal;
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++];
- }
+ v_t0[ii - nlocal][0] = v[ii][0] = buf[m++];
+ v_t0[ii - nlocal][1] = v[ii][1] = buf[m++];
+ v_t0[ii - nlocal][2] = v[ii][2] = buf[m++];
}
}
@@ -535,18 +465,20 @@ void FixShardlow::unpack_forward_comm(int n, int first, double *buf)
int FixShardlow::pack_reverse_comm(int n, int first, double *buf)
{
int i,m,last;
- double *duCond = atom->duCond;
- double *duMech = atom->duMech;
+ int nlocal = atom->nlocal;
+ double **v = atom->v;
+ double *uCond = atom->uCond;
+ double *uMech = atom->uMech;
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];
+ buf[m++] = v[i][0] - v_t0[i - nlocal][0];
+ buf[m++] = v[i][1] - v_t0[i - nlocal][1];
+ buf[m++] = v[i][2] - v_t0[i - nlocal][2];
if(pairDPDE){
- buf[m++] = duCond[i];
- buf[m++] = duMech[i];
+ buf[m++] = uCond[i]; // for ghosts, this is an accumulated delta
+ buf[m++] = uMech[i]; // for ghosts, this is an accumulated delta
}
}
return m;
@@ -557,19 +489,20 @@ int FixShardlow::pack_reverse_comm(int n, int first, double *buf)
void FixShardlow::unpack_reverse_comm(int n, int *list, double *buf)
{
int i,j,m;
- double *duCond = atom->duCond;
- double *duMech = atom->duMech;
+ double **v = atom->v;
+ double *uCond = atom->uCond;
+ double *uMech = atom->uMech;
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++];
+ v[j][0] += buf[m++];
+ v[j][1] += buf[m++];
+ v[j][2] += buf[m++];
if(pairDPDE){
- duCond[j] += buf[m++];
- duMech[j] += buf[m++];
+ uCond[j] += buf[m++]; // add in the accumulated delta
+ uMech[j] += buf[m++]; // add in the accumulated delta
}
}
}
diff --git a/src/USER-DPD/fix_shardlow.h b/src/USER-DPD/fix_shardlow.h
index 327e6357dc813c6c93a1af148d4a1056a485dd04..6421737b3d077a3ec160caaab118b0d802ab425f 100644
--- a/src/USER-DPD/fix_shardlow.h
+++ b/src/USER-DPD/fix_shardlow.h
@@ -1,4 +1,4 @@
-/* ----------------------------------------------------------------------
+/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
@@ -31,12 +31,8 @@ class FixShardlow : public Fix {
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 *);
@@ -45,7 +41,7 @@ class FixShardlow : public Fix {
class PairDPDfdt *pairDPD;
class PairDPDfdtEnergy *pairDPDE;
- double **dvSSA;
+ double (*v_t0)[3];
private:
class NeighList *list;
@@ -71,7 +67,7 @@ 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
+E: A deterministic integrator must be specified after fix shardlow in input
file (e.g. fix nve or fix nph).
Self-explanatory.
@@ -84,11 +80,11 @@ E: Fix shardlow does not yet support triclinic geometries
Self-explanatory.
-E: Shardlow algorithm requires sub-domain length > 2*(rcut+skin). Either
+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
+The Shardlow splitting algorithm requires the size of the sub-domain lengths
+to be are larger than twice the cutoff+skin. Generally, the domain decomposition
is dependant on the number of processors requested.
*/
diff --git a/src/USER-DPD/pair_dpd_fdt.cpp b/src/USER-DPD/pair_dpd_fdt.cpp
index de313b820e9fa8d4e567a11d004de423ed540626..21c5e420fa179bee600476f984118e5828195597 100644
--- a/src/USER-DPD/pair_dpd_fdt.cpp
+++ b/src/USER-DPD/pair_dpd_fdt.cpp
@@ -114,12 +114,12 @@ void PairDPDfdt::compute(int eflag, int vflag)
if (r < EPSILON) continue; // r can be 0.0 in DPD systems
rinv = 1.0/r;
wr = 1.0 - r/cut[itype][jtype];
- wd = wr*wr;
+ wd = wr*wr;
// 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;
@@ -179,7 +179,9 @@ void PairDPDfdt::settings(int narg, char **arg)
temperature = force->numeric(FLERR,arg[0]);
cut_global = force->numeric(FLERR,arg[1]);
seed = force->inumeric(FLERR,arg[2]);
-
+ if (atom->dpd_flag != 1)
+ error->all(FLERR,"pair_style dpd/fdt requires atom_style dpd");
+
// initialize Marsaglia RNG with processor-unique seed
if (seed <= 0) error->all(FLERR,"Illegal pair_style command");
@@ -212,7 +214,7 @@ void PairDPDfdt::coeff(int narg, char **arg)
double a0_one = force->numeric(FLERR,arg[2]);
double sigma_one = force->numeric(FLERR,arg[3]);
double cut_one = cut_global;
-
+
if (narg == 5) cut_one = force->numeric(FLERR,arg[4]);
int count = 0;
diff --git a/src/USER-DPD/pair_dpd_fdt.h b/src/USER-DPD/pair_dpd_fdt.h
index a11d963ef9fdf2787d4ff248cdf9f317559df66f..f7721853befaf61d810c7b09c8302f1c4f46d310 100644
--- a/src/USER-DPD/pair_dpd_fdt.h
+++ b/src/USER-DPD/pair_dpd_fdt.h
@@ -1,4 +1,4 @@
-/* ----------------------------------------------------------------------
+/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
diff --git a/src/USER-DPD/pair_dpd_fdt_energy.cpp b/src/USER-DPD/pair_dpd_fdt_energy.cpp
index 2fd3c0bf55864bd213d2fc5b1c81a484d080b729..34e21918d1d36ad6ac32427367e6f1c5af0818ec 100644
--- a/src/USER-DPD/pair_dpd_fdt_energy.cpp
+++ b/src/USER-DPD/pair_dpd_fdt_energy.cpp
@@ -115,12 +115,12 @@ void PairDPDfdtEnergy::compute(int eflag, int vflag)
if (r < EPSILON) continue; // r can be 0.0 in DPD systems
rinv = 1.0/r;
wr = 1.0 - r/cut[itype][jtype];
- wd = wr*wr;
+ wd = wr*wr;
// 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;
@@ -180,9 +180,9 @@ void PairDPDfdtEnergy::settings(int narg, char **arg)
cut_global = force->numeric(FLERR,arg[0]);
seed = force->inumeric(FLERR,arg[1]);
- if (atom->dpd_flag != 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");
@@ -259,7 +259,7 @@ void PairDPDfdtEnergy::init_style()
bool eos_flag = false;
for (int i = 0; i < modify->nfix; i++)
if (strncmp(modify->fix[i]->style,"eos",3) == 0) eos_flag = true;
- if(!eos_flag) error->all(FLERR,"pair_style dpd/fdt/energy requires an EOS to be specified");
+ if(!eos_flag) error->all(FLERR,"pair_style dpd/fdt/energy requires an EOS to be specified");
}
/* ----------------------------------------------------------------------
diff --git a/src/USER-DPD/pair_dpd_fdt_energy.h b/src/USER-DPD/pair_dpd_fdt_energy.h
index ef1afed4e820f899a3b641427877579b76eb54d0..c296e8f1a3a3d14da4a9f2bb066b6f55ff3ad6c6 100644
--- a/src/USER-DPD/pair_dpd_fdt_energy.h
+++ b/src/USER-DPD/pair_dpd_fdt_energy.h
@@ -1,4 +1,4 @@
-/* ----------------------------------------------------------------------
+/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
diff --git a/src/USER-DPD/pair_exp6_rx.cpp b/src/USER-DPD/pair_exp6_rx.cpp
index e6aeb9ac7d7329acd681233f38a5825ec2a1cb1c..d4710e1a8d09bf60fde3b9327b833992b73d77a7 100644
--- a/src/USER-DPD/pair_exp6_rx.cpp
+++ b/src/USER-DPD/pair_exp6_rx.cpp
@@ -43,12 +43,17 @@ PairExp6rx::PairExp6rx(LAMMPS *lmp) : Pair(lmp)
nparams = maxparam = 0;
params = NULL;
mol2param = NULL;
+ site1 = site2 = NULL;
}
/* ---------------------------------------------------------------------- */
PairExp6rx::~PairExp6rx()
{
+ for (int i=0; i < nparams; ++i) {
+ delete[] params[i].name;
+ delete[] params[i].potential;
+ }
memory->destroy(params);
memory->destroy(mol2param);
@@ -58,6 +63,8 @@ PairExp6rx::~PairExp6rx()
memory->destroy(cut);
}
+ delete[] site1;
+ delete[] site2;
}
/* ---------------------------------------------------------------------- */
@@ -66,7 +73,7 @@ void PairExp6rx::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum,itype,jtype;
double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,evdwlOld,fpair;
- double rsq,rinv,r2inv,r6inv,forceExp6,factor_lj;
+ double rsq,r2inv,r6inv,forceExp6,factor_lj;
double rCut,rCutInv,rCut2inv,rCut6inv,rCutExp,urc,durc;
double rm2ij,rm6ij;
double r,rexp;
@@ -106,8 +113,8 @@ void PairExp6rx::compute(int eflag, int vflag)
double *uCG = atom->uCG;
double *uCGnew = atom->uCGnew;
- const double nRep = double(12.0);
- const double shift = double(1.05);
+ const double nRep = 12.0;
+ const double shift = 1.05;
double rin1, aRep, uin1, win1, uin1rep, rin1exp, rin6, rin6inv;
inum = list->inum;
@@ -125,7 +132,7 @@ void PairExp6rx::compute(int eflag, int vflag)
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
-
+
getParamsEXP6(i,epsilon1_i,alpha1_i,rm1_i,fraction1_i,epsilon2_i,alpha2_i,rm2_i,fraction2_i,epsilonOld1_i,alphaOld1_i,rmOld1_i,fractionOld1_i,epsilonOld2_i,alphaOld2_i,rmOld2_i,fractionOld2_i);
for (jj = 0; jj < jnum; jj++) {
@@ -141,236 +148,233 @@ void PairExp6rx::compute(int eflag, int vflag)
jtype = type[j];
if (rsq < cutsq[itype][jtype]) {
- r2inv = double(1.0)/rsq;
+ r2inv = 1.0/rsq;
r6inv = r2inv*r2inv*r2inv;
r = sqrt(rsq);
- rinv = double(1.0)/r;
- rCut2inv = double(1.0)/cutsq[itype][jtype];
- rCut6inv = rCut2inv*rCut2inv*rCut2inv;
- rCut = sqrt(cutsq[itype][jtype]);
- rCutInv = double(1.0)/rCut;
-
- //
- // A. Compute the exp-6 potential
- //
-
- // A1. Get alpha, epsilon and rm for particle j
- getParamsEXP6(j,epsilon1_j,alpha1_j,rm1_j,fraction1_j,epsilon2_j,alpha2_j,rm2_j,fraction2_j,epsilonOld1_j,alphaOld1_j,rmOld1_j,fractionOld1_j,epsilonOld2_j,alphaOld2_j,rmOld2_j,fractionOld2_j);
-
- // A2. Apply Lorentz-Berthelot mixing rules for the i-j pair
- alphaOld12_ij = sqrt(alphaOld1_i*alphaOld2_j);
- rmOld12_ij = 0.5*(rmOld1_i + rmOld2_j);
- epsilonOld12_ij = sqrt(epsilonOld1_i*epsilonOld2_j);
- alphaOld21_ij = sqrt(alphaOld2_i*alphaOld1_j);
- rmOld21_ij = 0.5*(rmOld2_i + rmOld1_j);
- epsilonOld21_ij = sqrt(epsilonOld2_i*epsilonOld1_j);
-
- alpha12_ij = sqrt(alpha1_i*alpha2_j);
- rm12_ij = 0.5*(rm1_i + rm2_j);
- epsilon12_ij = sqrt(epsilon1_i*epsilon2_j);
- alpha21_ij = sqrt(alpha2_i*alpha1_j);
- rm21_ij = 0.5*(rm2_i + rm1_j);
- epsilon21_ij = sqrt(epsilon2_i*epsilon1_j);
-
- if(rmOld12_ij!=double(0.0) && rmOld21_ij!=double(0.0)){
- if(alphaOld21_ij == double(6.0) || alphaOld12_ij == double(6.0))
- error->all(FLERR,"alpha_ij is 6.0 in pair exp6");
-
- // A3. Compute some convenient quantities for evaluating the force
- rminv = 1.0/rmOld12_ij;
- buck1 = epsilonOld12_ij / (alphaOld12_ij - 6.0);
- rexp = expValue(alphaOld12_ij*(1.0-r*rminv));
- rm2ij = rmOld12_ij*rmOld12_ij;
- rm6ij = rm2ij*rm2ij*rm2ij;
-
- // Compute the shifted potential
- rCutExp = expValue(alphaOld12_ij*(1.0-rCut*rminv));
- buck2 = 6.0*alphaOld12_ij;
- urc = buck1*(6.0*rCutExp - alphaOld12_ij*rm6ij*rCut6inv);
- durc = -buck1*buck2*(rCutExp* rminv - rCutInv*rm6ij*rCut6inv);
- rin1 = shift*rmOld12_ij*func_rin(alphaOld12_ij);
- if(r < rin1){
- rin6 = rin1*rin1*rin1*rin1*rin1*rin1;
- rin6inv = double(1.0)/rin6;
-
- rin1exp = expValue(alphaOld12_ij*(1.0-rin1*rminv));
-
- uin1 = buck1*(6.0*rin1exp - alphaOld12_ij*rm6ij*rin6inv) - urc - durc*(rin1-rCut);
-
- win1 = buck1*buck2*(rin1*rin1exp*rminv - rm6ij*rin6inv) - rin1*durc;
-
- aRep = double(-1.0)*win1*pow(rin1,nRep)/nRep;
-
- uin1rep = aRep/pow(rin1,nRep);
-
- evdwlOldEXP6_12 = uin1 - uin1rep + aRep/pow(r,nRep);
-
- } else {
- evdwlOldEXP6_12 = buck1*(6.0*rexp - alphaOld12_ij*rm6ij*r6inv) - urc - durc*(r-rCut);
- }
-
- // A3. Compute some convenient quantities for evaluating the force
- rminv = 1.0/rmOld21_ij;
- buck1 = epsilonOld21_ij / (alphaOld21_ij - 6.0);
- buck2 = 6.0*alphaOld21_ij;
- rexp = expValue(alphaOld21_ij*(1.0-r*rminv));
- rm2ij = rmOld21_ij*rmOld21_ij;
- rm6ij = rm2ij*rm2ij*rm2ij;
-
- // Compute the shifted potential
- rCutExp = expValue(alphaOld21_ij*(1.0-rCut*rminv));
- buck2 = 6.0*alphaOld21_ij;
- urc = buck1*(6.0*rCutExp - alphaOld21_ij*rm6ij*rCut6inv);
- durc = -buck1*buck2*(rCutExp* rminv - rCutInv*rm6ij*rCut6inv);
- rin1 = shift*rmOld21_ij*func_rin(alphaOld21_ij);
-
- if(r < rin1){
- rin6 = rin1*rin1*rin1*rin1*rin1*rin1;
- rin6inv = double(1.0)/rin6;
-
- rin1exp = expValue(alphaOld21_ij*(1.0-rin1*rminv));
+ rCut2inv = 1.0/cutsq[itype][jtype];
+ rCut6inv = rCut2inv*rCut2inv*rCut2inv;
+ rCut = sqrt(cutsq[itype][jtype]);
+ rCutInv = 1.0/rCut;
+
+ //
+ // A. Compute the exp-6 potential
+ //
+
+ // A1. Get alpha, epsilon and rm for particle j
+ getParamsEXP6(j,epsilon1_j,alpha1_j,rm1_j,fraction1_j,epsilon2_j,alpha2_j,rm2_j,fraction2_j,epsilonOld1_j,alphaOld1_j,rmOld1_j,fractionOld1_j,epsilonOld2_j,alphaOld2_j,rmOld2_j,fractionOld2_j);
+
+ // A2. Apply Lorentz-Berthelot mixing rules for the i-j pair
+ alphaOld12_ij = sqrt(alphaOld1_i*alphaOld2_j);
+ rmOld12_ij = 0.5*(rmOld1_i + rmOld2_j);
+ epsilonOld12_ij = sqrt(epsilonOld1_i*epsilonOld2_j);
+ alphaOld21_ij = sqrt(alphaOld2_i*alphaOld1_j);
+ rmOld21_ij = 0.5*(rmOld2_i + rmOld1_j);
+ epsilonOld21_ij = sqrt(epsilonOld2_i*epsilonOld1_j);
+
+ alpha12_ij = sqrt(alpha1_i*alpha2_j);
+ rm12_ij = 0.5*(rm1_i + rm2_j);
+ epsilon12_ij = sqrt(epsilon1_i*epsilon2_j);
+ alpha21_ij = sqrt(alpha2_i*alpha1_j);
+ rm21_ij = 0.5*(rm2_i + rm1_j);
+ epsilon21_ij = sqrt(epsilon2_i*epsilon1_j);
+
+ if(rmOld12_ij!=0.0 && rmOld21_ij!=0.0){
+ if(alphaOld21_ij == 6.0 || alphaOld12_ij == 6.0)
+ error->all(FLERR,"alpha_ij is 6.0 in pair exp6");
+
+ // A3. Compute some convenient quantities for evaluating the force
+ rminv = 1.0/rmOld12_ij;
+ buck1 = epsilonOld12_ij / (alphaOld12_ij - 6.0);
+ rexp = expValue(alphaOld12_ij*(1.0-r*rminv));
+ rm2ij = rmOld12_ij*rmOld12_ij;
+ rm6ij = rm2ij*rm2ij*rm2ij;
+
+ // Compute the shifted potential
+ rCutExp = expValue(alphaOld12_ij*(1.0-rCut*rminv));
+ buck2 = 6.0*alphaOld12_ij;
+ urc = buck1*(6.0*rCutExp - alphaOld12_ij*rm6ij*rCut6inv);
+ durc = -buck1*buck2*(rCutExp* rminv - rCutInv*rm6ij*rCut6inv);
+ rin1 = shift*rmOld12_ij*func_rin(alphaOld12_ij);
+ if(r < rin1){
+ rin6 = rin1*rin1*rin1*rin1*rin1*rin1;
+ rin6inv = 1.0/rin6;
+
+ rin1exp = expValue(alphaOld12_ij*(1.0-rin1*rminv));
+
+ uin1 = buck1*(6.0*rin1exp - alphaOld12_ij*rm6ij*rin6inv) - urc - durc*(rin1-rCut);
+
+ win1 = buck1*buck2*(rin1*rin1exp*rminv - rm6ij*rin6inv) - rin1*durc;
+
+ aRep = -1.0*win1*pow(rin1,nRep)/nRep;
+
+ uin1rep = aRep/pow(rin1,nRep);
+
+ evdwlOldEXP6_12 = uin1 - uin1rep + aRep/pow(r,nRep);
+
+ } else {
+ evdwlOldEXP6_12 = buck1*(6.0*rexp - alphaOld12_ij*rm6ij*r6inv) - urc - durc*(r-rCut);
+ }
+
+ // A3. Compute some convenient quantities for evaluating the force
+ rminv = 1.0/rmOld21_ij;
+ buck1 = epsilonOld21_ij / (alphaOld21_ij - 6.0);
+ buck2 = 6.0*alphaOld21_ij;
+ rexp = expValue(alphaOld21_ij*(1.0-r*rminv));
+ rm2ij = rmOld21_ij*rmOld21_ij;
+ rm6ij = rm2ij*rm2ij*rm2ij;
+
+ // Compute the shifted potential
+ rCutExp = expValue(alphaOld21_ij*(1.0-rCut*rminv));
+ buck2 = 6.0*alphaOld21_ij;
+ urc = buck1*(6.0*rCutExp - alphaOld21_ij*rm6ij*rCut6inv);
+ durc = -buck1*buck2*(rCutExp* rminv - rCutInv*rm6ij*rCut6inv);
+ rin1 = shift*rmOld21_ij*func_rin(alphaOld21_ij);
+
+ if(r < rin1){
+ rin6 = rin1*rin1*rin1*rin1*rin1*rin1;
+ rin6inv = 1.0/rin6;
+
+ rin1exp = expValue(alphaOld21_ij*(1.0-rin1*rminv));
+
+ uin1 = buck1*(6.0*rin1exp - alphaOld21_ij*rm6ij*rin6inv) - urc - durc*(rin1-rCut);
+
+ win1 = buck1*buck2*(rin1*rin1exp*rminv - rm6ij*rin6inv) - rin1*durc;
+
+ aRep = -1.0*win1*pow(rin1,nRep)/nRep;
+
+ uin1rep = aRep/pow(rin1,nRep);
+
+ evdwlOldEXP6_21 = uin1 - uin1rep + aRep/pow(r,nRep);
+
+ } else {
+ evdwlOldEXP6_21 = buck1*(6.0*rexp - alphaOld21_ij*rm6ij*r6inv) - urc - durc*(r-rCut);
+ }
+
+ if (strcmp(site1,site2) == 0)
+ evdwlOld = sqrt(fractionOld1_i*fractionOld2_j)*evdwlOldEXP6_12;
+ else
+ evdwlOld = sqrt(fractionOld1_i*fractionOld2_j)*evdwlOldEXP6_12 + sqrt(fractionOld2_i*fractionOld1_j)*evdwlOldEXP6_21;
+
+ evdwlOld *= factor_lj;
+
+ uCG[i] += 0.5*evdwlOld;
+ uCG[j] += 0.5*evdwlOld;
+ }
+
+ if(rm12_ij!=0.0 && rm21_ij!=0.0){
+ if(alpha21_ij == 6.0 || alpha12_ij == 6.0)
+ error->all(FLERR,"alpha_ij is 6.0 in pair exp6");
+
+ // A3. Compute some convenient quantities for evaluating the force
+ rminv = 1.0/rm12_ij;
+ buck1 = epsilon12_ij / (alpha12_ij - 6.0);
+ buck2 = 6.0*alpha12_ij;
+ rexp = expValue(alpha12_ij*(1.0-r*rminv));
+ rm2ij = rm12_ij*rm12_ij;
+ rm6ij = rm2ij*rm2ij*rm2ij;
+
+ // Compute the shifted potential
+ rCutExp = expValue(alpha12_ij*(1.0-rCut*rminv));
+ urc = buck1*(6.0*rCutExp - alpha12_ij*rm6ij*rCut6inv);
+ durc = -buck1*buck2*(rCutExp*rminv - rCutInv*rm6ij*rCut6inv);
+ rin1 = shift*rm12_ij*func_rin(alpha12_ij);
- uin1 = buck1*(6.0*rin1exp - alphaOld21_ij*rm6ij*rin6inv) - urc - durc*(rin1-rCut);
-
- win1 = buck1*buck2*(rin1*rin1exp*rminv - rm6ij*rin6inv) - rin1*durc;
+ if(r < rin1){
+ rin6 = rin1*rin1*rin1*rin1*rin1*rin1;
+ rin6inv = 1.0/rin6;
- aRep = double(-1.0)*win1*pow(rin1,nRep)/nRep;
+ rin1exp = expValue(alpha12_ij*(1.0-rin1*rminv));
- uin1rep = aRep/pow(rin1,nRep);
+ uin1 = buck1*(6.0*rin1exp - alpha12_ij*rm6ij*rin6inv) - urc - durc*(rin1-rCut);
- evdwlOldEXP6_21 = uin1 - uin1rep + aRep/pow(r,nRep);
-
- } else {
- evdwlOldEXP6_21 = buck1*(6.0*rexp - alphaOld21_ij*rm6ij*r6inv) - urc - durc*(r-rCut);
- }
-
- if (strcmp(site1,site2) == 0)
- evdwlOld = sqrt(fractionOld1_i*fractionOld2_j)*evdwlOldEXP6_12;
- else
- evdwlOld = sqrt(fractionOld1_i*fractionOld2_j)*evdwlOldEXP6_12 + sqrt(fractionOld2_i*fractionOld1_j)*evdwlOldEXP6_21;
-
- evdwlOld *= factor_lj;
-
- uCG[i] += double(0.5)*evdwlOld;
- uCG[j] += double(0.5)*evdwlOld;
- }
-
- if(rm12_ij!=double(0.0) && rm21_ij!=double(0.0)){
- if(alpha21_ij == double(6.0) || alpha12_ij == double(6.0))
- error->all(FLERR,"alpha_ij is 6.0 in pair exp6");
-
- // A3. Compute some convenient quantities for evaluating the force
- rminv = 1.0/rm12_ij;
- buck1 = epsilon12_ij / (alpha12_ij - 6.0);
- buck2 = 6.0*alpha12_ij;
- rexp = expValue(alpha12_ij*(1.0-r*rminv));
- rm2ij = rm12_ij*rm12_ij;
- rm6ij = rm2ij*rm2ij*rm2ij;
-
- // Compute the shifted potential
- rCutExp = expValue(alpha12_ij*(1.0-rCut*rminv));
- urc = buck1*(6.0*rCutExp - alpha12_ij*rm6ij*rCut6inv);
- durc = -buck1*buck2*(rCutExp*rminv - rCutInv*rm6ij*rCut6inv);
- rin1 = shift*rm12_ij*func_rin(alpha12_ij);
-
- if(r < rin1){
- rin6 = rin1*rin1*rin1*rin1*rin1*rin1;
- rin6inv = double(1.0)/rin6;
-
- rin1exp = expValue(alpha12_ij*(1.0-rin1*rminv));
-
- uin1 = buck1*(6.0*rin1exp - alpha12_ij*rm6ij*rin6inv) - urc - durc*(rin1-rCut);
-
- win1 = buck1*buck2*(rin1*rin1exp*rminv - rm6ij*rin6inv) - rin1*durc;
-
- aRep = double(-1.0)*win1*pow(rin1,nRep)/nRep;
-
- uin1rep = aRep/pow(rin1,nRep);
-
- evdwlEXP6_12 = uin1 - uin1rep + aRep/pow(r,nRep);
-
- forceExp6 = double(-1.0)*nRep*aRep/pow(r,nRep);
- fpairEXP6_12 = factor_lj*forceExp6*r2inv;
-
- } else {
-
- // A4. Compute the exp-6 force and energy
- forceExp6 = buck1*buck2*(r*rexp*rminv - rm6ij*r6inv) + r*durc;
- fpairEXP6_12 = factor_lj*forceExp6*r2inv;
-
- evdwlEXP6_12 = buck1*(6.0*rexp - alpha12_ij*rm6ij*r6inv) - urc - durc*(r-rCut);
- }
-
- rminv = 1.0/rm21_ij;
- buck1 = epsilon21_ij / (alpha21_ij - 6.0);
- buck2 = 6.0*alpha21_ij;
- rexp = expValue(alpha21_ij*(1.0-r*rminv));
- rm2ij = rm21_ij*rm21_ij;
- rm6ij = rm2ij*rm2ij*rm2ij;
-
- // Compute the shifted potential
- rCutExp = expValue(alpha21_ij*(1.0-rCut*rminv));
- urc = buck1*(6.0*rCutExp - alpha21_ij*rm6ij*rCut6inv);
- durc = -buck1*buck2*(rCutExp*rminv - rCutInv*rm6ij*rCut6inv);
- rin1 = shift*rm21_ij*func_rin(alpha21_ij);
-
- if(r < rin1){
- rin6 = rin1*rin1*rin1*rin1*rin1*rin1;
- rin6inv = double(1.0)/rin6;
-
- rin1exp = expValue(alpha21_ij*(1.0-rin1*rminv));
-
- uin1 = buck1*(6.0*rin1exp - alpha21_ij*rm6ij*rin6inv) - urc - durc*(rin1-rCut);
-
- win1 = buck1*buck2*(rin1*rin1exp*rminv - rm6ij*rin6inv) - rin1*durc;
-
- aRep = double(-1.0)*win1*pow(rin1,nRep)/nRep;
-
- uin1rep = aRep/pow(rin1,nRep);
-
- evdwlEXP6_21 = uin1 - uin1rep + aRep/pow(r,nRep);
-
- forceExp6 = double(-1.0)*nRep*aRep/pow(r,nRep);
- fpairEXP6_21 = factor_lj*forceExp6*r2inv;
-
- } else {
-
- // A4. Compute the exp-6 force and energy
- forceExp6 = buck1*buck2*(r*rexp*rminv - rm6ij*r6inv) + r*durc;
- fpairEXP6_21 = factor_lj*forceExp6*r2inv;
-
- evdwlEXP6_21 = buck1*(6.0*rexp - alpha21_ij*rm6ij*r6inv) - urc - durc*(r-rCut);
- }
-
- //
- // Apply Mixing Rule to get the overall force for the CG pair
- //
- if (strcmp(site1,site2) == 0) fpair = sqrt(fractionOld1_i*fractionOld2_j)*fpairEXP6_12;
- else fpair = sqrt(fractionOld1_i*fractionOld2_j)*fpairEXP6_12 + sqrt(fractionOld2_i*fractionOld1_j)*fpairEXP6_21;
-
- 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 (strcmp(site1,site2) == 0)
- evdwl = sqrt(fraction1_i*fraction2_j)*evdwlEXP6_12;
- else
- evdwl = sqrt(fraction1_i*fraction2_j)*evdwlEXP6_12 + sqrt(fraction2_i*fraction1_j)*evdwlEXP6_21;
- evdwl *= factor_lj;
-
- uCGnew[i] += double(0.5)*evdwl;
- uCGnew[j] += double(0.5)*evdwl;
- evdwl = evdwlOld;
- if (evflag) ev_tally(i,j,nlocal,newton_pair,
- evdwl,0.0,fpair,delx,dely,delz);
- }
+ win1 = buck1*buck2*(rin1*rin1exp*rminv - rm6ij*rin6inv) - rin1*durc;
+
+ aRep = -1.0*win1*pow(rin1,nRep)/nRep;
+
+ uin1rep = aRep/pow(rin1,nRep);
+
+ evdwlEXP6_12 = uin1 - uin1rep + aRep/pow(r,nRep);
+
+ forceExp6 = -1.0*nRep*aRep/pow(r,nRep);
+ fpairEXP6_12 = factor_lj*forceExp6*r2inv;
+
+ } else {
+
+ // A4. Compute the exp-6 force and energy
+ forceExp6 = buck1*buck2*(r*rexp*rminv - rm6ij*r6inv) + r*durc;
+ fpairEXP6_12 = factor_lj*forceExp6*r2inv;
+ evdwlEXP6_12 = buck1*(6.0*rexp - alpha12_ij*rm6ij*r6inv) - urc - durc*(r-rCut);
+ }
+
+ rminv = 1.0/rm21_ij;
+ buck1 = epsilon21_ij / (alpha21_ij - 6.0);
+ buck2 = 6.0*alpha21_ij;
+ rexp = expValue(alpha21_ij*(1.0-r*rminv));
+ rm2ij = rm21_ij*rm21_ij;
+ rm6ij = rm2ij*rm2ij*rm2ij;
+
+ // Compute the shifted potential
+ rCutExp = expValue(alpha21_ij*(1.0-rCut*rminv));
+ urc = buck1*(6.0*rCutExp - alpha21_ij*rm6ij*rCut6inv);
+ durc = -buck1*buck2*(rCutExp*rminv - rCutInv*rm6ij*rCut6inv);
+ rin1 = shift*rm21_ij*func_rin(alpha21_ij);
+
+ if(r < rin1){
+ rin6 = rin1*rin1*rin1*rin1*rin1*rin1;
+ rin6inv = 1.0/rin6;
+
+ rin1exp = expValue(alpha21_ij*(1.0-rin1*rminv));
+
+ uin1 = buck1*(6.0*rin1exp - alpha21_ij*rm6ij*rin6inv) - urc - durc*(rin1-rCut);
+
+ win1 = buck1*buck2*(rin1*rin1exp*rminv - rm6ij*rin6inv) - rin1*durc;
+
+ aRep = -1.0*win1*pow(rin1,nRep)/nRep;
+
+ uin1rep = aRep/pow(rin1,nRep);
+
+ evdwlEXP6_21 = uin1 - uin1rep + aRep/pow(r,nRep);
+
+ forceExp6 = -1.0*nRep*aRep/pow(r,nRep);
+ fpairEXP6_21 = factor_lj*forceExp6*r2inv;
+
+ } else {
+
+ // A4. Compute the exp-6 force and energy
+ forceExp6 = buck1*buck2*(r*rexp*rminv - rm6ij*r6inv) + r*durc;
+ fpairEXP6_21 = factor_lj*forceExp6*r2inv;
+ evdwlEXP6_21 = buck1*(6.0*rexp - alpha21_ij*rm6ij*r6inv) - urc - durc*(r-rCut);
+ }
+
+ //
+ // Apply Mixing Rule to get the overall force for the CG pair
+ //
+ if (strcmp(site1,site2) == 0) fpair = sqrt(fractionOld1_i*fractionOld2_j)*fpairEXP6_12;
+ else fpair = sqrt(fractionOld1_i*fractionOld2_j)*fpairEXP6_12 + sqrt(fractionOld2_i*fractionOld1_j)*fpairEXP6_21;
+
+ 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 (strcmp(site1,site2) == 0)
+ evdwl = sqrt(fraction1_i*fraction2_j)*evdwlEXP6_12;
+ else
+ evdwl = sqrt(fraction1_i*fraction2_j)*evdwlEXP6_12 + sqrt(fraction2_i*fraction1_j)*evdwlEXP6_21;
+ evdwl *= factor_lj;
+
+ uCGnew[i] += 0.5*evdwl;
+ uCGnew[j] += 0.5*evdwl;
+ evdwl = evdwlOld;
+ if (evflag) ev_tally(i,j,nlocal,newton_pair,
+ evdwl,0.0,fpair,delx,dely,delz);
+ }
}
}
}
@@ -451,11 +455,11 @@ void PairExp6rx::coeff(int narg, char **arg)
}
if (ispecies == nspecies && strcmp(site1,"1fluid") != 0)
error->all(FLERR,"Site1 name not recognized in pair coefficients");
-
+
n = strlen(arg[4]) + 1;
site2 = new char[n];
strcpy(site2,arg[4]);
-
+
for (ispecies = 0; ispecies < nspecies; ispecies++){
if (strcmp(site2,&atom->dname[ispecies][0]) == 0) break;
}
@@ -597,8 +601,8 @@ void PairExp6rx::read_file(char *file)
params[nparams].epsilon = atof(words[3]);
params[nparams].rm = atof(words[4]);
if (params[nparams].epsilon <= 0.0 || params[nparams].rm <= 0.0 ||
- params[nparams].alpha < 0.0)
- error->all(FLERR,"Illegal exp6/rx parameters. Rm and Epsilon must be greater than zero. Alpha cannot be negative.");
+ params[nparams].alpha < 0.0)
+ error->all(FLERR,"Illegal exp6/rx parameters. Rm and Epsilon must be greater than zero. Alpha cannot be negative.");
} else {
error->all(FLERR,"Illegal exp6/rx parameters. Interaction potential does not exist.");
}
@@ -624,8 +628,8 @@ void PairExp6rx::setup()
n = -1;
for (j = 0; j < nparams; j++) {
if (i == params[j].ispecies) {
- if (n >= 0) error->all(FLERR,"Potential file has duplicate entry");
- n = j;
+ if (n >= 0) error->all(FLERR,"Potential file has duplicate entry");
+ n = j;
}
}
mol2param[i] = n;
@@ -713,31 +717,29 @@ void PairExp6rx::getParamsEXP6(int id,double &epsilon1,double &alpha1,double &rm
double rmi, rmj, rmij, rm3ij;
double epsiloni, epsilonj, epsilonij;
double alphai, alphaj, alphaij;
- double epsilon_old, rm3_old, alpha_old, fraction_old;
- double epsilon, rm3, alpha, fraction;
+ double epsilon_old, rm3_old, alpha_old;
+ double epsilon, rm3, alpha;
double fractionOFA, fractionOFA_old;
double nTotalOFA, nTotalOFA_old;
double nTotal, nTotal_old;
double xMolei, xMolej, xMolei_old, xMolej_old;
- rm3 = double(0.0);
- epsilon = double(0.0);
- alpha = double(0.0);
- epsilon_old = double(0.0);
- rm3_old = double(0.0);
- alpha_old = double(0.0);
- fraction_old = double(0.0);
- fraction = double(0.0);
- fractionOFA = double(0.0);
- fractionOFA_old = double(0.0);
- nTotalOFA = double(0.0);
- nTotalOFA_old = double(0.0);
- nTotal = double(0.0);
- nTotal_old = double(0.0);
+ rm3 = 0.0;
+ epsilon = 0.0;
+ alpha = 0.0;
+ epsilon_old = 0.0;
+ rm3_old = 0.0;
+ alpha_old = 0.0;
+ fractionOFA = 0.0;
+ fractionOFA_old = 0.0;
+ nTotalOFA = 0.0;
+ nTotalOFA_old = 0.0;
+ nTotal = 0.0;
+ nTotal_old = 0.0;
// Compute the total number of molecules in the old and new CG particle as well as the total number of molecules in the fluid portion of the old and new CG particle
for (int ispecies = 0; ispecies < nspecies; ispecies++){
- nTotal += atom->dvector[ispecies][id];
+ nTotal += atom->dvector[ispecies][id];
nTotal_old += atom->dvector[ispecies+nspecies][id];
iparam = mol2param[ispecies];
@@ -748,8 +750,8 @@ void PairExp6rx::getParamsEXP6(int id,double &epsilon1,double &alpha1,double &rm
nTotalOFA += atom->dvector[ispecies][id];
}
}
- if(nTotal < 1e-8 || nTotal_old < 1e-8)
- error->all(FLERR,"The number of molecules in CG particle is less than 1e-8.");
+ if(nTotal < 1e-8 || nTotal_old < 1e-8)
+ error->all(FLERR,"The number of molecules in CG particle is less than 1e-8.");
// Compute the mole fraction of molecules within the fluid portion of the particle (One Fluid Approximation)
fractionOFA_old = nTotalOFA_old / nTotal_old;
@@ -758,7 +760,7 @@ void PairExp6rx::getParamsEXP6(int id,double &epsilon1,double &alpha1,double &rm
for (int ispecies = 0; ispecies < nspecies; ispecies++) {
iparam = mol2param[ispecies];
if (iparam < 0 || strcmp(params[iparam].potential,"exp6") != 0) continue;
-
+
// If Site1 matches a pure species, then grab the parameters
if (strcmp(site1,params[iparam].name) == 0){
rm1_old = params[iparam].rm;
@@ -771,7 +773,7 @@ void PairExp6rx::getParamsEXP6(int id,double &epsilon1,double &alpha1,double &rm
// Compute the mole fraction of Site1
fraction1_old = atom->dvector[ispecies+nspecies][id]/nTotal_old;
fraction1 = atom->dvector[ispecies][id]/nTotal;
- }
+ }
// If Site2 matches a pure species, then grab the parameters
if (strcmp(site2,params[iparam].name) == 0){
@@ -785,7 +787,7 @@ void PairExp6rx::getParamsEXP6(int id,double &epsilon1,double &alpha1,double &rm
// Compute the mole fraction of Site2
fraction2_old = atom->dvector[ispecies+nspecies][id]/nTotal_old;
fraction2 = atom->dvector[ispecies][id]/nTotal;
- }
+ }
// If Site1 or Site2 matches is a fluid, then compute the paramters
if (strcmp(site1,"1fluid") == 0 || strcmp(site2,"1fluid") == 0) {
@@ -797,30 +799,30 @@ void PairExp6rx::getParamsEXP6(int id,double &epsilon1,double &alpha1,double &rm
xMolei_old = atom->dvector[ispecies+nspecies][id]/nTotalOFA_old;
for (int jspecies = 0; jspecies < nspecies; jspecies++) {
- jparam = mol2param[jspecies];
- if (jparam < 0 || strcmp(params[jparam].potential,"exp6") != 0) continue;
- if (strcmp(site1,params[jparam].name) == 0 || strcmp(site2,params[jparam].name) == 0) continue;
- rmj = params[jparam].rm;
- epsilonj = params[jparam].epsilon;
- alphaj = params[jparam].alpha;
- xMolej = atom->dvector[jspecies][id]/nTotalOFA;
- xMolej_old = atom->dvector[jspecies+nspecies][id]/nTotalOFA_old;
-
- rmij = (rmi+rmj)/2.0;
- rm3ij = rmij*rmij*rmij;
- epsilonij = sqrt(epsiloni*epsilonj);
- alphaij = sqrt(alphai*alphaj);
-
- if(fractionOFA_old > double(0.0)){
- rm3_old += xMolei_old*xMolej_old*rm3ij;
- epsilon_old += xMolei_old*xMolej_old*rm3ij*epsilonij;
- alpha_old += xMolei_old*xMolej_old*rm3ij*epsilonij*alphaij;
- }
- if(fractionOFA > double(0.0)){
- rm3 += xMolei*xMolej*rm3ij;
- epsilon += xMolei*xMolej*rm3ij*epsilonij;
- alpha += xMolei*xMolej*rm3ij*epsilonij*alphaij;
- }
+ jparam = mol2param[jspecies];
+ if (jparam < 0 || strcmp(params[jparam].potential,"exp6") != 0) continue;
+ if (strcmp(site1,params[jparam].name) == 0 || strcmp(site2,params[jparam].name) == 0) continue;
+ rmj = params[jparam].rm;
+ epsilonj = params[jparam].epsilon;
+ alphaj = params[jparam].alpha;
+ xMolej = atom->dvector[jspecies][id]/nTotalOFA;
+ xMolej_old = atom->dvector[jspecies+nspecies][id]/nTotalOFA_old;
+
+ rmij = (rmi+rmj)/2.0;
+ rm3ij = rmij*rmij*rmij;
+ epsilonij = sqrt(epsiloni*epsilonj);
+ alphaij = sqrt(alphai*alphaj);
+
+ if(fractionOFA_old > 0.0){
+ rm3_old += xMolei_old*xMolej_old*rm3ij;
+ epsilon_old += xMolei_old*xMolej_old*rm3ij*epsilonij;
+ alpha_old += xMolei_old*xMolej_old*rm3ij*epsilonij*alphaij;
+ }
+ if(fractionOFA > 0.0){
+ rm3 += xMolei*xMolej*rm3ij;
+ epsilon += xMolei*xMolej*rm3ij*epsilonij;
+ alpha += xMolei*xMolej*rm3ij*epsilonij*alphaij;
+ }
}
}
}
@@ -828,9 +830,9 @@ void PairExp6rx::getParamsEXP6(int id,double &epsilon1,double &alpha1,double &rm
if(strcmp(site1,"1fluid") == 0){
rm1 = cbrt(rm3);
if(rm1 < 1e-16) {
- rm1 = double(0.0);
- epsilon1 = double(0.0);
- alpha1 = double(0.0);
+ rm1 = 0.0;
+ epsilon1 = 0.0;
+ alpha1 = 0.0;
} else {
epsilon1 = epsilon / rm3;
alpha1 = alpha / epsilon1 / rm3;
@@ -840,9 +842,9 @@ void PairExp6rx::getParamsEXP6(int id,double &epsilon1,double &alpha1,double &rm
rm1_old = cbrt(rm3_old);
if(rm1_old < 1e-16) {
- rm1_old = double(0.0);
- epsilon1_old = double(0.0);
- alpha1_old = double(0.0);
+ rm1_old = 0.0;
+ epsilon1_old = 0.0;
+ alpha1_old = 0.0;
} else {
epsilon1_old = epsilon_old / rm3_old;
alpha1_old = alpha_old / epsilon1_old / rm3_old;
@@ -880,14 +882,14 @@ void PairExp6rx::getParamsEXP6(int id,double &epsilon1,double &alpha1,double &rm
rm1 *= pow(nTotalOFA,fuchslinR);
rm1_old *= pow(nTotalOFA_old,fuchslinR);
}
- }
+ }
if(strcmp(site2,"1fluid") == 0){
rm2 = cbrt(rm3);
if(rm2 < 1e-16) {
- rm2 = double(0.0);
- epsilon2 = double(0.0);
- alpha2 = double(0.0);
+ rm2 = 0.0;
+ epsilon2 = 0.0;
+ alpha2 = 0.0;
} else {
epsilon2 = epsilon / rm3;
alpha2 = alpha / epsilon2 / rm3;
@@ -896,9 +898,9 @@ void PairExp6rx::getParamsEXP6(int id,double &epsilon1,double &alpha1,double &rm
rm2_old = cbrt(rm3_old);
if(rm2_old < 1e-16) {
- rm2_old = double(0.0);
- epsilon2_old = double(0.0);
- alpha2_old = double(0.0);
+ rm2_old = 0.0;
+ epsilon2_old = 0.0;
+ alpha2_old = 0.0;
} else {
epsilon2_old = epsilon_old / rm3_old;
alpha2_old = alpha_old / epsilon2_old / rm3_old;
@@ -945,13 +947,11 @@ double PairExp6rx::func_rin(double &alpha)
{
double function;
- const double alpha_min = double(11.0);
- const double alpha_max = double(16.0);
- const double a = double(3.7682065);
- const double b = double(-1.4308614);
-
- function = a+b*pow(alpha,0.5);
- function = expValue(function);
+ const double a = 3.7682065;
+ const double b = -1.4308614;
+
+ function = a+b*sqrt(alpha);
+ function = expValue(function);
return function;
}
@@ -961,8 +961,8 @@ double PairExp6rx::func_rin(double &alpha)
double PairExp6rx::expValue(double value)
{
double returnValue;
- if(value < DBL_MIN_EXP) returnValue = double(0.0);
+ if(value < DBL_MIN_EXP) returnValue = 0.0;
else returnValue = exp(value);
-
+
return returnValue;
}
diff --git a/src/USER-DPD/pair_exp6_rx.h b/src/USER-DPD/pair_exp6_rx.h
index ece29dcda74a32f35287f669b2f130d981decd09..45abfab77aa81fa3966d2d2c5b8b7faad29dcb9a 100644
--- a/src/USER-DPD/pair_exp6_rx.h
+++ b/src/USER-DPD/pair_exp6_rx.h
@@ -1,4 +1,4 @@
-/* ----------------------------------------------------------------------
+/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
@@ -130,7 +130,7 @@ Self-explanatory
E: The number of molecules in CG particle is less than 1e-8.
Self-explanatory. Check the species concentrations have been properly set
-and check the reaction kinetic solver parameters in fix rx to more for
+and check the reaction kinetic solver parameters in fix rx to more for
sufficient accuracy.
diff --git a/src/USER-DPD/pair_multi_lucy.cpp b/src/USER-DPD/pair_multi_lucy.cpp
index 37f235b8e3d527cd4ee41c1c3115990b610e10e8..501ef2d04c51779bc6dd202d1aef40d13fcd3fef 100644
--- a/src/USER-DPD/pair_multi_lucy.cpp
+++ b/src/USER-DPD/pair_multi_lucy.cpp
@@ -12,11 +12,11 @@
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------------------------
- Contributing authors:
+ Contributing authors:
James Larentzos and Joshua Moore (U.S. Army Research Laboratory)
Please cite the related publications:
- J.D. Moore, B.C. Barnes, S. Izvekov, M. Lisal, M.S. Sellers, D.E. Taylor & J.K. Brennan
+ J.D. Moore, B.C. Barnes, S. Izvekov, M. Lisal, M.S. Sellers, D.E. Taylor & J.K. Brennan
"A coarse-grain force field for RDX: Density dependent and energy conserving"
The Journal of Chemical Physics, 2016, 144, 104501.
------------------------------------------------------------------------------------------- */
@@ -87,12 +87,10 @@ PairMultiLucy::~PairMultiLucy()
void PairMultiLucy::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum,itype,jtype,itable;
- double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
- double rsq,factor_lj,fraction,value,a,b;
+ double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair,rsq;
int *ilist,*jlist,*numneigh,**firstneigh;
Table *tb;
- union_int_float_t rsq_lookup;
int tlm1 = tablength - 1;
evdwl = 0.0;
@@ -103,14 +101,12 @@ void PairMultiLucy::compute(int eflag, int vflag)
double **f = atom->f;
int *type = atom->type;
int nlocal = atom->nlocal;
- double *special_lj = force->special_lj;
int newton_pair = force->newton_pair;
double pi = MathConst::MY_PI;
double A_i;
double A_j;
double fraction_i,fraction_j;
- double a_i,a_j,b_i,b_j;
int jtable;
double *rho = atom->rho;
@@ -134,7 +130,6 @@ void PairMultiLucy::compute(int eflag, int vflag)
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
- factor_lj = special_lj[sbmask(j)];
j &= NEIGHMASK;
delx = xtmp - x[j][0];
@@ -151,19 +146,19 @@ void PairMultiLucy::compute(int eflag, int vflag)
if (tabstyle == LOOKUP) {
itable = static_cast<int> (((rho[i]*rho[i]) - tb->innersq) * tb->invdelta);
jtable = static_cast<int> (((rho[j]*rho[j]) - tb->innersq) * tb->invdelta);
- if (itable >= tlm1 || jtable >= tlm1)
- error->one(FLERR,"Density > table outer cutoff");
+ if (itable >= tlm1 || jtable >= tlm1)
+ error->one(FLERR,"Density > table outer cutoff");
A_i = tb->f[itable];
A_j = tb->f[jtable];
- fpair = double(0.5)*(A_i + A_j)*(double(1.0)+double(3.0)*sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(double(1.0) - sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(double(1.0) - sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(double(1.0) - sqrt(rsq)/sqrt(cutsq[itype][jtype]));
+ fpair = 0.5*(A_i + A_j)*(1.0+3.0*sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(1.0 - sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(1.0 - sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(1.0 - sqrt(rsq)/sqrt(cutsq[itype][jtype]));
fpair = fpair/sqrt(rsq);
} else if (tabstyle == LINEAR) {
- itable = static_cast<int> ((rho[i]*rho[i] - tb->innersq) * tb->invdelta);
+ itable = static_cast<int> ((rho[i]*rho[i] - tb->innersq) * tb->invdelta);
jtable = static_cast<int> (((rho[j]*rho[j]) - tb->innersq) * tb->invdelta);
- if (itable >= tlm1 || jtable >= tlm1)
- error->one(FLERR,"Density > table outer cutoff");
+ if (itable >= tlm1 || jtable >= tlm1)
+ error->one(FLERR,"Density > table outer cutoff");
fraction_i = (((rho[i]*rho[i]) - tb->rsq[itable]) * tb->invdelta);
fraction_j = (((rho[j]*rho[j]) - tb->rsq[jtable]) * tb->invdelta);
@@ -171,7 +166,7 @@ void PairMultiLucy::compute(int eflag, int vflag)
A_i = tb->f[itable] + fraction_i*tb->df[itable];
A_j = tb->f[jtable] + fraction_j*tb->df[jtable];
- fpair = double(0.5)*(A_i + A_j)*(double(1.0)+double(3.0)*sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(double(1.0) - sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(double(1.0) - sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(double(1.0) - sqrt(rsq)/sqrt(cutsq[itype][jtype]));
+ fpair = 0.5*(A_i + A_j)*(1.0+3.0*sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(1.0 - sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(1.0 - sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(1.0 - sqrt(rsq)/sqrt(cutsq[itype][jtype]));
fpair = fpair / sqrt(rsq);
@@ -185,8 +180,8 @@ void PairMultiLucy::compute(int eflag, int vflag)
f[j][1] -= dely*fpair;
f[j][2] -= delz*fpair;
}
- if (evflag) ev_tally(i,j,nlocal,newton_pair,
- 0.0,0.0,fpair,delx,dely,delz);
+ if (evflag) ev_tally(i,j,nlocal,newton_pair,
+ 0.0,0.0,fpair,delx,dely,delz);
}
}
@@ -195,17 +190,17 @@ void PairMultiLucy::compute(int eflag, int vflag)
error->one(FLERR,"Density < table inner cutoff");
itable = static_cast<int> (((rho[i]*rho[i]) - tb->innersq) * tb->invdelta);
if (tabstyle == LOOKUP) evdwl = tb->e[itable];
- else if (tabstyle == LINEAR){
+ else if (tabstyle == LINEAR){
if (itable >= tlm1) error->one(FLERR,"Density > table outer cutoff");
- if(itable==0) fraction_i=double(0.0);
+ if(itable==0) fraction_i=0.0;
else fraction_i = (((rho[i]*rho[i]) - tb->rsq[itable]) * tb->invdelta);
evdwl = tb->e[itable] + fraction_i*tb->de[itable];
} else error->one(FLERR,"Only LOOKUP and LINEAR table styles have been implemented for pair multi/lucy");
- evdwl *=(pi*cutsq[itype][itype]*cutsq[itype][itype])/double(84.0);
+ evdwl *=(pi*cutsq[itype][itype]*cutsq[itype][itype])/84.0;
if (evflag) ev_tally(0,0,nlocal,newton_pair,
- evdwl,0.0,0.0,0.0,0.0,0.0);
+ evdwl,0.0,0.0,0.0,0.0,0.0);
}
if (vflag_fdotr) virial_fdotr_compute();
@@ -293,13 +288,11 @@ void PairMultiLucy::coeff(int narg, char **arg)
// insure cutoff is within table
if (tb->ninput <= 1) error->one(FLERR,"Invalid pair table length");
- double rlo,rhi;
+ double rlo;
if (tb->rflag == 0) {
rlo = tb->rfile[0];
- rhi = tb->rfile[tb->ninput-1];
} else {
rlo = tb->rlo;
- rhi = tb->rhi;
}
rho_0 = rlo;
@@ -390,7 +383,6 @@ void PairMultiLucy::read_table(Table *tb, char *file, char *keyword)
int itmp;
double rtmp;
- union_int_float_t rsq_lookup;
fgets(line,MAXLINE,fp);
for (int i = 0; i < tb->ninput; i++) {
@@ -733,7 +725,7 @@ void PairMultiLucy::computeLocalDensity()
int *type = atom->type;
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
- double factor, volume;
+ double factor;
inum = list->inum;
ilist = list->ilist;
@@ -742,9 +734,9 @@ void PairMultiLucy::computeLocalDensity()
double pi = MathConst::MY_PI;
double *rho = atom->rho;
-
+
// zero out density
-
+
if (newton_pair) {
m = nlocal + atom->nghost;
for (i = 0; i < m; i++) rho[i] = 0.0;
@@ -773,7 +765,7 @@ void PairMultiLucy::computeLocalDensity()
if (rsq < cutsq[itype][jtype]) {
double rcut = sqrt(cutsq[itype][jtype]);
- factor= (double(84.0)/(double(5.0)*pi*rcut*rcut*rcut))*(double(1.0)+double(3.0)*sqrt(rsq)/(double(2.0)*rcut))*(double(1.0)-sqrt(rsq)/rcut)*(double(1.0)-sqrt(rsq)/rcut)*(double(1.0)-sqrt(rsq)/rcut)*(double(1.0)-sqrt(rsq)/rcut);
+ factor= (84.0/(5.0*pi*rcut*rcut*rcut))*(1.0+3.0*sqrt(rsq)/(2.0*rcut))*(1.0-sqrt(rsq)/rcut)*(1.0-sqrt(rsq)/rcut)*(1.0-sqrt(rsq)/rcut)*(1.0-sqrt(rsq)/rcut);
rho[i] += factor;
if (newton_pair || j < nlocal) {
rho[j] += factor;
diff --git a/src/USER-DPD/pair_multi_lucy.h b/src/USER-DPD/pair_multi_lucy.h
index ee6f8c443178144619ac86057f01e5a730d0aed8..f3c67e4fa41124a824a2c7e846c6e1b61b33df46 100644
--- a/src/USER-DPD/pair_multi_lucy.h
+++ b/src/USER-DPD/pair_multi_lucy.h
@@ -1,4 +1,4 @@
-/* ----------------------------------------------------------------------
+/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
diff --git a/src/USER-DPD/pair_multi_lucy_rx.cpp b/src/USER-DPD/pair_multi_lucy_rx.cpp
index 6a349101d5d8f77530bbdc8ea12a0edfbe519723..b90252f66f3bbec14eca254bfb4741255e36456e 100644
--- a/src/USER-DPD/pair_multi_lucy_rx.cpp
+++ b/src/USER-DPD/pair_multi_lucy_rx.cpp
@@ -12,15 +12,15 @@
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------------------------
- Contributing authors:
+ Contributing authors:
James Larentzos and Joshua Moore (U.S. Army Research Laboratory)
Please cite the related publications:
- J.D. Moore, B.C. Barnes, S. Izvekov, M. Lisal, M.S. Sellers, D.E. Taylor & J.K. Brennan
+ J.D. Moore, B.C. Barnes, S. Izvekov, M. Lisal, M.S. Sellers, D.E. Taylor & J.K. Brennan
"A coarse-grain force field for RDX: Density dependent and energy conserving"
The Journal of Chemical Physics, 2016, 144, 104501.
------------------------------------------------------------------------------------------- */
-
+
#include "mpi.h"
#include <math.h>
#include "math_const.h"
@@ -90,7 +90,7 @@ void PairMultiLucyRX::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum,itype,jtype,itable;
double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,evdwlOld,fpair;
- double rsq,factor_lj;
+ double rsq;
int *ilist,*jlist,*numneigh,**firstneigh;
Table *tb;
@@ -105,7 +105,6 @@ void PairMultiLucyRX::compute(int eflag, int vflag)
double **f = atom->f;
int *type = atom->type;
int nlocal = atom->nlocal;
- double *special_lj = force->special_lj;
int newton_pair = force->newton_pair;
double fractionOld1_i,fractionOld1_j;
@@ -143,7 +142,6 @@ void PairMultiLucyRX::compute(int eflag, int vflag)
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
- factor_lj = special_lj[sbmask(j)];
j &= NEIGHMASK;
delx = xtmp - x[j][0];
@@ -153,63 +151,63 @@ void PairMultiLucyRX::compute(int eflag, int vflag)
jtype = type[j];
if (rsq < cutsq[itype][jtype]) {
- fpair = double(0.0);
- getParams(j,fractionOld1_j,fractionOld2_j,fraction1_j,fraction2_j);
+ fpair = 0.0;
+ getParams(j,fractionOld1_j,fractionOld2_j,fraction1_j,fraction2_j);
tb = &tables[tabindex[itype][jtype]];
if (rho[i]*rho[i] < tb->innersq || rho[j]*rho[j] < tb->innersq){
- printf("Table inner cutoff = %lf\n",sqrt(tb->innersq));
- printf("rho[%d]=%lf\n",i,rho[i]);
- printf("rho[%d]=%lf\n",j,rho[j]);
+ printf("Table inner cutoff = %lf\n",sqrt(tb->innersq));
+ printf("rho[%d]=%lf\n",i,rho[i]);
+ printf("rho[%d]=%lf\n",j,rho[j]);
error->one(FLERR,"Density < table inner cutoff");
- }
+ }
if (tabstyle == LOOKUP) {
- itable = static_cast<int> (((rho[i]*rho[i]) - tb->innersq) * tb->invdelta);
- jtable = static_cast<int> (((rho[j]*rho[j]) - tb->innersq) * tb->invdelta);
- if (itable >= tlm1 || jtable >= tlm1){
- printf("Table outer index = %d\n",tlm1);
- printf("itableIndex=%d rho[%d]=%lf\n",itable,i,rho[i]);
- printf("jtableIndex=%d rho[%d]=%lf\n",jtable,j,rho[j]);
- error->one(FLERR,"Density > table outer cutoff");
- }
+ itable = static_cast<int> (((rho[i]*rho[i]) - tb->innersq) * tb->invdelta);
+ jtable = static_cast<int> (((rho[j]*rho[j]) - tb->innersq) * tb->invdelta);
+ if (itable >= tlm1 || jtable >= tlm1){
+ printf("Table outer index = %d\n",tlm1);
+ printf("itableIndex=%d rho[%d]=%lf\n",itable,i,rho[i]);
+ printf("jtableIndex=%d rho[%d]=%lf\n",jtable,j,rho[j]);
+ error->one(FLERR,"Density > table outer cutoff");
+ }
A_i = tb->f[itable];
A_j = tb->f[jtable];
- fpair = double(0.5)*(A_i + A_j)*(double(1.0)+double(3.0)*sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(double(1.0) - sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(double(1.0) - sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(double(1.0) - sqrt(rsq)/sqrt(cutsq[itype][jtype]));
+ fpair = 0.5*(A_i + A_j)*(1.0+3.0*sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(1.0 - sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(1.0 - sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(1.0 - sqrt(rsq)/sqrt(cutsq[itype][jtype]));
fpair = fpair/sqrt(rsq);
} else if (tabstyle == LINEAR) {
- itable = static_cast<int> ((rho[i]*rho[i] - tb->innersq) * tb->invdelta);
- jtable = static_cast<int> (((rho[j]*rho[j]) - tb->innersq) * tb->invdelta);
- if (itable >= tlm1 || jtable >= tlm1){
- printf("Table outer index = %d\n",tlm1);
- printf("itableIndex=%d rho[%d]=%lf\n",itable,i,rho[i]);
- printf("jtableIndex=%d rho[%d]=%lf\n",jtable,j,rho[j]);
- error->one(FLERR,"Density > table outer cutoff");
- }
- if(itable<0) itable=0;
- if(itable>=tlm1) itable=tlm1;
- if(jtable<0) jtable=0;
- if(jtable>=tlm1)jtable=tlm1;
+ itable = static_cast<int> ((rho[i]*rho[i] - tb->innersq) * tb->invdelta);
+ jtable = static_cast<int> (((rho[j]*rho[j]) - tb->innersq) * tb->invdelta);
+ if (itable >= tlm1 || jtable >= tlm1){
+ printf("Table outer index = %d\n",tlm1);
+ printf("itableIndex=%d rho[%d]=%lf\n",itable,i,rho[i]);
+ printf("jtableIndex=%d rho[%d]=%lf\n",jtable,j,rho[j]);
+ error->one(FLERR,"Density > table outer cutoff");
+ }
+ if(itable<0) itable=0;
+ if(itable>=tlm1) itable=tlm1;
+ if(jtable<0) jtable=0;
+ if(jtable>=tlm1)jtable=tlm1;
fraction_i = (((rho[i]*rho[i]) - tb->rsq[itable]) * tb->invdelta);
fraction_j = (((rho[j]*rho[j]) - tb->rsq[jtable]) * tb->invdelta);
- if(itable==0) fraction_i=double(0.0);
- if(itable==tlm1) fraction_i=double(0.0);
- if(jtable==0) fraction_j=double(0.0);
- if(jtable==tlm1) fraction_j=double(0.0);
+ if(itable==0) fraction_i=0.0;
+ if(itable==tlm1) fraction_i=0.0;
+ if(jtable==0) fraction_j=0.0;
+ if(jtable==tlm1) fraction_j=0.0;
A_i = tb->f[itable] + fraction_i*tb->df[itable];
A_j = tb->f[jtable] + fraction_j*tb->df[jtable];
- fpair = double(0.5)*(A_i + A_j)*(double(1.0)+double(3.0)*sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(double(1.0) - sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(double(1.0) - sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(double(1.0) - sqrt(rsq)/sqrt(cutsq[itype][jtype]));
+ fpair = 0.5*(A_i + A_j)*(1.0+3.0*sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(1.0 - sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(1.0 - sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(1.0 - sqrt(rsq)/sqrt(cutsq[itype][jtype]));
fpair = fpair / sqrt(rsq);
} else error->one(FLERR,"Only LOOKUP and LINEAR table styles have been implemented for pair multi/lucy/rx");
-
- if (strcmp(site1,site2) == 0) fpair = sqrt(fractionOld1_i*fractionOld2_j)*fpair;
- else fpair = (sqrt(fractionOld1_i*fractionOld2_j) + sqrt(fractionOld2_i*fractionOld1_j))*fpair;
+
+ if (strcmp(site1,site2) == 0) fpair = sqrt(fractionOld1_i*fractionOld2_j)*fpair;
+ else fpair = (sqrt(fractionOld1_i*fractionOld2_j) + sqrt(fractionOld2_i*fractionOld1_j))*fpair;
f[i][0] += delx*fpair;
f[i][1] += dely*fpair;
@@ -219,35 +217,35 @@ void PairMultiLucyRX::compute(int eflag, int vflag)
f[j][1] -= dely*fpair;
f[j][2] -= delz*fpair;
}
- if (evflag) ev_tally(i,j,nlocal,newton_pair,
- 0.0,0.0,fpair,delx,dely,delz);
+ if (evflag) ev_tally(i,j,nlocal,newton_pair,
+ 0.0,0.0,fpair,delx,dely,delz);
}
}
tb = &tables[tabindex[itype][itype]];
itable = static_cast<int> (((rho[i]*rho[i]) - tb->innersq) * tb->invdelta);
if (tabstyle == LOOKUP) evdwl = tb->e[itable];
- else if (tabstyle == LINEAR){
+ else if (tabstyle == LINEAR){
if (itable >= tlm1){
- printf("itableIndex=%d rho[%d]=%lf\n",itable,i,rho[i]);
- error->one(FLERR,"Density > table outer cutoff");
+ printf("itableIndex=%d rho[%d]=%lf\n",itable,i,rho[i]);
+ error->one(FLERR,"Density > table outer cutoff");
}
- if(itable==0) fraction_i=double(0.0);
+ if(itable==0) fraction_i=0.0;
else fraction_i = (((rho[i]*rho[i]) - tb->rsq[itable]) * tb->invdelta);
evdwl = tb->e[itable] + fraction_i*tb->de[itable];
} else error->one(FLERR,"Only LOOKUP and LINEAR table styles have been implemented for pair multi/lucy/rx");
-
- evdwl *=(pi*cutsq[itype][itype]*cutsq[itype][itype])/double(84.0);
+
+ evdwl *=(pi*cutsq[itype][itype]*cutsq[itype][itype])/84.0;
evdwlOld = fractionOld1_i*evdwl;
evdwl = fraction1_i*evdwl;
-
+
uCG[i] += evdwlOld;
uCGnew[i] += evdwl;
evdwl = evdwlOld;
if (evflag) ev_tally(0,0,nlocal,newton_pair,
- evdwl,0.0,0.0,0.0,0.0,0.0);
+ evdwl,0.0,0.0,0.0,0.0,0.0);
}
if (vflag_fdotr) virial_fdotr_compute();
@@ -337,7 +335,7 @@ void PairMultiLucyRX::coeff(int narg, char **arg)
n = strlen(arg[3]) + 1;
site1 = new char[n];
strcpy(site1,arg[4]);
-
+
n = strlen(arg[4]) + 1;
site2 = new char[n];
strcpy(site2,arg[5]);
@@ -352,13 +350,11 @@ void PairMultiLucyRX::coeff(int narg, char **arg)
// insure cutoff is within table
if (tb->ninput <= 1) error->one(FLERR,"Invalid pair table length");
- double rlo,rhi;
+ double rlo;
if (tb->rflag == 0) {
rlo = tb->rfile[0];
- rhi = tb->rfile[tb->ninput-1];
} else {
rlo = tb->rlo;
- rhi = tb->rhi;
}
rho_0 = rlo;
@@ -412,7 +408,7 @@ void PairMultiLucyRX::read_table(Table *tb, char *file, char *keyword)
// open file
- FILE *fp = fopen(file,"r");
+ FILE *fp = force->open_potential(file);
if (fp == NULL) {
char str[128];
sprintf(str,"Cannot open file %s",file);
@@ -461,7 +457,7 @@ void PairMultiLucyRX::read_table(Table *tb, char *file, char *keyword)
rtmp = tb->rlo*tb->rlo +
(tb->rhi*tb->rhi - tb->rlo*tb->rlo)*i/(tb->ninput-1);
rtmp = sqrt(rtmp);
- }
+ }
tb->rfile[i] = rtmp;
}
@@ -800,9 +796,9 @@ void PairMultiLucyRX::computeLocalDensity()
double pi = MathConst::MY_PI;
double *rho = atom->rho;
-
+
// zero out density
-
+
if (newton_pair) {
m = nlocal + atom->nghost;
for (i = 0; i < m; i++) rho[i] = 0.0;
@@ -831,9 +827,9 @@ void PairMultiLucyRX::computeLocalDensity()
if (rsq < cutsq[itype][jtype]) {
double rcut = sqrt(cutsq[itype][jtype]);
- double tmpFactor = double(1.0)-sqrt(rsq)/rcut;
- double tmpFactor4 = tmpFactor*tmpFactor*tmpFactor*tmpFactor;
- factor = (double(84.0)/(double(5.0)*pi*rcut*rcut*rcut))*(double(1.0)+double(3.0)*sqrt(rsq)/(double(2.0)*rcut))*tmpFactor4;
+ double tmpFactor = 1.0-sqrt(rsq)/rcut;
+ double tmpFactor4 = tmpFactor*tmpFactor*tmpFactor*tmpFactor;
+ factor = (84.0/(5.0*pi*rcut*rcut*rcut))*(1.0+3.0*sqrt(rsq)/(2.0*rcut))*tmpFactor4;
rho[i] += factor;
if (newton_pair || j < nlocal) {
rho[j] += factor;
@@ -854,17 +850,17 @@ void PairMultiLucyRX::getParams(int id, double &fractionOld1, double &fractionOl
double fractionOld, fraction;
double nTotal, nTotalOld;
- fractionOld = double(0.0);
- fraction = double(0.0);
- fractionOld1 = double(0.0);
- fractionOld2 = double(0.0);
- fraction1 = double(0.0);
- fraction2 = double(0.0);
+ fractionOld = 0.0;
+ fraction = 0.0;
+ fractionOld1 = 0.0;
+ fractionOld2 = 0.0;
+ fraction1 = 0.0;
+ fraction2 = 0.0;
- nTotal = double(0.0);
- nTotalOld = double(0.0);
+ nTotal = 0.0;
+ nTotalOld = 0.0;
for(int ispecies=0;ispecies<nspecies;ispecies++){
- nTotal += atom->dvector[ispecies][id];
+ nTotal += atom->dvector[ispecies][id];
nTotalOld += atom->dvector[ispecies+nspecies][id];
}
@@ -880,7 +876,7 @@ void PairMultiLucyRX::getParams(int id, double &fractionOld1, double &fractionOl
}
for (int ispecies = 0; ispecies < nspecies; ispecies++) {
- if (strcmp(site1,&atom->dname[ispecies][0]) == 0){
+ if (strcmp(site1,&atom->dname[ispecies][0]) == 0){
fractionOld1 = atom->dvector[ispecies+nspecies][id]/nTotalOld;
fraction1 = atom->dvector[ispecies][id]/nTotal;
}
diff --git a/src/USER-DPD/pair_multi_lucy_rx.h b/src/USER-DPD/pair_multi_lucy_rx.h
index 8600a3d2c67e3d95fa273cda04fdb5c88abad268..a85bc6cd75f1955c8534cd45b213e6990ad3b122 100644
--- a/src/USER-DPD/pair_multi_lucy_rx.h
+++ b/src/USER-DPD/pair_multi_lucy_rx.h
@@ -1,4 +1,4 @@
-/* ----------------------------------------------------------------------
+/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
diff --git a/src/USER-DPD/pair_table_rx.cpp b/src/USER-DPD/pair_table_rx.cpp
index 71b48b277677dbce30a65299a1d9bed8f9da1475..6b824b415a71e99d268de7b74f0c918eed65a839 100644
--- a/src/USER-DPD/pair_table_rx.cpp
+++ b/src/USER-DPD/pair_table_rx.cpp
@@ -70,9 +70,9 @@ void PairTableRX::compute(int eflag, int vflag)
union_int_float_t rsq_lookup;
int tlm1 = tablength - 1;
- fraction = double(0.0);
- a = double(0.0);
- b = double(0.0);
+ fraction = 0.0;
+ a = 0.0;
+ b = 0.0;
evdwlOld = 0.0;
evdwl = 0.0;
@@ -122,7 +122,7 @@ void PairTableRX::compute(int eflag, int vflag)
jtype = type[j];
if (rsq < cutsq[itype][jtype]) {
- getParams(j,fractionOld1_j,fractionOld2_j,fraction1_j,fraction2_j);
+ getParams(j,fractionOld1_j,fractionOld2_j,fraction1_j,fraction2_j);
tb = &tables[tabindex[itype][jtype]];
if (rsq < tb->innersq)
@@ -158,8 +158,8 @@ void PairTableRX::compute(int eflag, int vflag)
value = tb->f[itable] + fraction*tb->df[itable];
fpair = factor_lj * value;
}
- if (strcmp(site1,site2) == 0) fpair = sqrt(fractionOld1_i*fractionOld2_j)*fpair;
- else fpair = (sqrt(fractionOld1_i*fractionOld2_j) + sqrt(fractionOld2_i*fractionOld1_j))*fpair;
+ if (strcmp(site1,site2) == 0) fpair = sqrt(fractionOld1_i*fractionOld2_j)*fpair;
+ else fpair = (sqrt(fractionOld1_i*fractionOld2_j) + sqrt(fractionOld2_i*fractionOld1_j))*fpair;
f[i][0] += delx*fpair;
f[i][1] += dely*fpair;
@@ -170,29 +170,29 @@ void PairTableRX::compute(int eflag, int vflag)
f[j][2] -= delz*fpair;
}
- if (tabstyle == LOOKUP)
- evdwl = tb->e[itable];
- else if (tabstyle == LINEAR || tabstyle == BITMAP){
- evdwl = tb->e[itable] + fraction*tb->de[itable];
- }
- else
- evdwl = a * tb->e[itable] + b * tb->e[itable+1] +
- ((a*a*a-a)*tb->e2[itable] + (b*b*b-b)*tb->e2[itable+1]) *
- tb->deltasq6;
- if (strcmp(site1,site2) == 0){
- evdwlOld = sqrt(fractionOld1_i*fractionOld2_j)*evdwl;
- evdwl = sqrt(fraction1_i*fraction2_j)*evdwl;
- } else {
- evdwlOld = (sqrt(fractionOld1_i*fractionOld2_j) + sqrt(fractionOld2_i*fractionOld1_j))*evdwl;
- evdwl = (sqrt(fraction1_i*fraction2_j) + sqrt(fraction2_i*fraction1_j))*evdwl;
- }
- evdwlOld *= factor_lj;
- evdwl *= factor_lj;
- uCG[i] += double(0.5)*evdwlOld;
- uCG[j] += double(0.5)*evdwlOld;
- uCGnew[i] += double(0.5)*evdwl;
- uCGnew[j] += double(0.5)*evdwl;
- evdwl = evdwlOld;
+ if (tabstyle == LOOKUP)
+ evdwl = tb->e[itable];
+ else if (tabstyle == LINEAR || tabstyle == BITMAP){
+ evdwl = tb->e[itable] + fraction*tb->de[itable];
+ }
+ else
+ evdwl = a * tb->e[itable] + b * tb->e[itable+1] +
+ ((a*a*a-a)*tb->e2[itable] + (b*b*b-b)*tb->e2[itable+1]) *
+ tb->deltasq6;
+ if (strcmp(site1,site2) == 0){
+ evdwlOld = sqrt(fractionOld1_i*fractionOld2_j)*evdwl;
+ evdwl = sqrt(fraction1_i*fraction2_j)*evdwl;
+ } else {
+ evdwlOld = (sqrt(fractionOld1_i*fractionOld2_j) + sqrt(fractionOld2_i*fractionOld1_j))*evdwl;
+ evdwl = (sqrt(fraction1_i*fraction2_j) + sqrt(fraction2_i*fraction1_j))*evdwl;
+ }
+ evdwlOld *= factor_lj;
+ evdwl *= factor_lj;
+ uCG[i] += 0.5*evdwlOld;
+ uCG[j] += 0.5*evdwlOld;
+ uCGnew[i] += 0.5*evdwl;
+ uCGnew[j] += 0.5*evdwl;
+ evdwl = evdwlOld;
if (evflag) ev_tally(i,j,nlocal,newton_pair,
evdwl,0.0,fpair,delx,dely,delz);
@@ -309,7 +309,7 @@ void PairTableRX::coeff(int narg, char **arg)
}
if (ispecies == nspecies && strcmp(site1,"1fluid") != 0)
error->all(FLERR,"Site1 name not recognized in pair coefficients");
-
+
n = strlen(arg[4]) + 1;
site2 = new char[n];
strcpy(site2,arg[5]);
@@ -599,7 +599,7 @@ void PairTableRX::compute_table(Table *tb)
else inner = tb->rfile[0];
tb->innersq = double(inner)*double(inner);
tb->delta = double(tb->cut*tb->cut - double(tb->innersq)) / double(tlm1);
- tb->invdelta = double(1.0)/double(tb->delta);
+ tb->invdelta = 1.0/double(tb->delta);
// direct lookup tables
// N-1 evenly spaced bins in rsq from inner to cut
@@ -988,9 +988,9 @@ double PairTableRX::single(int i, int j, int itype, int jtype, double rsq,
double fractionOld1_i, fractionOld1_j;
double fractionOld2_i, fractionOld2_j;
- fraction = double(0.0);
- a = double(0.0);
- b = double(0.0);
+ fraction = 0.0;
+ a = 0.0;
+ b = 0.0;
getParams(i,fractionOld1_i,fractionOld2_i,fraction1_i,fraction2_i);
getParams(j,fractionOld1_j,fractionOld2_j,fraction1_j,fraction2_j);
@@ -1026,7 +1026,7 @@ double PairTableRX::single(int i, int j, int itype, int jtype, double rsq,
fforce = factor_lj * value;
}
- if (strcmp(site1,site2) == 0) fforce = sqrt(fraction1_i*fraction2_j)*fforce;
+ if (strcmp(site1,site2) == 0) fforce = sqrt(fraction1_i*fraction2_j)*fforce;
else fforce = (sqrt(fraction1_i*fraction2_j) + sqrt(fraction2_i*fraction1_j))*fforce;
if (tabstyle == LOOKUP)
@@ -1071,39 +1071,39 @@ void PairTableRX::getParams(int id, double &fractionOld1, double &fractionOld2,
double fractionOld, fraction;
double nTotal, nTotalOld;
- fractionOld = double(0.0);
- fraction = double(0.0);
- fractionOld1 = double(0.0);
- fractionOld2 = double(0.0);
- fraction1 = double(0.0);
- fraction2 = double(0.0);
+ fractionOld = 0.0;
+ fraction = 0.0;
+ fractionOld1 = 0.0;
+ fractionOld2 = 0.0;
+ fraction1 = 0.0;
+ fraction2 = 0.0;
- nTotal = double(0.0);
- nTotalOld = double(0.0);
+ nTotal = 0.0;
+ nTotalOld = 0.0;
for(int ispecies=0;ispecies<nspecies;ispecies++){
- nTotal += atom->dvector[ispecies][id];
+ nTotal += atom->dvector[ispecies][id];
nTotalOld += atom->dvector[ispecies+nspecies][id];
}
- if(nTotal < 1e-8 || nTotalOld < 1e-8)
- error->all(FLERR,"The number of molecules in CG particle is less than 1e-8.");
+ if(nTotal < 1e-8 || nTotalOld < 1e-8)
+ error->all(FLERR,"The number of molecules in CG particle is less than 1e-8.");
for (int ispecies = 0; ispecies < nspecies; ispecies++) {
- if (strcmp(site1,&atom->dname[ispecies][0]) == 0){
+ if (strcmp(site1,&atom->dname[ispecies][0]) == 0){
fractionOld1 = atom->dvector[ispecies+nspecies][id]/nTotalOld;
fraction1 = atom->dvector[ispecies][id]/nTotal;
}
- if (strcmp(site2,&atom->dname[ispecies][0]) == 0){
+ if (strcmp(site2,&atom->dname[ispecies][0]) == 0){
fractionOld2 = atom->dvector[ispecies+nspecies][id]/nTotalOld;
fraction2 = atom->dvector[ispecies][id]/nTotal;
}
}
for (int ispecies = 0; ispecies < nspecies; ispecies++) {
- if (strcmp(site1,&atom->dname[ispecies][0]) == 0){
+ if (strcmp(site1,&atom->dname[ispecies][0]) == 0){
fractionOld1 = atom->dvector[ispecies+nspecies][id]/nTotalOld;
fraction1 = atom->dvector[ispecies][id]/nTotal;
}
- if (strcmp(site2,&atom->dname[ispecies][0]) == 0){
+ if (strcmp(site2,&atom->dname[ispecies][0]) == 0){
fractionOld2 = atom->dvector[ispecies+nspecies][id]/nTotalOld;
fraction2 = atom->dvector[ispecies][id]/nTotal;
}
diff --git a/src/USER-DPD/pair_table_rx.h b/src/USER-DPD/pair_table_rx.h
index 00b83bafb7abda9a7f4ec845035ee22aa04f5f36..42e27390dea2b500387c874c658ac5d1a53235e7 100644
--- a/src/USER-DPD/pair_table_rx.h
+++ b/src/USER-DPD/pair_table_rx.h
@@ -1,4 +1,4 @@
-/* ----------------------------------------------------------------------
+/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
@@ -166,7 +166,7 @@ long-range solver starts at that cutoff.
E: The number of molecules in CG particle is less than 1e-8
Self-explanatory. Check the species concentrations have been properly set
-and check the reaction kinetic solver parameters in fix rx to more for
+and check the reaction kinetic solver parameters in fix rx to more for
sufficient accuracy.
*/
diff --git a/src/USER-INTEL/improper_cvff_intel.cpp b/src/USER-INTEL/improper_cvff_intel.cpp
index 6c0735a75504f501621fc5c7a90f6e4ce3583d40..0fb02420b92684c2e3a6fd0277d574b85d5e7461 100644
--- a/src/USER-INTEL/improper_cvff_intel.cpp
+++ b/src/USER-INTEL/improper_cvff_intel.cpp
@@ -364,10 +364,8 @@ void ImproperCvffIntel::eval(const int vflag,
/* ---------------------------------------------------------------------- */
-void ImproperCvffIntel::init()
+void ImproperCvffIntel::init_style()
{
- ImproperCvff::init();
-
int ifix = modify->find_fix("package_intel");
if (ifix < 0)
error->all(FLERR,
diff --git a/src/USER-INTEL/improper_cvff_intel.h b/src/USER-INTEL/improper_cvff_intel.h
index 154bb5bd1879f972ac0dc18a89bf19d09368b429..95ccd8f9d212c8ebd4c174ec3bd376086b752a52 100644
--- a/src/USER-INTEL/improper_cvff_intel.h
+++ b/src/USER-INTEL/improper_cvff_intel.h
@@ -35,7 +35,7 @@ class ImproperCvffIntel : public ImproperCvff {
ImproperCvffIntel(class LAMMPS *);
virtual ~ImproperCvffIntel();
virtual void compute(int, int);
- virtual void init();
+ virtual void init_style();
protected:
FixIntel *fix;
diff --git a/src/USER-INTEL/improper_harmonic_intel.cpp b/src/USER-INTEL/improper_harmonic_intel.cpp
index 7693793f02228c837e0231dde9b24b987d7a6e56..071ff548ea8419fd90c8c1c3a235e0a53856e412 100644
--- a/src/USER-INTEL/improper_harmonic_intel.cpp
+++ b/src/USER-INTEL/improper_harmonic_intel.cpp
@@ -325,10 +325,8 @@ void ImproperHarmonicIntel::eval(const int vflag,
/* ---------------------------------------------------------------------- */
-void ImproperHarmonicIntel::init()
+void ImproperHarmonicIntel::init_style()
{
- ImproperHarmonic::init();
-
int ifix = modify->find_fix("package_intel");
if (ifix < 0)
error->all(FLERR,
diff --git a/src/USER-INTEL/improper_harmonic_intel.h b/src/USER-INTEL/improper_harmonic_intel.h
index d5b0e4ee0fc7709d7f2944bf5164f42fe7b9b6aa..4e3838386347809afd3b228da7d9fc57151105e8 100644
--- a/src/USER-INTEL/improper_harmonic_intel.h
+++ b/src/USER-INTEL/improper_harmonic_intel.h
@@ -35,7 +35,7 @@ class ImproperHarmonicIntel : public ImproperHarmonic {
ImproperHarmonicIntel(class LAMMPS *);
virtual ~ImproperHarmonicIntel();
virtual void compute(int, int);
- virtual void init();
+ virtual void init_style();
protected:
FixIntel *fix;
diff --git a/src/USER-MISC/fix_addtorque.cpp b/src/USER-MISC/fix_addtorque.cpp
index 0e1ffaf0e1d7cd5400de486d931cad65193d2b64..ed5d8875f7c37d0b1fb57ce2d2136c42f73e8083 100644
--- a/src/USER-MISC/fix_addtorque.cpp
+++ b/src/USER-MISC/fix_addtorque.cpp
@@ -48,6 +48,8 @@ FixAddTorque::FixAddTorque(LAMMPS *lmp, int narg, char **arg) :
global_freq = 1;
extscalar = 1;
extvector = 1;
+ respa_level_support = 1;
+ ilevel_respa = 0;
xstr = ystr = zstr = NULL;
@@ -133,8 +135,10 @@ void FixAddTorque::init()
varflag = EQUAL;
else varflag = CONSTANT;
- if (strcmp(update->integrate_style,"respa") == 0)
- nlevels_respa = ((Respa *) update->integrate)->nlevels;
+ if (strstr(update->integrate_style,"respa")) {
+ ilevel_respa = ((Respa *) update->integrate)->nlevels-1;
+ if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa);
+ }
}
/* ---------------------------------------------------------------------- */
@@ -144,9 +148,9 @@ void FixAddTorque::setup(int vflag)
if (strcmp(update->integrate_style,"verlet") == 0)
post_force(vflag);
else {
- ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
- post_force_respa(vflag,nlevels_respa-1,0);
- ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
+ ((Respa *) update->integrate)->copy_flevel_f(ilevel_respa);
+ post_force_respa(vflag,ilevel_respa,0);
+ ((Respa *) update->integrate)->copy_f_flevel(ilevel_respa);
}
}
@@ -249,7 +253,7 @@ void FixAddTorque::post_force(int vflag)
void FixAddTorque::post_force_respa(int vflag, int ilevel, int iloop)
{
- if (ilevel == nlevels_respa-1) post_force(vflag);
+ if (ilevel == ilevel_respa) post_force(vflag);
}
/* ---------------------------------------------------------------------- */
diff --git a/src/USER-MISC/fix_addtorque.h b/src/USER-MISC/fix_addtorque.h
index 513b658c87ac4e6d3b96b6a48c74063cd85ead84..49631ec9163ec32239b12c0fe341e83418c9430a 100644
--- a/src/USER-MISC/fix_addtorque.h
+++ b/src/USER-MISC/fix_addtorque.h
@@ -45,7 +45,7 @@ class FixAddTorque : public Fix {
int xvar,yvar,zvar,xstyle,ystyle,zstyle;
double foriginal[4],foriginal_all[4];
int force_flag;
- int nlevels_respa;
+ int ilevel_respa;
};
}
diff --git a/src/USER-MISC/fix_smd.cpp b/src/USER-MISC/fix_smd.cpp
index 3f6ca3f58a053d1c6830231459ba7f1b7c57bb1c..b527e4b788dacb60f2566bf725e789a9b1fd005c 100644
--- a/src/USER-MISC/fix_smd.cpp
+++ b/src/USER-MISC/fix_smd.cpp
@@ -56,6 +56,8 @@ FixSMD::FixSMD(LAMMPS *lmp, int narg, char **arg) :
size_vector = 7;
global_freq = 1;
extvector = 1;
+ respa_level_support = 1;
+ ilevel_respa = 0;
int argoffs=3;
if (strcmp(arg[argoffs],"cvel") == 0) {
@@ -156,8 +158,10 @@ void FixSMD::init()
zn = dz/r_old;
}
- if (strstr(update->integrate_style,"respa"))
- nlevels_respa = ((Respa *) update->integrate)->nlevels;
+ if (strstr(update->integrate_style,"respa")) {
+ ilevel_respa = ((Respa *) update->integrate)->nlevels-1;
+ if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa);
+ }
}
/* ---------------------------------------------------------------------- */
@@ -167,9 +171,9 @@ void FixSMD::setup(int vflag)
if (strstr(update->integrate_style,"verlet"))
post_force(vflag);
else {
- ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
- post_force_respa(vflag,nlevels_respa-1,0);
- ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
+ ((Respa *) update->integrate)->copy_flevel_f(ilevel_respa);
+ post_force_respa(vflag,ilevel_respa,0);
+ ((Respa *) update->integrate)->copy_f_flevel(ilevel_respa);
}
}
@@ -180,7 +184,12 @@ void FixSMD::post_force(int vflag)
if (styleflag & SMD_TETHER) smd_tether();
else smd_couple();
- if (styleflag & SMD_CVEL) r_old += v_smd * update->dt;
+ if (styleflag & SMD_CVEL) {
+ if (strstr(update->integrate_style,"verlet"))
+ r_old += v_smd * update->dt;
+ else
+ r_old += v_smd * ((Respa *) update->integrate)->step[ilevel_respa];
+ }
}
/* ---------------------------------------------------------------------- */
@@ -190,6 +199,10 @@ void FixSMD::smd_tether()
double xcm[3];
group->xcm(igroup,masstotal,xcm);
+ double dt = update->dt;
+ if (strstr(update->integrate_style,"respa"))
+ dt = ((Respa *) update->integrate)->step[ilevel_respa];
+
// fx,fy,fz = components of k * (r-r0)
double dx,dy,dz,fx,fy,fz,r,dr;
@@ -209,7 +222,7 @@ void FixSMD::smd_tether()
fx = k_smd*dx*dr/r;
fy = k_smd*dy*dr/r;
fz = k_smd*dz*dr/r;
- pmf += (fx*xn + fy*yn + fz*zn) * v_smd * update->dt;
+ pmf += (fx*xn + fy*yn + fz*zn) * v_smd * dt;
} else {
fx = 0.0;
fy = 0.0;
@@ -269,6 +282,10 @@ void FixSMD::smd_couple()
group->xcm(igroup,masstotal,xcm);
group->xcm(igroup2,masstotal2,xcm2);
+ double dt = update->dt;
+ if (strstr(update->integrate_style,"respa"))
+ dt = ((Respa *) update->integrate)->step[ilevel_respa];
+
// renormalize direction of spring
double dx,dy,dz,r,dr;
if (styleflag & SMD_AUTOX) dx = xcm2[0] - xcm[0];
@@ -309,7 +326,7 @@ void FixSMD::smd_couple()
fx = k_smd*dx*dr/r;
fy = k_smd*dy*dr/r;
fz = k_smd*dz*dr/r;
- pmf += (fx*xn + fy*yn + fz*zn) * fsign * v_smd * update->dt;
+ pmf += (fx*xn + fy*yn + fz*zn) * fsign * v_smd * dt;
} else {
fx = 0.0;
fy = 0.0;
@@ -417,7 +434,7 @@ void FixSMD::restart(char *buf)
void FixSMD::post_force_respa(int vflag, int ilevel, int iloop)
{
- if (ilevel == nlevels_respa-1) post_force(vflag);
+ if (ilevel == ilevel_respa) post_force(vflag);
}
/* ----------------------------------------------------------------------
diff --git a/src/USER-MISC/fix_smd.h b/src/USER-MISC/fix_smd.h
index 687c1900b2dbc62346b98a055a8028618b9dc3a6..4f86456f994f5b3e70f9e7ac53777bc67a5693f1 100644
--- a/src/USER-MISC/fix_smd.h
+++ b/src/USER-MISC/fix_smd.h
@@ -46,7 +46,7 @@ class FixSMD : public Fix {
int igroup2,group2bit;
double masstotal,masstotal2;
- int nlevels_respa;
+ int ilevel_respa;
double ftotal[3],ftotal_all[7];
int force_flag;
diff --git a/src/USER-MISC/fix_ttm_mod.cpp b/src/USER-MISC/fix_ttm_mod.cpp
index 5a58dd8489efe1ce48cc1df74207a75a458c5b86..f7b81c72061e9e75a55b761b17754efef8f20cf5 100644
--- a/src/USER-MISC/fix_ttm_mod.cpp
+++ b/src/USER-MISC/fix_ttm_mod.cpp
@@ -78,24 +78,33 @@ FixTTMMod::FixTTMMod(LAMMPS *lmp, int narg, char **arg) :
nevery = 1;
restart_peratom = 1;
restart_global = 1;
- seed = atoi(arg[3]);
- fpr_2 = fopen(arg[4],"r");
- nxnodes = atoi(arg[5]);
- nynodes = atoi(arg[6]);
- nznodes = atoi(arg[7]);
- fpr = fopen(arg[8],"r");
- if (fpr == NULL) {
+
+ seed = force->inumeric(FLERR,arg[3]);
+ if (seed <= 0)
+ error->all(FLERR,"Invalid random number seed in fix ttm/mod command");
+
+ FILE *fpr_2 = force->open_potential(arg[4]);
+ if (fpr_2 == NULL) {
char str[128];
- sprintf(str,"Cannot open file %s",arg[7]);
- error->one(FLERR,str);
+ sprintf(str,"Cannot open file %s",arg[4]);
+ error->all(FLERR,str);
}
+
+ nxnodes = force->inumeric(FLERR,arg[5]);
+ nynodes = force->inumeric(FLERR,arg[6]);
+ nznodes = force->inumeric(FLERR,arg[7]);
+ if (nxnodes <= 0 || nynodes <= 0 || nznodes <= 0)
+ error->all(FLERR,"Fix ttm/mod number of nodes must be > 0");
+
+ FILE *fpr = force->open_potential(arg[8]);
if (fpr == NULL) {
char str[128];
sprintf(str,"Cannot open file %s",arg[8]);
- error->one(FLERR,str);
+ error->all(FLERR,str);
}
- nfileevery = atoi(arg[9]);
- if (nfileevery) {
+
+ nfileevery = force->inumeric(FLERR,arg[9]);
+ if (nfileevery > 0) {
if (narg != 11) error->all(FLERR,"Illegal fix ttm/mod command");
MPI_Comm_rank(world,&me);
if (me == 0) {
@@ -226,19 +235,15 @@ FixTTMMod::FixTTMMod(LAMMPS *lmp, int narg, char **arg) :
mult_factor = intensity;
duration = 0.0;
v_0_sq = v_0*v_0;
- // error checks
- if (nxnodes <= 0 || nynodes <= 0 || nznodes <= 0)
- error->all(FLERR,"Fix ttm number of nodes must be > 0");
surface_double = double(t_surface_l)*(domain->xprd/nxnodes);
if ((C_limit+esheat_0) < 0.0)
- error->all(FLERR,"Fix ttm electronic_specific_heat must be >= 0.0");
+ error->all(FLERR,"Fix ttm/mod electronic_specific_heat must be >= 0.0");
if (electronic_density <= 0.0)
- error->all(FLERR,"Fix ttm electronic_density must be > 0.0");
- if (gamma_p < 0.0) error->all(FLERR,"Fix ttm gamma_p must be >= 0.0");
- if (gamma_s < 0.0) error->all(FLERR,"Fix ttm gamma_s must be >= 0.0");
- if (v_0 < 0.0) error->all(FLERR,"Fix ttm v_0 must be >= 0.0");
- if (ionic_density <= 0.0) error->all(FLERR,"Fix ttm ionic_density must be > 0.0");
- if (seed <= 0) error->all(FLERR,"Invalid random number seed in fix ttm command");
+ error->all(FLERR,"Fix ttm/mod electronic_density must be > 0.0");
+ if (gamma_p < 0.0) error->all(FLERR,"Fix ttm/mod gamma_p must be >= 0.0");
+ if (gamma_s < 0.0) error->all(FLERR,"Fix ttm/mod gamma_s must be >= 0.0");
+ if (v_0 < 0.0) error->all(FLERR,"Fix ttm/mod v_0 must be >= 0.0");
+ if (ionic_density <= 0.0) error->all(FLERR,"Fix ttm/mod ionic_density must be > 0.0");
if (surface_l < 0) error->all(FLERR,"Surface coordinates must be >= 0");
if (surface_l >= surface_r) error->all(FLERR, "Left surface coordinate must be less than right surface coordinate");
// initialize Marsaglia RNG with processor-unique seed
@@ -274,8 +279,10 @@ FixTTMMod::FixTTMMod(LAMMPS *lmp, int narg, char **arg) :
atom->add_callback(0);
atom->add_callback(1);
// set initial electron temperatures from user input file
- if (me == 0) read_initial_electron_temperatures();
+ if (me == 0) read_initial_electron_temperatures(fpr);
MPI_Bcast(&T_electron[0][0][0],total_nnodes,MPI_DOUBLE,0,world);
+ fclose(fpr);
+ fclose(fpr_2);
}
/* ---------------------------------------------------------------------- */
@@ -317,11 +324,11 @@ int FixTTMMod::setmask()
void FixTTMMod::init()
{
if (domain->dimension == 2)
- error->all(FLERR,"Cannot use fix ttm with 2d simulation");
+ error->all(FLERR,"Cannot use fix ttm/mod with 2d simulation");
if (domain->nonperiodic != 0)
- error->all(FLERR,"Cannot use nonperiodic boundares with fix ttm");
+ error->all(FLERR,"Cannot use nonperiodic boundares with fix ttm/mod");
if (domain->triclinic)
- error->all(FLERR,"Cannot use fix ttm with triclinic box");
+ error->all(FLERR,"Cannot use fix ttm/mod with triclinic box");
// set force prefactors
for (int i = 1; i <= atom->ntypes; i++) {
gfactor1[i] = - gamma_p / force->ftm2v;
@@ -488,7 +495,7 @@ void FixTTMMod::reset_dt()
only called by proc 0
------------------------------------------------------------------------- */
-void FixTTMMod::read_initial_electron_temperatures()
+void FixTTMMod::read_initial_electron_temperatures(FILE *fpr)
{
char line[MAXLINE];
for (int ixnode = 0; ixnode < nxnodes; ixnode++)
@@ -501,7 +508,7 @@ void FixTTMMod::read_initial_electron_temperatures()
while (1) {
if (fgets(line,MAXLINE,fpr) == NULL) break;
sscanf(line,"%d %d %d %lg",&ixnode,&iynode,&iznode,&T_tmp);
- if (T_tmp < 0.0) error->one(FLERR,"Fix ttm electron temperatures must be >= 0.0");
+ if (T_tmp < 0.0) error->one(FLERR,"Fix ttm/mod electron temperatures must be >= 0.0");
T_electron[ixnode][iynode][iznode] = T_tmp;
T_initial_set[ixnode][iynode][iznode] = 1;
}
@@ -509,9 +516,7 @@ void FixTTMMod::read_initial_electron_temperatures()
for (int iynode = 0; iynode < nynodes; iynode++)
for (int iznode = 0; iznode < nznodes; iznode++)
if (T_initial_set[ixnode][iynode][iznode] == 0)
- error->one(FLERR,"Initial temperatures not all set in fix ttm");
- // close file
- fclose(fpr);
+ error->one(FLERR,"Initial temperatures not all set in fix ttm/mod");
}
/* ---------------------------------------------------------------------- */
@@ -634,7 +639,7 @@ void FixTTMMod::end_of_step()
num_inner_timesteps = static_cast<unsigned int>(update->dt/inner_dt) + 1;
inner_dt = update->dt/double(num_inner_timesteps);
if (num_inner_timesteps > 1000000)
- error->warning(FLERR,"Too many inner timesteps in fix ttm",0);
+ error->warning(FLERR,"Too many inner timesteps in fix ttm/mod",0);
for (int ith_inner_timestep = 0; ith_inner_timestep < num_inner_timesteps;
ith_inner_timestep++) {
for (int ixnode = 0; ixnode < nxnodes; ixnode++)
diff --git a/src/USER-MISC/fix_ttm_mod.h b/src/USER-MISC/fix_ttm_mod.h
index 11a68e4485869850b3fa226fbb3629d652cadc89..3415baf5b10f13f8cfc1cbffe13cb68cfd0da5f4 100644
--- a/src/USER-MISC/fix_ttm_mod.h
+++ b/src/USER-MISC/fix_ttm_mod.h
@@ -58,7 +58,7 @@ class FixTTMMod : public Fix {
int nlevels_respa;
int seed;
class RanMars *random;
- FILE *fp,*fpr,*fpr_2;
+ FILE *fp;
int nxnodes,nynodes,nznodes,total_nnodes;
int ***nsum;
int ***nsum_all,***T_initial_set;
@@ -79,7 +79,7 @@ class FixTTMMod : public Fix {
double electron_temperature_min;
el_heat_capacity_thermal_conductivity el_properties(double);
double el_sp_heat_integral(double);
- void read_initial_electron_temperatures();
+ void read_initial_electron_temperatures(FILE *);
};
}
diff --git a/src/USER-OMP/fix_gravity_omp.cpp b/src/USER-OMP/fix_gravity_omp.cpp
index b853d844f0994c25c1969febe515478edd6e7b4b..1e65d1f7641abd558040f8cc27f74788ea2f435f 100644
--- a/src/USER-OMP/fix_gravity_omp.cpp
+++ b/src/USER-OMP/fix_gravity_omp.cpp
@@ -108,5 +108,5 @@ void FixGravityOMP::post_force(int vflag)
void FixGravityOMP::post_force_respa(int vflag, int ilevel, int iloop)
{
- if (ilevel == nlevels_respa-1) post_force(vflag);
+ if (ilevel == ilevel_respa) post_force(vflag);
}
diff --git a/src/USER-OMP/fix_qeq_comb_omp.cpp b/src/USER-OMP/fix_qeq_comb_omp.cpp
index f8de18498da9649be8c712c1bdc6edee93c5053b..0c69876b86fd3af322bee1f5cca6d9aee8e869ec 100644
--- a/src/USER-OMP/fix_qeq_comb_omp.cpp
+++ b/src/USER-OMP/fix_qeq_comb_omp.cpp
@@ -63,8 +63,10 @@ void FixQEQCombOMP::init()
error->all(FLERR,"Must use pair_style comb or "
"comb/omp with fix qeq/comb/omp");
- if (strstr(update->integrate_style,"respa"))
- nlevels_respa = ((Respa *) update->integrate)->nlevels;
+ if (strstr(update->integrate_style,"respa")) {
+ ilevel_respa = ((Respa *) update->integrate)->nlevels-1;
+ if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa);
+ }
ngroup = group->count(igroup);
if (ngroup == 0) error->all(FLERR,"Fix qeq/comb group has no atoms");
diff --git a/src/USER-OMP/respa_omp.cpp b/src/USER-OMP/respa_omp.cpp
index d925088edaa9afc1200c3c7e2764d28cc49188c9..738538a20935ffee2d8077a7713c6bab34b5f0c5 100644
--- a/src/USER-OMP/respa_omp.cpp
+++ b/src/USER-OMP/respa_omp.cpp
@@ -73,7 +73,18 @@ void RespaOMP::setup()
fprintf(screen,"Setting up r-RESPA/omp run ...\n");
fprintf(screen," Unit style : %s\n", update->unit_style);
fprintf(screen," Current step : " BIGINT_FORMAT "\n", update->ntimestep);
- fprintf(screen," OuterTime step: %g\n", update->dt);
+ fprintf(screen," Time steps :");
+ for (int ilevel=0; ilevel < nlevels; ++ilevel)
+ fprintf(screen," %d:%g",ilevel+1, step[ilevel]);
+ fprintf(screen,"\n r-RESPA fixes :");
+ for (int l=0; l < modify->n_post_force_respa; ++l) {
+ Fix *f = modify->fix[modify->list_post_force_respa[l]];
+ if (f->respa_level >= 0)
+ fprintf(screen," %d:%s[%s]",
+ MIN(f->respa_level+1,nlevels),f->style,f->id);
+ }
+ fprintf(screen,"\n");
+ timer->print_timeout(screen);
}
update->setupflag = 1;
diff --git a/src/USER-SMD/fix_smd_setvel.cpp b/src/USER-SMD/fix_smd_setvel.cpp
index 568ffb79775fffaf66f17f809edf27b05dd1bd97..92c4a005818c8da49d7c5446bb271acbd52b9c48 100644
--- a/src/USER-SMD/fix_smd_setvel.cpp
+++ b/src/USER-SMD/fix_smd_setvel.cpp
@@ -31,7 +31,6 @@
#include "modify.h"
#include "domain.h"
#include "region.h"
-#include "respa.h"
#include "input.h"
#include "variable.h"
#include "memory.h"
@@ -192,9 +191,6 @@ void FixSMDSetVel::init() {
else
varflag = CONSTANT;
- if (strstr(update->integrate_style, "respa"))
- nlevels_respa = ((Respa *) update->integrate)->nlevels;
-
// cannot use non-zero forces for a minimization since no energy is integrated
// use fix addforce instead
@@ -223,11 +219,7 @@ void FixSMDSetVel::setup(int vflag) {
if (strstr(update->integrate_style, "verlet"))
post_force(vflag);
else
- for (int ilevel = 0; ilevel < nlevels_respa; ilevel++) {
- ((Respa *) update->integrate)->copy_flevel_f(ilevel);
- post_force_respa(vflag, ilevel, 0);
- ((Respa *) update->integrate)->copy_f_flevel(ilevel);
- }
+ error->all(FLERR,"Fix smd/setvel does not support RESPA");
}
/* ---------------------------------------------------------------------- */
diff --git a/src/atom.cpp b/src/atom.cpp
index ac428b9a95c74ba0c65104cf2267a3dbcbdbde02..acb3cbce7cd0eac75dd40b46b3c368b70da0b212 100644
--- a/src/atom.cpp
+++ b/src/atom.cpp
@@ -94,8 +94,9 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp)
// USER-DPD
uCond = uMech = uChem = uCG = uCGnew = NULL;
- duCond = duMech = duChem = NULL;
+ duChem = NULL;
dpdTheta = NULL;
+ ssaAIR = NULL;
// USER-SMD
@@ -284,9 +285,8 @@ Atom::~Atom()
memory->destroy(uChem);
memory->destroy(uCG);
memory->destroy(uCGnew);
- memory->destroy(duCond);
- memory->destroy(duMech);
memory->destroy(duChem);
+ memory->destroy(ssaAIR);
memory->destroy(nspecial);
memory->destroy(special);
diff --git a/src/atom.h b/src/atom.h
index 4168df32cfcf6e98223525a5554c610c05d75158..5e6f257789e089aae5b13d6758d1eeb0409888d2 100644
--- a/src/atom.h
+++ b/src/atom.h
@@ -88,9 +88,10 @@ class Atom : protected Pointers {
// USER-DPD package
double *uCond,*uMech,*uChem,*uCGnew,*uCG;
- double *duCond,*duMech,*duChem;
+ double *duChem;
double *dpdTheta;
int nspecies_dpd;
+ int *ssaAIR; // Shardlow Splitting Algorithm Active Interaction Region number
// molecular info
diff --git a/src/balance.cpp b/src/balance.cpp
index c615db61129b35f223d47c541f10a2bee8a50cf6..db7361703c4cdb426dcd95ece9cd70707bf746db 100644
--- a/src/balance.cpp
+++ b/src/balance.cpp
@@ -325,10 +325,6 @@ void Balance::command(int narg, char **arg)
bisection(1);
}
- // output of final result
-
- if (outflag) dumpout(update->ntimestep,fp);
-
// reset proc sub-domains
// for either brick or tiled comm style
@@ -344,6 +340,10 @@ void Balance::command(int narg, char **arg)
delete irregular;
if (domain->triclinic) domain->lamda2x(atom->nlocal);
+ // output of final result
+
+ if (outflag) dumpout(update->ntimestep,fp);
+
// check if any atoms were lost
bigint natoms;
diff --git a/src/dump_image.cpp b/src/dump_image.cpp
index b9ff534c6616023dbdb8780f8486b10fcfc241c0..5e6701ff344d6e9d3f457ea6766b2ca14c0d646a 100644
--- a/src/dump_image.cpp
+++ b/src/dump_image.cpp
@@ -893,7 +893,6 @@ void DumpImage::create_image()
int *molatom = atom->molatom;
int *type = atom->type;
int nlocal = atom->nlocal;
- int nall = atom->nlocal + atom->nghost;
int newton_bond = force->newton_bond;
int molecular = atom->molecular;
Molecule **onemols = atom->avec->onemols;
diff --git a/src/fix.cpp b/src/fix.cpp
index 165536d8342c172a215e64b4993425bc77159b4f..1a3b6a42c98052f1b6051024216984e66a772436 100644
--- a/src/fix.cpp
+++ b/src/fix.cpp
@@ -16,6 +16,7 @@
#include "fix.h"
#include "atom.h"
#include "group.h"
+#include "force.h"
#include "atom_masks.h"
#include "memory.h"
#include "error.h"
@@ -68,6 +69,8 @@ Fix::Fix(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
dof_flag = 0;
special_alter_flag = 0;
enforce2d_flag = 0;
+ respa_level_support = 0;
+ respa_level = -1;
scalar_flag = vector_flag = array_flag = 0;
peratom_flag = local_flag = 0;
@@ -130,6 +133,13 @@ void Fix::modify_params(int narg, char **arg)
else if (strcmp(arg[iarg+1],"yes") == 0) thermo_energy = 1;
else error->all(FLERR,"Illegal fix_modify command");
iarg += 2;
+ } else if (strcmp(arg[iarg],"respa") == 0) {
+ if (iarg+2 > narg) error->all(FLERR,"Illegal fix_modify command");
+ if (!respa_level_support) error->all(FLERR,"Illegal fix_modify command");
+ int lvl = force->inumeric(FLERR,arg[1]);
+ if (lvl < 0) error->all(FLERR,"Illegal fix_modify command");
+ respa_level = lvl-1;
+ iarg += 2;
} else {
int n = modify_param(narg-iarg,&arg[iarg]);
if (n == 0) error->all(FLERR,"Illegal fix_modify command");
diff --git a/src/fix.h b/src/fix.h
index 4bf2609effa021c008c16b153db1a98e3be51103..a2aa3782cc658121d43342a111c9c7dddf09f0c4 100644
--- a/src/fix.h
+++ b/src/fix.h
@@ -52,6 +52,8 @@ class Fix : protected Pointers {
int dof_flag; // 1 if has dof() method (not min_dof())
int special_alter_flag; // 1 if has special_alter() meth for spec lists
int enforce2d_flag; // 1 if has enforce2d method
+ int respa_level_support; // 1 if fix supports fix_modify respa
+ int respa_level; // which respa level to apply fix (1-Nrespa)
int scalar_flag; // 0/1 if compute_scalar() function exists
int vector_flag; // 0/1 if compute_vector() function exists
diff --git a/src/fix_addforce.cpp b/src/fix_addforce.cpp
index 2b5c682266aba970424e5341027bfd9ab7c0b5f5..ac8384f2c0ac8a51c164f0870cbfd12083ae75e0 100644
--- a/src/fix_addforce.cpp
+++ b/src/fix_addforce.cpp
@@ -47,6 +47,8 @@ FixAddForce::FixAddForce(LAMMPS *lmp, int narg, char **arg) :
global_freq = 1;
extscalar = 1;
extvector = 1;
+ respa_level_support = 1;
+ ilevel_respa = 0;
xstr = ystr = zstr = NULL;
@@ -201,8 +203,13 @@ void FixAddForce::init()
update->whichflag == 2 && estyle == NONE)
error->all(FLERR,"Must use variable energy with fix addforce");
+ int max_respa = 0;
+
if (strstr(update->integrate_style,"respa"))
- nlevels_respa = ((Respa *) update->integrate)->nlevels;
+ max_respa = ((Respa *) update->integrate)->nlevels-1;
+
+ if (respa_level >= 0)
+ ilevel_respa = MIN(respa_level,max_respa);
}
/* ---------------------------------------------------------------------- */
@@ -212,9 +219,9 @@ void FixAddForce::setup(int vflag)
if (strstr(update->integrate_style,"verlet"))
post_force(vflag);
else {
- ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
- post_force_respa(vflag,nlevels_respa-1,0);
- ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
+ ((Respa *) update->integrate)->copy_flevel_f(ilevel_respa);
+ post_force_respa(vflag,ilevel_respa,0);
+ ((Respa *) update->integrate)->copy_f_flevel(ilevel_respa);
}
}
@@ -324,7 +331,7 @@ void FixAddForce::post_force(int vflag)
void FixAddForce::post_force_respa(int vflag, int ilevel, int iloop)
{
- if (ilevel == nlevels_respa-1) post_force(vflag);
+ if (ilevel == ilevel_respa) post_force(vflag);
}
/* ---------------------------------------------------------------------- */
diff --git a/src/fix_addforce.h b/src/fix_addforce.h
index d590463a5efa52fa2d6f168b73525dd6724a9e5d..d19c0db0d1d170206a0e6873f4e8807a69d00a3e 100644
--- a/src/fix_addforce.h
+++ b/src/fix_addforce.h
@@ -47,7 +47,7 @@ class FixAddForce : public Fix {
int xvar,yvar,zvar,evar,xstyle,ystyle,zstyle,estyle;
double foriginal[4],foriginal_all[4];
int force_flag;
- int nlevels_respa;
+ int ilevel_respa;
int maxatom;
double **sforce;
diff --git a/src/fix_aveforce.cpp b/src/fix_aveforce.cpp
index 443ae6ea57110d47a7425d703ce2670524d37813..09a685cb1a150d11d0c0c880beb7648ee404b071 100644
--- a/src/fix_aveforce.cpp
+++ b/src/fix_aveforce.cpp
@@ -43,6 +43,8 @@ FixAveForce::FixAveForce(LAMMPS *lmp, int narg, char **arg) :
size_vector = 3;
global_freq = 1;
extvector = 1;
+ respa_level_support = 1;
+ ilevel_respa = nlevels_respa = 0;
xstr = ystr = zstr = NULL;
@@ -161,8 +163,11 @@ void FixAveForce::init()
if (xstyle == EQUAL || ystyle == EQUAL || zstyle == EQUAL) varflag = EQUAL;
else varflag = CONSTANT;
- if (strstr(update->integrate_style,"respa"))
+ if (strstr(update->integrate_style,"respa")) {
nlevels_respa = ((Respa *) update->integrate)->nlevels;
+ if (respa_level >= 0) ilevel_respa = MIN(respa_level,nlevels_respa-1);
+ else ilevel_respa = nlevels_respa-1;
+ }
}
/* ---------------------------------------------------------------------- */
@@ -255,10 +260,10 @@ void FixAveForce::post_force(int vflag)
void FixAveForce::post_force_respa(int vflag, int ilevel, int iloop)
{
- // ave + extra force on outermost level
- // just ave on inner levels
+ // ave + extra force on selected RESPA level
+ // just ave on all other levels
- if (ilevel == nlevels_respa-1) post_force(vflag);
+ if (ilevel == ilevel_respa) post_force(vflag);
else {
Region *region = NULL;
if (iregion >= 0) {
diff --git a/src/fix_aveforce.h b/src/fix_aveforce.h
index 2ec311034367a8eef63919f4ea4edaa3d5da98be..a6228be7dfb473303ed57f9fd0bf8236e90ec3c4 100644
--- a/src/fix_aveforce.h
+++ b/src/fix_aveforce.h
@@ -45,7 +45,7 @@ class FixAveForce : public Fix {
int xvar,yvar,zvar,xstyle,ystyle,zstyle;
int iregion;
double foriginal_all[4];
- int nlevels_respa;
+ int nlevels_respa,ilevel_respa;
};
}
diff --git a/src/fix_drag.cpp b/src/fix_drag.cpp
index 6a28e89806955e7d16df3da567eb8c223e7dbe04..3fda219126cbccc3bc8463bb07eea26ae96d3259 100644
--- a/src/fix_drag.cpp
+++ b/src/fix_drag.cpp
@@ -36,6 +36,8 @@ FixDrag::FixDrag(LAMMPS *lmp, int narg, char **arg) :
size_vector = 3;
global_freq = 1;
extvector = 1;
+ respa_level_support = 1;
+ ilevel_respa = 0;
xflag = yflag = zflag = 1;
@@ -67,8 +69,10 @@ int FixDrag::setmask()
void FixDrag::init()
{
- if (strstr(update->integrate_style,"respa"))
- nlevels_respa = ((Respa *) update->integrate)->nlevels;
+ if (strstr(update->integrate_style,"respa")) {
+ ilevel_respa = ((Respa *) update->integrate)->nlevels-1;
+ if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa);
+ }
}
/* ---------------------------------------------------------------------- */
@@ -78,9 +82,9 @@ void FixDrag::setup(int vflag)
if (strstr(update->integrate_style,"verlet"))
post_force(vflag);
else {
- ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
- post_force_respa(vflag,nlevels_respa-1,0);
- ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
+ ((Respa *) update->integrate)->copy_flevel_f(ilevel_respa);
+ post_force_respa(vflag,ilevel_respa,0);
+ ((Respa *) update->integrate)->copy_f_flevel(ilevel_respa);
}
}
@@ -130,7 +134,7 @@ void FixDrag::post_force(int vflag)
void FixDrag::post_force_respa(int vflag, int ilevel, int iloop)
{
- if (ilevel == nlevels_respa-1) post_force(vflag);
+ if (ilevel == ilevel_respa) post_force(vflag);
}
/* ----------------------------------------------------------------------
diff --git a/src/fix_drag.h b/src/fix_drag.h
index dc0075143d06483905e649daecb9f62d1f91d8e7..fa444a4b35dd5794f2d59321e4c514594459273d 100644
--- a/src/fix_drag.h
+++ b/src/fix_drag.h
@@ -39,7 +39,7 @@ class FixDrag : public Fix {
double f_mag;
int xflag,yflag,zflag;
double delta;
- int nlevels_respa;
+ int ilevel_respa;
double ftotal[3],ftotal_all[3];
int force_flag;
};
diff --git a/src/fix_gravity.cpp b/src/fix_gravity.cpp
index a08213e3d027dc7c7d0905195239d2f618090297..95f899ae86f08cff92de84c33385607dc1922e6f 100644
--- a/src/fix_gravity.cpp
+++ b/src/fix_gravity.cpp
@@ -45,6 +45,8 @@ FixGravity::FixGravity(LAMMPS *lmp, int narg, char **arg) :
scalar_flag = 1;
global_freq = 1;
extscalar = 1;
+ respa_level_support = 1;
+ ilevel_respa = 0;
mstr = vstr = pstr = tstr = xstr = ystr = zstr = NULL;
mstyle = vstyle = pstyle = tstyle = xstyle = ystyle = zstyle = CONSTANT;
@@ -162,8 +164,10 @@ int FixGravity::setmask()
void FixGravity::init()
{
- if (strstr(update->integrate_style,"respa"))
- nlevels_respa = ((Respa *) update->integrate)->nlevels;
+ if (strstr(update->integrate_style,"respa")) {
+ ilevel_respa = ((Respa *) update->integrate)->nlevels-1;
+ if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa);
+ }
// check variables
@@ -234,9 +238,9 @@ void FixGravity::setup(int vflag)
if (strstr(update->integrate_style,"verlet"))
post_force(vflag);
else {
- ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
- post_force_respa(vflag,nlevels_respa-1,0);
- ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
+ ((Respa *) update->integrate)->copy_flevel_f(ilevel_respa);
+ post_force_respa(vflag,ilevel_respa,0);
+ ((Respa *) update->integrate)->copy_f_flevel(ilevel_respa);
}
}
@@ -297,7 +301,7 @@ void FixGravity::post_force(int vflag)
void FixGravity::post_force_respa(int vflag, int ilevel, int iloop)
{
- if (ilevel == nlevels_respa-1) post_force(vflag);
+ if (ilevel == ilevel_respa) post_force(vflag);
}
/* ---------------------------------------------------------------------- */
diff --git a/src/fix_gravity.h b/src/fix_gravity.h
index f37a5bb411168615fdf91cea45df2b0a9903ead0..35be23a40a12e37d35224f792426f34c6eafc538 100644
--- a/src/fix_gravity.h
+++ b/src/fix_gravity.h
@@ -44,7 +44,7 @@ class FixGravity : public Fix {
double xdir,ydir,zdir;
double xgrav,ygrav,zgrav,xacc,yacc,zacc;
double degree2rad;
- int nlevels_respa;
+ int ilevel_respa;
int time_origin;
int eflag;
double egrav,egrav_all;
diff --git a/src/fix_indent.cpp b/src/fix_indent.cpp
index 4b85dd2788fd302e35ffa39bc364ad2788035226..a6693e98a23e61e0cbe29a713586134ab1399b86 100644
--- a/src/fix_indent.cpp
+++ b/src/fix_indent.cpp
@@ -50,6 +50,8 @@ FixIndent::FixIndent(LAMMPS *lmp, int narg, char **arg) :
global_freq = 1;
extscalar = 1;
extvector = 1;
+ respa_level_support = 1;
+ ilevel_respa = 0;
k = force->numeric(FLERR,arg[3]);
k3 = k/3.0;
@@ -151,8 +153,10 @@ void FixIndent::init()
error->all(FLERR,"Variable for fix indent is not equal style");
}
- if (strstr(update->integrate_style,"respa"))
- nlevels_respa = ((Respa *) update->integrate)->nlevels;
+ if (strstr(update->integrate_style,"respa")) {
+ ilevel_respa = ((Respa *) update->integrate)->nlevels-1;
+ if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa);
+ }
}
/* ---------------------------------------------------------------------- */
@@ -162,9 +166,9 @@ void FixIndent::setup(int vflag)
if (strstr(update->integrate_style,"verlet"))
post_force(vflag);
else {
- ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
- post_force_respa(vflag,nlevels_respa-1,0);
- ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
+ ((Respa *) update->integrate)->copy_flevel_f(ilevel_respa);
+ post_force_respa(vflag,ilevel_respa,0);
+ ((Respa *) update->integrate)->copy_f_flevel(ilevel_respa);
}
}
@@ -354,7 +358,7 @@ void FixIndent::post_force(int vflag)
void FixIndent::post_force_respa(int vflag, int ilevel, int iloop)
{
- if (ilevel == nlevels_respa-1) post_force(vflag);
+ if (ilevel == ilevel_respa) post_force(vflag);
}
/* ---------------------------------------------------------------------- */
diff --git a/src/fix_indent.h b/src/fix_indent.h
index 98f8d081ea330ff1005fec279ff25df6ec4039bd..c3a8882734d27ec366a304fb720df065218c8bbc 100644
--- a/src/fix_indent.h
+++ b/src/fix_indent.h
@@ -47,7 +47,7 @@ class FixIndent : public Fix {
int indenter_flag,planeside;
double indenter[4],indenter_all[4];
int cdim,varflag;
- int nlevels_respa;
+ int ilevel_respa;
void options(int, char **);
};
diff --git a/src/fix_restrain.cpp b/src/fix_restrain.cpp
index 97e898494eab6be77b373d9b1433d4576c405687..a92352e56315bb3eef259c11adfd28915e3acc45 100644
--- a/src/fix_restrain.cpp
+++ b/src/fix_restrain.cpp
@@ -51,6 +51,8 @@ FixRestrain::FixRestrain(LAMMPS *lmp, int narg, char **arg) :
scalar_flag = 1;
global_freq = 1;
extscalar = 1;
+ respa_level_support = 1;
+ ilevel_respa = 0;
// parse args
@@ -147,8 +149,10 @@ int FixRestrain::setmask()
void FixRestrain::init()
{
- if (strcmp(update->integrate_style,"respa") == 0)
- nlevels_respa = ((Respa *) update->integrate)->nlevels;
+ if (strstr(update->integrate_style,"respa")) {
+ ilevel_respa = ((Respa *) update->integrate)->nlevels-1;
+ if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa);
+ }
}
/* ---------------------------------------------------------------------- */
@@ -158,9 +162,9 @@ void FixRestrain::setup(int vflag)
if (strcmp(update->integrate_style,"verlet") == 0)
post_force(vflag);
else {
- ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
- post_force_respa(vflag,nlevels_respa-1,0);
- ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
+ ((Respa *) update->integrate)->copy_flevel_f(ilevel_respa);
+ post_force_respa(vflag,ilevel_respa,0);
+ ((Respa *) update->integrate)->copy_f_flevel(ilevel_respa);
}
}
@@ -187,7 +191,7 @@ void FixRestrain::post_force(int vflag)
void FixRestrain::post_force_respa(int vflag, int ilevel, int iloop)
{
- if (ilevel == nlevels_respa-1) post_force(vflag);
+ if (ilevel == ilevel_respa) post_force(vflag);
}
/* ---------------------------------------------------------------------- */
diff --git a/src/fix_restrain.h b/src/fix_restrain.h
index 0a5669a7e83dcffc7051d44654180481b4aa2608..a6f23f86c6c311afeee63dac0b1da5b9456e20ce 100644
--- a/src/fix_restrain.h
+++ b/src/fix_restrain.h
@@ -38,7 +38,7 @@ class FixRestrain : public Fix {
double compute_scalar();
private:
- int nlevels_respa;
+ int ilevel_respa;
int nrestrain,maxrestrain;
int *rstyle;
int **ids;
diff --git a/src/fix_setforce.cpp b/src/fix_setforce.cpp
index d73ab0db63faf5b72169918dcceec6c136eaf02e..08f177ae8c642da35e76b29ebf9e738d9344a562 100644
--- a/src/fix_setforce.cpp
+++ b/src/fix_setforce.cpp
@@ -43,7 +43,8 @@ FixSetForce::FixSetForce(LAMMPS *lmp, int narg, char **arg) :
size_vector = 3;
global_freq = 1;
extvector = 1;
-
+ respa_level_support = 1;
+ ilevel_respa = nlevels_respa = 0;
xstr = ystr = zstr = NULL;
if (strstr(arg[3],"v_") == arg[3]) {
@@ -172,8 +173,11 @@ void FixSetForce::init()
varflag = EQUAL;
else varflag = CONSTANT;
- if (strstr(update->integrate_style,"respa"))
+ if (strstr(update->integrate_style,"respa")) {
nlevels_respa = ((Respa *) update->integrate)->nlevels;
+ if (respa_level >= 0) ilevel_respa = MIN(respa_level,nlevels_respa-1);
+ else ilevel_respa = nlevels_respa-1;
+ }
// cannot use non-zero forces for a minimization since no energy is integrated
// use fix addforce instead
@@ -290,9 +294,9 @@ void FixSetForce::post_force(int vflag)
void FixSetForce::post_force_respa(int vflag, int ilevel, int iloop)
{
- // set force to desired value on outermost level, 0.0 on other levels
+ // set force to desired value on requested level, 0.0 on other levels
- if (ilevel == nlevels_respa-1) post_force(vflag);
+ if (ilevel == ilevel_respa) post_force(vflag);
else {
Region *region = NULL;
if (iregion >= 0) {
diff --git a/src/fix_setforce.h b/src/fix_setforce.h
index a445f342ed896adae45bdc60f002d1ef63fe91fe..54d17049f741b162c518a6b15bddfc94d91b8e7a 100644
--- a/src/fix_setforce.h
+++ b/src/fix_setforce.h
@@ -36,6 +36,7 @@ class FixSetForce : public Fix {
void post_force_respa(int, int, int);
void min_post_force(int);
double compute_vector(int);
+
double memory_usage();
protected:
@@ -46,7 +47,7 @@ class FixSetForce : public Fix {
int xvar,yvar,zvar,xstyle,ystyle,zstyle;
double foriginal[3],foriginal_all[3];
int force_flag;
- int nlevels_respa;
+ int nlevels_respa,ilevel_respa;
int maxatom;
double **sforce;
diff --git a/src/fix_spring.cpp b/src/fix_spring.cpp
index c90feaaafed697237871bf0cc173d166b3867a62..5b96e1e95d482d49a4f9dcbb1a5a20c93f0650b6 100644
--- a/src/fix_spring.cpp
+++ b/src/fix_spring.cpp
@@ -48,6 +48,8 @@ FixSpring::FixSpring(LAMMPS *lmp, int narg, char **arg) :
extscalar = 1;
extvector = 1;
dynamic_group_allow = 1;
+ respa_level_support = 1;
+ ilevel_respa = 0;
group2 = NULL;
@@ -130,8 +132,10 @@ void FixSpring::init()
masstotal = group->mass(igroup);
if (styleflag == COUPLE) masstotal2 = group->mass(igroup2);
- if (strstr(update->integrate_style,"respa"))
- nlevels_respa = ((Respa *) update->integrate)->nlevels;
+ if (strstr(update->integrate_style,"respa")) {
+ ilevel_respa = ((Respa *) update->integrate)->nlevels-1;
+ if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa);
+ }
}
/* ---------------------------------------------------------------------- */
@@ -141,9 +145,9 @@ void FixSpring::setup(int vflag)
if (strstr(update->integrate_style,"verlet"))
post_force(vflag);
else {
- ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
- post_force_respa(vflag,nlevels_respa-1,0);
- ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
+ ((Respa *) update->integrate)->copy_flevel_f(ilevel_respa);
+ post_force_respa(vflag,ilevel_respa,0);
+ ((Respa *) update->integrate)->copy_f_flevel(ilevel_respa);
}
}
@@ -334,7 +338,7 @@ void FixSpring::spring_couple()
void FixSpring::post_force_respa(int vflag, int ilevel, int iloop)
{
- if (ilevel == nlevels_respa-1) post_force(vflag);
+ if (ilevel == ilevel_respa) post_force(vflag);
}
/* ---------------------------------------------------------------------- */
diff --git a/src/fix_spring.h b/src/fix_spring.h
index f784ad9d3450d94567dc5e75f17300282937cd29..1bb551068bd979aade8c0cab0944da62bb350ab5 100644
--- a/src/fix_spring.h
+++ b/src/fix_spring.h
@@ -46,7 +46,7 @@ class FixSpring : public Fix {
char *group2;
int igroup2,group2bit;
double masstotal,masstotal2;
- int nlevels_respa;
+ int ilevel_respa;
double espring,ftotal[4];
void spring_tether();
diff --git a/src/fix_spring_chunk.cpp b/src/fix_spring_chunk.cpp
index d6dc36a73c09c527bfe3b1696ddf94557133a2ad..16b3d7939e47f90780c99f183a2c8189d4e98755 100644
--- a/src/fix_spring_chunk.cpp
+++ b/src/fix_spring_chunk.cpp
@@ -41,6 +41,8 @@ FixSpringChunk::FixSpringChunk(LAMMPS *lmp, int narg, char **arg) :
scalar_flag = 1;
global_freq = 1;
extscalar = 1;
+ respa_level_support = 1;
+ ilevel_respa = 0;
k_spring = force->numeric(FLERR,arg[3]);
@@ -114,8 +116,10 @@ void FixSpringChunk::init()
if (strcmp(idchunk,ccom->idchunk) != 0)
error->all(FLERR,"Fix spring chunk chunkID not same as comID chunkID");
- if (strstr(update->integrate_style,"respa"))
- nlevels_respa = ((Respa *) update->integrate)->nlevels;
+ if (strstr(update->integrate_style,"respa")) {
+ ilevel_respa = ((Respa *) update->integrate)->nlevels-1;
+ if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa);
+ }
}
/* ---------------------------------------------------------------------- */
@@ -125,9 +129,9 @@ void FixSpringChunk::setup(int vflag)
if (strstr(update->integrate_style,"verlet"))
post_force(vflag);
else {
- ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
- post_force_respa(vflag,nlevels_respa-1,0);
- ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
+ ((Respa *) update->integrate)->copy_flevel_f(ilevel_respa);
+ post_force_respa(vflag,ilevel_respa,0);
+ ((Respa *) update->integrate)->copy_f_flevel(ilevel_respa);
}
}
@@ -229,7 +233,7 @@ void FixSpringChunk::post_force(int vflag)
void FixSpringChunk::post_force_respa(int vflag, int ilevel, int iloop)
{
- if (ilevel == nlevels_respa-1) post_force(vflag);
+ if (ilevel == ilevel_respa) post_force(vflag);
}
/* ---------------------------------------------------------------------- */
diff --git a/src/fix_spring_chunk.h b/src/fix_spring_chunk.h
index c95d9063cab1bd11dd3b5dc8e7c80703fa2f968c..1453a8e621834361cd00cf87e5b711576799d50c 100644
--- a/src/fix_spring_chunk.h
+++ b/src/fix_spring_chunk.h
@@ -38,7 +38,7 @@ class FixSpringChunk : public Fix {
double compute_scalar();
private:
- int nlevels_respa;
+ int ilevel_respa;
double k_spring;
double esprings;
char *idchunk,*idcom;
diff --git a/src/fix_spring_rg.cpp b/src/fix_spring_rg.cpp
index 4145e0104771b15cfc8b91396246f0cd59488955..b7fd33ce7749c2dfb1b0c445edf41273a1a18a1c 100644
--- a/src/fix_spring_rg.cpp
+++ b/src/fix_spring_rg.cpp
@@ -44,6 +44,8 @@ FixSpringRG::FixSpringRG(LAMMPS *lmp, int narg, char **arg) :
else rg0 = force->numeric(FLERR,arg[4]);
dynamic_group_allow = 1;
+ respa_level_support = 1;
+ ilevel_respa = 0;
}
/* ---------------------------------------------------------------------- */
@@ -72,8 +74,10 @@ void FixSpringRG::init()
rg0_flag = 0;
}
- if (strstr(update->integrate_style,"respa"))
- nlevels_respa = ((Respa *) update->integrate)->nlevels;
+ if (strstr(update->integrate_style,"respa")) {
+ ilevel_respa = ((Respa *) update->integrate)->nlevels-1;
+ if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa);
+ }
}
/* ---------------------------------------------------------------------- */
@@ -83,9 +87,9 @@ void FixSpringRG::setup(int vflag)
if (strstr(update->integrate_style,"verlet"))
post_force(vflag);
else {
- ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
- post_force_respa(vflag,nlevels_respa-1,0);
- ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
+ ((Respa *) update->integrate)->copy_flevel_f(ilevel_respa);
+ post_force_respa(vflag,ilevel_respa,0);
+ ((Respa *) update->integrate)->copy_f_flevel(ilevel_respa);
}
}
@@ -140,5 +144,5 @@ void FixSpringRG::post_force(int vflag)
void FixSpringRG::post_force_respa(int vflag, int ilevel, int iloop)
{
- if (ilevel == nlevels_respa-1) post_force(vflag);
+ if (ilevel == ilevel_respa) post_force(vflag);
}
diff --git a/src/fix_spring_rg.h b/src/fix_spring_rg.h
index c3a01298432d9ad0bb43ee65e5f89a23cc28390f..17337be6fdf57b421cc0276cc79e3a2abd18c94f 100644
--- a/src/fix_spring_rg.h
+++ b/src/fix_spring_rg.h
@@ -34,7 +34,7 @@ class FixSpringRG : public Fix {
void post_force_respa(int, int, int);
private:
- int nlevels_respa,rg0_flag;
+ int ilevel_respa,rg0_flag;
double rg0,k,masstotal;
};
diff --git a/src/fix_spring_self.cpp b/src/fix_spring_self.cpp
index 830af421621100bc451c926837901cdfcd9cc76c..9517cf0c316bdcb4d0bacb432937116717869dcb 100644
--- a/src/fix_spring_self.cpp
+++ b/src/fix_spring_self.cpp
@@ -41,6 +41,7 @@ FixSpringSelf::FixSpringSelf(LAMMPS *lmp, int narg, char **arg) :
scalar_flag = 1;
global_freq = 1;
extscalar = 1;
+ respa_level_support = 1;
k = force->numeric(FLERR,arg[3]);
if (k <= 0.0) error->all(FLERR,"Illegal fix spring/self command");
@@ -118,8 +119,10 @@ int FixSpringSelf::setmask()
void FixSpringSelf::init()
{
- if (strstr(update->integrate_style,"respa"))
- nlevels_respa = ((Respa *) update->integrate)->nlevels;
+ if (strstr(update->integrate_style,"respa")) {
+ ilevel_respa = ((Respa *) update->integrate)->nlevels-1;
+ if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa);
+ }
}
/* ---------------------------------------------------------------------- */
@@ -129,9 +132,9 @@ void FixSpringSelf::setup(int vflag)
if (strstr(update->integrate_style,"verlet"))
post_force(vflag);
else {
- ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
- post_force_respa(vflag,nlevels_respa-1,0);
- ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
+ ((Respa *) update->integrate)->copy_flevel_f(ilevel_respa);
+ post_force_respa(vflag,ilevel_respa,0);
+ ((Respa *) update->integrate)->copy_f_flevel(ilevel_respa);
}
}
@@ -179,7 +182,7 @@ void FixSpringSelf::post_force(int vflag)
void FixSpringSelf::post_force_respa(int vflag, int ilevel, int iloop)
{
- if (ilevel == nlevels_respa-1) post_force(vflag);
+ if (ilevel == ilevel_respa) post_force(vflag);
}
/* ---------------------------------------------------------------------- */
diff --git a/src/fix_spring_self.h b/src/fix_spring_self.h
index 3eb5945abd3bc015ba35ef6a2f1f6468c4d563fc..0e5ddda127ca64af60c1a838af1eca420c904c5a 100644
--- a/src/fix_spring_self.h
+++ b/src/fix_spring_self.h
@@ -51,7 +51,7 @@ class FixSpringSelf : public Fix {
double k,espring;
double **xoriginal; // original coords of atoms
int xflag, yflag, zflag;
- int nlevels_respa;
+ int ilevel_respa;
};
}
diff --git a/src/fix_viscous.cpp b/src/fix_viscous.cpp
index bc5059ae69aed2397416b87794ca1e7e6814d5fd..9c82b34d8d95d1b961ae57ffaf1d8c227b54911c 100644
--- a/src/fix_viscous.cpp
+++ b/src/fix_viscous.cpp
@@ -49,6 +49,9 @@ FixViscous::FixViscous(LAMMPS *lmp, int narg, char **arg) :
iarg += 3;
} else error->all(FLERR,"Illegal fix viscous command");
}
+
+ respa_level_support = 1;
+ ilevel_respa = 0;
}
/* ---------------------------------------------------------------------- */
@@ -73,8 +76,12 @@ int FixViscous::setmask()
void FixViscous::init()
{
- if (strstr(update->integrate_style,"respa"))
- nlevels_respa = ((Respa *) update->integrate)->nlevels;
+ int max_respa = 0;
+
+ if (strstr(update->integrate_style,"respa")) {
+ ilevel_respa = max_respa = ((Respa *) update->integrate)->nlevels-1;
+ if (respa_level >= 0) ilevel_respa = MIN(respa_level,max_respa);
+ }
}
/* ---------------------------------------------------------------------- */
@@ -84,9 +91,9 @@ void FixViscous::setup(int vflag)
if (strstr(update->integrate_style,"verlet"))
post_force(vflag);
else {
- ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
- post_force_respa(vflag,nlevels_respa-1,0);
- ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
+ ((Respa *) update->integrate)->copy_flevel_f(ilevel_respa);
+ post_force_respa(vflag,ilevel_respa,0);
+ ((Respa *) update->integrate)->copy_f_flevel(ilevel_respa);
}
}
@@ -126,7 +133,7 @@ void FixViscous::post_force(int vflag)
void FixViscous::post_force_respa(int vflag, int ilevel, int iloop)
{
- if (ilevel == nlevels_respa-1) post_force(vflag);
+ if (ilevel == ilevel_respa) post_force(vflag);
}
/* ---------------------------------------------------------------------- */
diff --git a/src/fix_viscous.h b/src/fix_viscous.h
index 8045e29fa468fe28d6788c22b44c70c622b2a004..7d42a66cba46c22df29283761aa9234982b04996 100644
--- a/src/fix_viscous.h
+++ b/src/fix_viscous.h
@@ -38,7 +38,7 @@ class FixViscous : public Fix {
protected:
double *gamma;
- int nlevels_respa;
+ int ilevel_respa;
};
}
diff --git a/src/fix_wall.cpp b/src/fix_wall.cpp
index fcebd77a9a08695b260c2dfec1c65a3114bb96e2..83d34c30697e25d4687a183d785c2157c98c9d65 100644
--- a/src/fix_wall.cpp
+++ b/src/fix_wall.cpp
@@ -42,6 +42,8 @@ FixWall::FixWall(LAMMPS *lmp, int narg, char **arg) :
global_freq = 1;
extscalar = 1;
extvector = 1;
+ respa_level_support = 1;
+ ilevel_respa = 0;
// parse args
@@ -227,8 +229,6 @@ int FixWall::setmask()
void FixWall::init()
{
- dt = update->dt;
-
for (int m = 0; m < nwall; m++) {
if (xstyle[m] == VARIABLE) {
xindex[m] = input->variable->find(xstr[m]);
@@ -257,8 +257,10 @@ void FixWall::init()
for (int m = 0; m < nwall; m++) precompute(m);
- if (strstr(update->integrate_style,"respa"))
- nlevels_respa = ((Respa *) update->integrate)->nlevels;
+ if (strstr(update->integrate_style,"respa")) {
+ ilevel_respa = ((Respa *) update->integrate)->nlevels-1;
+ if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa);
+ }
}
/* ---------------------------------------------------------------------- */
@@ -268,9 +270,9 @@ void FixWall::setup(int vflag)
if (strstr(update->integrate_style,"verlet")) {
if (!fldflag) post_force(vflag);
} else {
- ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
- post_force_respa(vflag,nlevels_respa-1,0);
- ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
+ ((Respa *) update->integrate)->copy_flevel_f(ilevel_respa);
+ post_force_respa(vflag,ilevel_respa,0);
+ ((Respa *) update->integrate)->copy_f_flevel(ilevel_respa);
}
}
@@ -335,7 +337,7 @@ void FixWall::post_force(int vflag)
void FixWall::post_force_respa(int vflag, int ilevel, int iloop)
{
- if (ilevel == nlevels_respa-1) post_force(vflag);
+ if (ilevel == ilevel_respa) post_force(vflag);
}
/* ---------------------------------------------------------------------- */
diff --git a/src/fix_wall.h b/src/fix_wall.h
index 65163dbe073bf4484f10af5290d4084d0c9ced52..82634c349fb6f540505b65d25a4a34f41a539dca 100644
--- a/src/fix_wall.h
+++ b/src/fix_wall.h
@@ -53,8 +53,7 @@ class FixWall : public Fix {
char *estr[6],*sstr[6];
int varflag; // 1 if any wall position,epsilon,sigma is a var
int eflag; // per-wall flag for energy summation
- int nlevels_respa;
- double dt;
+ int ilevel_respa;
int fldflag;
};
diff --git a/src/fix_wall_region.cpp b/src/fix_wall_region.cpp
index ca17304f3633663c711f2dc8f510953f9554d58d..824f92eb89e66825cd7f40f5289d95bde76cbc2b 100644
--- a/src/fix_wall_region.cpp
+++ b/src/fix_wall_region.cpp
@@ -44,6 +44,8 @@ FixWallRegion::FixWallRegion(LAMMPS *lmp, int narg, char **arg) :
global_freq = 1;
extscalar = 1;
extvector = 1;
+ respa_level_support = 1;
+ ilevel_respa = 0;
// parse args
@@ -151,8 +153,10 @@ void FixWallRegion::init()
offset = coeff3*r4inv*r4inv*rinv - coeff4*r2inv*rinv;
}
- if (strstr(update->integrate_style,"respa"))
- nlevels_respa = ((Respa *) update->integrate)->nlevels;
+ if (strstr(update->integrate_style,"respa")) {
+ ilevel_respa = ((Respa *) update->integrate)->nlevels-1;
+ if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa);
+ }
}
/* ---------------------------------------------------------------------- */
@@ -162,9 +166,9 @@ void FixWallRegion::setup(int vflag)
if (strstr(update->integrate_style,"verlet"))
post_force(vflag);
else {
- ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
- post_force_respa(vflag,nlevels_respa-1,0);
- ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
+ ((Respa *) update->integrate)->copy_flevel_f(ilevel_respa);
+ post_force_respa(vflag,ilevel_respa,0);
+ ((Respa *) update->integrate)->copy_f_flevel(ilevel_respa);
}
}
@@ -245,7 +249,7 @@ void FixWallRegion::post_force(int vflag)
void FixWallRegion::post_force_respa(int vflag, int ilevel, int iloop)
{
- if (ilevel == nlevels_respa-1) post_force(vflag);
+ if (ilevel == ilevel_respa) post_force(vflag);
}
/* ---------------------------------------------------------------------- */
diff --git a/src/fix_wall_region.h b/src/fix_wall_region.h
index 96dd6c17dc4d05a2ac3acfecabbbf6531e70c4e9..cecd7faa33ead64f3fc18412e7bb4a2848661dce 100644
--- a/src/fix_wall_region.h
+++ b/src/fix_wall_region.h
@@ -43,7 +43,7 @@ class FixWallRegion : public Fix {
double epsilon,sigma,cutoff;
int eflag;
double ewall[4],ewall_all[4];
- int nlevels_respa;
+ int ilevel_respa;
char *idregion;
double coeff1,coeff2,coeff3,coeff4,offset;
diff --git a/src/improper.cpp b/src/improper.cpp
index 3b5453de59e100e136ee52f96f2386d20e0a5c7e..0b75ea3a229eb8717875ae146e0217ed5abab09f 100644
--- a/src/improper.cpp
+++ b/src/improper.cpp
@@ -68,6 +68,8 @@ void Improper::init()
error->all(FLERR,"Improper coeffs are not set");
for (int i = 1; i <= atom->nimpropertypes; i++)
if (setflag[i] == 0) error->all(FLERR,"All improper coeffs are not set");
+
+ init_style();
}
/* ----------------------------------------------------------------------
diff --git a/src/improper.h b/src/improper.h
index 5e0b525b6d7c8445b230b671ea5e0fcacc2fc7b1..e029233ae9d152459745a90eb094754f1210ed4f 100644
--- a/src/improper.h
+++ b/src/improper.h
@@ -40,6 +40,7 @@ class Improper : protected Pointers {
Improper(class LAMMPS *);
virtual ~Improper();
virtual void init();
+ virtual void init_style() {}
virtual void compute(int, int) = 0;
virtual void settings(int, char **) {}
virtual void coeff(int, char **) = 0;
diff --git a/src/improper_hybrid.cpp b/src/improper_hybrid.cpp
index 2b7d0cde8d07a1365f8d6de6bdd55e63ba150ef0..d819365b74103a43b8ff3445d350a9f386eb6191 100644
--- a/src/improper_hybrid.cpp
+++ b/src/improper_hybrid.cpp
@@ -156,6 +156,15 @@ void ImproperHybrid::allocate()
for (int m = 0; m < nstyles; m++) improperlist[m] = NULL;
}
+/* ---------------------------------------------------------------------- */
+
+void ImproperHybrid::init_style()
+{
+ for (int i = 0; i < nstyles; i++)
+ styles[i]->init_style();
+}
+
+
/* ----------------------------------------------------------------------
create one improper style for each arg in list
------------------------------------------------------------------------- */
diff --git a/src/improper_hybrid.h b/src/improper_hybrid.h
index cab5184891af7aece62a4e86039af15a9c91d554..dbba347cf947e96238e3b587ef933f1380627ce8 100644
--- a/src/improper_hybrid.h
+++ b/src/improper_hybrid.h
@@ -33,6 +33,7 @@ class ImproperHybrid : public Improper {
ImproperHybrid(class LAMMPS *);
~ImproperHybrid();
+ void init_style();
void compute(int, int);
void settings(int, char **);
void coeff(int, char **);
diff --git a/src/min.cpp b/src/min.cpp
index 5cdd5428de2b69dbb8320ad512977e5ec110b7eb..cce07e1b8bdec233a8b0007c922251b5d2d4d0ce 100644
--- a/src/min.cpp
+++ b/src/min.cpp
@@ -185,7 +185,8 @@ void Min::setup()
if (comm->me == 0 && screen) {
fprintf(screen,"Setting up %s style minimization ...\n",
update->minimize_style);
- fprintf(screen," Unit style: %s\n", update->unit_style);
+ fprintf(screen," Unit style : %s\n", update->unit_style);
+ timer->print_timeout(screen);
}
update->setupflag = 1;
@@ -397,7 +398,7 @@ void Min::run(int n)
// add ntimestep to all computes that store invocation times
// since are hardwiring call to thermo/dumps and computes may not be ready
- if (stop_condition) {
+ if (stop_condition != MAXITER) {
update->nsteps = niter;
if (update->restrict_output == 0) {
@@ -798,6 +799,7 @@ char *Min::stopstrings(int n)
"forces are zero",
"quadratic factors are zero",
"trust region too small",
- "HFTN minimizer error"};
+ "HFTN minimizer error",
+ "walltime limit reached"};
return (char *) strings[n];
}
diff --git a/src/min.h b/src/min.h
index 417cef1c79cfa2273c9108c2c2f125121de6616d..639f87ed66f2aba9d546e440ead3fde35d4d3740 100644
--- a/src/min.h
+++ b/src/min.h
@@ -46,6 +46,10 @@ class Min : protected Pointers {
virtual void reset_vectors() = 0;
virtual int iterate(int) = 0;
+ // possible return values of iterate() method
+ enum{MAXITER,MAXEVAL,ETOL,FTOL,DOWNHILL,ZEROALPHA,ZEROFORCE,
+ ZEROQUAD,TRSMALL,INTERROR,TIMEOUT};
+
protected:
int eflag,vflag; // flags for energy/virial computation
int virial_style; // compute virial explicitly or implicitly
diff --git a/src/min_cg.cpp b/src/min_cg.cpp
index 9bf41d16877dfe90e06c33247e7cd844825dbf40..eda59bd3d4fd19fe586ccb459b327d337ff4ebc8 100644
--- a/src/min_cg.cpp
+++ b/src/min_cg.cpp
@@ -26,10 +26,6 @@ using namespace LAMMPS_NS;
#define EPS_ENERGY 1.0e-8
-// same as in other min classes
-
-enum{MAXITER,MAXEVAL,ETOL,FTOL,DOWNHILL,ZEROALPHA,ZEROFORCE,ZEROQUAD};
-
/* ---------------------------------------------------------------------- */
MinCG::MinCG(LAMMPS *lmp) : MinLineSearch(lmp) {}
@@ -67,6 +63,9 @@ int MinCG::iterate(int maxiter)
for (int iter = 0; iter < maxiter; iter++) {
+ if (timer->check_timeout(niter))
+ return TIMEOUT;
+
ntimestep = ++update->ntimestep;
niter++;
diff --git a/src/min_fire.cpp b/src/min_fire.cpp
index 186c8232a826e45bf707fa8706f213e202587b7e..698b00b7a6781a021ce133e41e88ebe5ef3d8eff 100644
--- a/src/min_fire.cpp
+++ b/src/min_fire.cpp
@@ -27,10 +27,6 @@ using namespace LAMMPS_NS;
#define EPS_ENERGY 1.0e-8
-// same as in other min classes
-
-enum{MAXITER,MAXEVAL,ETOL,FTOL,DOWNHILL,ZEROALPHA,ZEROFORCE,ZEROQUAD};
-
#define DELAYSTEP 5
#define DT_GROW 1.1
#define DT_SHRINK 0.5
@@ -93,6 +89,9 @@ int MinFire::iterate(int maxiter)
for (int iter = 0; iter < maxiter; iter++) {
+ if (timer->check_timeout(niter))
+ return TIMEOUT;
+
ntimestep = ++update->ntimestep;
niter++;
diff --git a/src/min_hftn.cpp b/src/min_hftn.cpp
index 78dbe2836b122b0c0b88f3b6ddba82d5b2b54396..d08df65057f0028f59e1981a965939b790430a88 100644
--- a/src/min_hftn.cpp
+++ b/src/min_hftn.cpp
@@ -42,12 +42,12 @@ using namespace LAMMPS_NS;
//---- CONSTANTS MAP TO stopstrings DECLARED IN Min.run (min.cpp).
-static const int STOP_MAX_ITER = 0; //-- MAX ITERATIONS EXCEEDED
-static const int STOP_MAX_FORCE_EVALS = 1; //-- MAX FORCE EVALUATIONS EXCEEDED
-static const int STOP_ENERGY_TOL = 2; //-- STEP DID NOT CHANGE ENERGY
-static const int STOP_FORCE_TOL = 3; //-- CONVERGED TO DESIRED FORCE TOL
-static const int STOP_TR_TOO_SMALL = 8; //-- TRUST REGION TOO SMALL
-static const int STOP_ERROR = 9; //-- INTERNAL ERROR
+static const int STOP_MAX_ITER = Min::MAXITER; //-- MAX ITERATIONS EXCEEDED
+static const int STOP_MAX_FORCE_EVALS = Min::MAXEVAL; //-- MAX FORCE EVALUATIONS EXCEEDED
+static const int STOP_ENERGY_TOL = Min::ETOL; //-- STEP DID NOT CHANGE ENERGY
+static const int STOP_FORCE_TOL = Min::FTOL; //-- CONVERGED TO DESIRED FORCE TOL
+static const int STOP_TR_TOO_SMALL = Min::TRSMALL; //-- TRUST REGION TOO SMALL
+static const int STOP_ERROR = Min::INTERROR; //-- INTERNAL ERROR
static const int NO_CGSTEP_BECAUSE_F_TOL_SATISFIED = 0;
static const int CGSTEP_NEWTON = 1;
@@ -292,6 +292,10 @@ int MinHFTN::execute_hftn_(const bool bPrintProgress,
double dCurrentEnergy = dInitialEnergy;
double dCurrentForce2 = dInitialForce2;
for (niter = 0; niter < update->nsteps; niter++) {
+
+ if (timer->check_timeout(niter))
+ return(Min::TIMEOUT);
+
(update->ntimestep)++;
//---- CALL THE INNER LOOP TO GET THE NEXT TRUST REGION STEP.
diff --git a/src/min_linesearch.cpp b/src/min_linesearch.cpp
index 3b5566e8971d4732a2d90cf0b3784d1081fba292..ddc0c64321ccd84372e8bb2a2420bd732667b955 100644
--- a/src/min_linesearch.cpp
+++ b/src/min_linesearch.cpp
@@ -52,10 +52,6 @@ using namespace LAMMPS_NS;
#define EMACH 1.0e-8
#define EPS_QUAD 1.0e-28
-// same as in other min classes
-
-enum{MAXITER,MAXEVAL,ETOL,FTOL,DOWNHILL,ZEROALPHA,ZEROFORCE,ZEROQUAD};
-
/* ---------------------------------------------------------------------- */
MinLineSearch::MinLineSearch(LAMMPS *lmp) : Min(lmp)
diff --git a/src/min_quickmin.cpp b/src/min_quickmin.cpp
index 48aee919e3a7031509ce6ef813008ad17ea06e5a..491387406d10a6f1407cf0b970bacfdf5dc25865 100644
--- a/src/min_quickmin.cpp
+++ b/src/min_quickmin.cpp
@@ -28,10 +28,6 @@ using namespace LAMMPS_NS;
#define EPS_ENERGY 1.0e-8
-// same as in other min classes
-
-enum{MAXITER,MAXEVAL,ETOL,FTOL,DOWNHILL,ZEROALPHA,ZEROFORCE,ZEROQUAD};
-
#define DELAYSTEP 5
/* ---------------------------------------------------------------------- */
@@ -88,6 +84,9 @@ int MinQuickMin::iterate(int maxiter)
for (int iter = 0; iter < maxiter; iter++) {
+ if (timer->check_timeout(niter))
+ return TIMEOUT;
+
ntimestep = ++update->ntimestep;
niter++;
diff --git a/src/min_sd.cpp b/src/min_sd.cpp
index 5f97a117e3f14e54a7c6d9ac4dc44b6e0b2be670..73a3867d8a29cce5667ace0a367a2689b9b6f3ec 100644
--- a/src/min_sd.cpp
+++ b/src/min_sd.cpp
@@ -25,10 +25,6 @@ using namespace LAMMPS_NS;
#define EPS_ENERGY 1.0e-8
-// same as in other min classes
-
-enum{MAXITER,MAXEVAL,ETOL,FTOL,DOWNHILL,ZEROALPHA,ZEROFORCE,ZEROQUAD};
-
/* ---------------------------------------------------------------------- */
MinSD::MinSD(LAMMPS *lmp) : MinLineSearch(lmp) {}
@@ -58,6 +54,9 @@ int MinSD::iterate(int maxiter)
for (int iter = 0; iter < maxiter; iter++) {
+ if (timer->check_timeout(niter))
+ return TIMEOUT;
+
ntimestep = ++update->ntimestep;
niter++;
diff --git a/src/minimize.cpp b/src/minimize.cpp
index 89d63f5e72cca4c99d2bc3dba5b4b86c168d1770..e4e89804bf6e07ef9d63533d910b40b654423def 100644
--- a/src/minimize.cpp
+++ b/src/minimize.cpp
@@ -36,6 +36,9 @@ void Minimize::command(int narg, char **arg)
if (domain->box_exist == 0)
error->all(FLERR,"Minimize command before simulation box is defined");
+ // ignore minimize command, if walltime limit was already reached
+ if (timer->is_timeout()) return;
+
update->etol = force->numeric(FLERR,arg[0]);
update->ftol = force->numeric(FLERR,arg[1]);
update->nsteps = force->inumeric(FLERR,arg[2]);
@@ -54,6 +57,7 @@ void Minimize::command(int narg, char **arg)
error->all(FLERR,"Cannot yet use minimize with Kokkos");
lmp->init();
+ timer->init_timeout();
update->minimize->setup();
timer->init();
diff --git a/src/modify.h b/src/modify.h
index 17d409a46e2ea36f3958377b181d58fddcc05315..5ca70428103338dc9c1adcd7307948d92f0bdf71 100644
--- a/src/modify.h
+++ b/src/modify.h
@@ -24,6 +24,8 @@ namespace LAMMPS_NS {
class Modify : protected Pointers {
friend class Info;
friend class FixSRP;
+ friend class Respa;
+ friend class RespaOMP;
public:
int nfix,maxfix;
int n_initial_integrate,n_post_integrate,n_pre_exchange,n_pre_neighbor;
diff --git a/src/neigh_shardlow.cpp b/src/neigh_shardlow.cpp
index 0163e6eb4d6dbdf3e98d113181cf0f834d5d96d5..3a1246684b2d9c46eea40e966aa59ed6f98391de 100644
--- a/src/neigh_shardlow.cpp
+++ b/src/neigh_shardlow.cpp
@@ -30,93 +30,6 @@
using namespace LAMMPS_NS;
-/* ----------------------------------------------------------------------
- convert atom coords into the ssa active interaction region number
-------------------------------------------------------------------------- */
-
-int Neighbor::coord2ssa_airnum(double *x)
-{
- int ix, iy, iz;
-
- ix = iy = iz = 0;
- if (x[2] < domain->sublo[2]) iz = -1;
- if (x[2] > domain->subhi[2]) iz = 1;
- if (x[1] < domain->sublo[1]) iy = -1;
- if (x[1] > domain->subhi[1]) iy = 1;
- if (x[0] < domain->sublo[0]) ix = -1;
- if (x[0] > domain->subhi[0]) ix = 1;
-
- if(iz < 0) return 0;
-
- if(iz == 0){
- if( iy<0 ) return 0; // bottom left/middle/right
- if( (iy==0) && (ix<0) ) return 0; // left atoms
- if( (iy==0) && (ix==0) ) return 1; // Locally owned atoms
- if( (iy==0) && (ix>0) ) return 3; // Right atoms
- if( (iy>0) && (ix==0) ) return 2; // Top-middle atoms
- if( (iy>0) && (ix!=0) ) return 4; // Top-right and top-left atoms
- } else if(iz > 0) {
- if((ix==0) && (iy==0)) return 5; // Back atoms
- if((ix==0) && (iy!=0)) return 6; // Top-back and bottom-back atoms
- if((ix!=0) && (iy==0)) return 7; // Left-back and right-back atoms
- if((ix!=0) && (iy!=0)) return 8; // Back corner atoms
- }
-
- return 0;
-}
-
-/* ----------------------------------------------------------------------
- * assign owned and ghost atoms their ssa active interaction region numbers
- * Called in the pre_neighbor and setup_pre_neighbor fix stages
-------------------------------------------------------------------------- */
-
-void Neighbor::assign_ssa_airnums()
-{
- int i,ibin;
- double **x = atom->x;
- int *mask = atom->mask;
- int nlocal = atom->nlocal;
- int nall = nlocal + atom->nghost;
-
- if (nall > len_ssa_airnum) {
- len_ssa_airnum = nall;
- memory->destroy(ssa_airnum);
- memory->create(ssa_airnum,len_ssa_airnum,"ssa_airnum");
- }
-
- // bin in reverse order so linked list will be in forward order
-
- if (includegroup) {
- int bitmask = group->bitmask[includegroup];
- for (i = nall-1; i >= nlocal; i--) {
- if (mask[i] & bitmask) {
- ibin = coord2ssa_airnum(x[i]);
- } else {
- ibin = 0;
- }
- ssa_airnum[i] = ibin;
- }
- // All the local excluded atoms are in the zero airnum
- ibin = 0;
- for (i = atom->nlocal-1; i >= atom->nfirst; i--) {
- ssa_airnum[i] = ibin;
- }
- } else {
- for (i = nall-1; i >= nlocal; i--) {
- ibin = coord2ssa_airnum(x[i]);
- ssa_airnum[i] = ibin;
- }
- }
-
- // All the local included atoms are in the same airnum (#1)
- if (i >= 0) {
- ibin = coord2ssa_airnum(x[i]);
- do {
- ssa_airnum[i] = ibin;
- } while (--i >= 0);
- }
-}
-
/* ----------------------------------------------------------------------
routines to create a stencil = list of bin offsets
stencil = bins whose closest corner to central bin is within cutoff
@@ -207,6 +120,7 @@ void Neighbor::half_from_full_newton_ssa(NeighList *list)
int *neighptr,*jlist;
int nlocal = atom->nlocal;
+ int *ssaAIR = atom->ssaAIR;
int *ilist = list->ilist;
int *numneigh = list->numneigh;
@@ -240,7 +154,7 @@ void Neighbor::half_from_full_newton_ssa(NeighList *list)
if (j < nlocal) {
if (i > j) continue;
} else {
- if (ssa_airnum[j] <= 0) continue;
+ if (ssaAIR[j] <= 0) continue;
}
neighptr[n++] = joriginal;
}
@@ -279,6 +193,7 @@ void Neighbor::half_bin_newton_ssa(NeighList *list)
int **nspecial = atom->nspecial;
int nlocal = atom->nlocal;
int nall = nlocal + atom->nghost;
+ int *ssaAIR = atom->ssaAIR;
int *molindex = atom->molindex;
int *molatom = atom->molatom;
@@ -297,55 +212,57 @@ void Neighbor::half_bin_newton_ssa(NeighList *list)
int inum = 0;
+ if (binatomflag) { /* only false in Neighbor::build_one */
/* ----------------------------------------------------------------------
bin owned and ghost atoms for use by Shardlow Splitting Algorithm
exclude ghost atoms that are not in the Active Interaction Regions (AIR)
------------------------------------------------------------------------- */
- if (mbins > maxhead_ssa) {
- maxhead_ssa = mbins;
- memory->destroy(gbinhead_ssa);
- memory->destroy(binhead_ssa);
- memory->create(binhead_ssa,maxhead_ssa,"binhead_ssa");
- memory->create(gbinhead_ssa,maxhead_ssa,"gbinhead_ssa");
- }
- for (i = 0; i < mbins; i++) {
- gbinhead_ssa[i] = -1;
- binhead_ssa[i] = -1;
- }
+ if (mbins > maxhead_ssa) {
+ maxhead_ssa = mbins;
+ memory->destroy(gbinhead_ssa);
+ memory->destroy(binhead_ssa);
+ memory->create(binhead_ssa,maxhead_ssa,"binhead_ssa");
+ memory->create(gbinhead_ssa,maxhead_ssa,"gbinhead_ssa");
+ }
+ for (i = 0; i < mbins; i++) {
+ gbinhead_ssa[i] = -1;
+ binhead_ssa[i] = -1;
+ }
- if (maxbin > maxbin_ssa) {
- maxbin_ssa = maxbin;
- memory->destroy(bins_ssa);
- memory->create(bins_ssa,maxbin_ssa,"bins_ssa");
- }
+ if (maxbin > maxbin_ssa) {
+ maxbin_ssa = maxbin;
+ memory->destroy(bins_ssa);
+ memory->create(bins_ssa,maxbin_ssa,"bins_ssa");
+ }
- // bin in reverse order so linked list will be in forward order
+ // bin in reverse order so linked list will be in forward order
- if (includegroup) {
- int bitmask = group->bitmask[includegroup];
- for (i = nall-1; i >= nlocal; i--) {
- if (ssa_airnum[i] <= 0) continue; // skip ghost atoms not in AIR
- if (mask[i] & bitmask) {
+ if (includegroup) {
+ int bitmask = group->bitmask[includegroup];
+ for (i = nall-1; i >= nlocal; i--) {
+ if (ssaAIR[i] <= 0) continue; // skip ghost atoms not in AIR
+ if (mask[i] & bitmask) {
+ ibin = coord2bin(x[i]);
+ bins_ssa[i] = gbinhead_ssa[ibin];
+ gbinhead_ssa[ibin] = i;
+ }
+ }
+ nlocal = atom->nfirst; // This is important for the code that follows!
+ } else {
+ for (i = nall-1; i >= nlocal; i--) {
+ if (ssaAIR[i] <= 0) continue; // skip ghost atoms not in AIR
ibin = coord2bin(x[i]);
bins_ssa[i] = gbinhead_ssa[ibin];
gbinhead_ssa[ibin] = i;
}
}
- nlocal = atom->nfirst; // This is important for the code that follows!
- } else {
- for (i = nall-1; i >= nlocal; i--) {
- if (ssa_airnum[i] <= 0) continue; // skip ghost atoms not in AIR
+ for (i = nlocal-1; i >= 0; i--) {
ibin = coord2bin(x[i]);
- bins_ssa[i] = gbinhead_ssa[ibin];
- gbinhead_ssa[ibin] = i;
+ bins_ssa[i] = binhead_ssa[ibin];
+ binhead_ssa[ibin] = i;
}
- }
- for (i = nlocal-1; i >= 0; i--) {
- ibin = coord2bin(x[i]);
- bins_ssa[i] = binhead_ssa[ibin];
- binhead_ssa[ibin] = i;
- }
+ } /* else reuse previous binning. See Neighbor::build_one comment. */
ipage->reset();
@@ -429,6 +346,7 @@ void Neighbor::half_bin_newton_ssa(NeighList *list)
// loop over AIR ghost atoms in all bins in "full" stencil
// Note: the non-AIR ghost atoms have already been filtered out
+ // That is a significant time savings because of the "full" stencil
for (k = 0; k < maxstencil; k++) {
if (stencil[k] > mbins) break; /* Check if ghost stencil bins are exhausted */
for (j = gbinhead_ssa[ibin+stencil[k]]; j >= 0; j = bins_ssa[j]) {
diff --git a/src/neighbor.cpp b/src/neighbor.cpp
index f1766dd81e3abb1eb5d24f8306027108985acf59..f82a20acd944e205b3091bdf161d97b879ece2e9 100644
--- a/src/neighbor.cpp
+++ b/src/neighbor.cpp
@@ -104,10 +104,8 @@ Neighbor::Neighbor(LAMMPS *lmp) : Pointers(lmp)
maxbin = 0;
bins = NULL;
- // SSA AIR binning
+ // USER-DPD SSA AIR binning
- len_ssa_airnum = 0;
- ssa_airnum = NULL;
maxbin_ssa = 0;
bins_ssa = NULL;
binhead_ssa = NULL;
@@ -189,7 +187,6 @@ Neighbor::~Neighbor()
memory->destroy(gbinhead_ssa);
memory->destroy(binhead_ssa);
memory->destroy(bins_ssa);
- memory->destroy(ssa_airnum);
memory->destroy(ex1_type);
memory->destroy(ex2_type);
@@ -390,6 +387,15 @@ void Neighbor::init()
maxbin = maxhead = 0;
binhead = NULL;
bins = NULL;
+
+ // for USER-DPD Shardlow Splitting Algorithm (SSA)
+ memory->destroy(bins_ssa);
+ memory->destroy(binhead_ssa);
+ memory->destroy(gbinhead_ssa);
+ maxbin_ssa = maxhead_ssa = 0;
+ bins_ssa = NULL;
+ binhead_ssa = NULL;
+ gbinhead_ssa = NULL;
}
// 1st time allocation of xhold and bins
@@ -2202,10 +2208,11 @@ bigint Neighbor::memory_usage()
if (style != NSQ) {
bytes += memory->usage(bins,maxbin);
bytes += memory->usage(binhead,maxhead);
+ bytes += memory->usage(bins_ssa,maxbin_ssa);
+ bytes += memory->usage(binhead_ssa,maxhead_ssa);
+ bytes += memory->usage(gbinhead_ssa,maxhead_ssa);
}
- bytes += memory->usage(ssa_airnum,len_ssa_airnum);
-
for (int i = 0; i < nrequest; i++)
if (lists[i]) bytes += lists[i]->memory_usage();
diff --git a/src/neighbor.h b/src/neighbor.h
index 9496db647bc565ef036ff5dad00a1764b1a96ccb..b44e0fde00f5a43033e6c40faffeaeeef1295f3f 100644
--- a/src/neighbor.h
+++ b/src/neighbor.h
@@ -71,10 +71,6 @@ class Neighbor : protected Pointers {
int cluster_check; // 1 if check bond/angle/etc satisfies minimg
- // USER-DPD package
-
- int *ssa_airnum; // AIR number of each atom for SSA in USER-DPD
-
// methods
Neighbor(class LAMMPS *);
@@ -95,10 +91,6 @@ class Neighbor : protected Pointers {
int exclude_setting();
void exclusion_group_group_delete(int, int); // rm a group-group exclusion
- // USER-DPD package
-
- void assign_ssa_airnums(); // set ssa_airnum values
-
protected:
int me,nprocs;
@@ -184,7 +176,6 @@ class Neighbor : protected Pointers {
// USER-DPD package
- int len_ssa_airnum; // length of ssa_airnum array
int *bins_ssa; // ptr to next atom in each bin used by SSA
int maxbin_ssa; // size of bins array used by SSA
int *binhead_ssa; // ptr to 1st atom in each bin used by SSA
@@ -329,8 +320,6 @@ class Neighbor : protected Pointers {
// SSA neighboring for USER-DPD
- int coord2ssa_airnum(double *); // map atom coord to an AIR number
-
void half_bin_newton_ssa(NeighList *);
void half_from_full_newton_ssa(class NeighList *);
void stencil_half_bin_2d_ssa(class NeighList *, int, int, int);
diff --git a/src/random_mars.cpp b/src/random_mars.cpp
index 96398b170ccbf672f1f413cc6481fcd73ed918a8..fbac44501dc7e10e8dd0789fca4f3c153ff41984 100644
--- a/src/random_mars.cpp
+++ b/src/random_mars.cpp
@@ -97,13 +97,11 @@ double RanMars::gaussian()
double first,v1,v2,rsq,fac;
if (!save) {
- int again = 1;
- while (again) {
+ do {
v1 = 2.0*uniform()-1.0;
v2 = 2.0*uniform()-1.0;
rsq = v1*v1 + v2*v2;
- if (rsq < 1.0 && rsq != 0.0) again = 0;
- }
+ } while ((rsq >= 1.0) || (rsq == 0.0));
fac = sqrt(-2.0*log(rsq)/rsq);
second = v1*fac;
first = v2*fac;
diff --git a/src/random_park.cpp b/src/random_park.cpp
index 28580c5cf5b99b05517059ac9a4da82bff17920c..80460a39fd26168a8e24b034795e87f2e97d325f 100644
--- a/src/random_park.cpp
+++ b/src/random_park.cpp
@@ -57,13 +57,11 @@ double RanPark::gaussian()
double first,v1,v2,rsq,fac;
if (!save) {
- int again = 1;
- while (again) {
+ do {
v1 = 2.0*uniform()-1.0;
v2 = 2.0*uniform()-1.0;
rsq = v1*v1 + v2*v2;
- if (rsq < 1.0 && rsq != 0.0) again = 0;
- }
+ } while ((rsq >= 1.0) || (rsq == 0.0));
fac = sqrt(-2.0*log(rsq)/rsq);
second = v1*fac;
first = v2*fac;
diff --git a/src/respa.cpp b/src/respa.cpp
index b96e5e79cb5e2b34d51d2268ee829883637d2e5c..6c1be33921f5a457fdf6df03568854a56e61f50a 100644
--- a/src/respa.cpp
+++ b/src/respa.cpp
@@ -394,7 +394,18 @@ void Respa::setup()
fprintf(screen,"Setting up r-RESPA run ...\n");
fprintf(screen," Unit style : %s\n", update->unit_style);
fprintf(screen," Current step : " BIGINT_FORMAT "\n", update->ntimestep);
- fprintf(screen," OuterTime step: %g\n", update->dt);
+ fprintf(screen," Time steps :");
+ for (int ilevel=0; ilevel < nlevels; ++ilevel)
+ fprintf(screen," %d:%g",ilevel+1, step[ilevel]);
+ fprintf(screen,"\n r-RESPA fixes :");
+ for (int l=0; l < modify->n_post_force_respa; ++l) {
+ Fix *f = modify->fix[modify->list_post_force_respa[l]];
+ if (f->respa_level >= 0)
+ fprintf(screen," %d:%s[%s]",
+ MIN(f->respa_level+1,nlevels),f->style,f->id);
+ }
+ fprintf(screen,"\n");
+ timer->print_timeout(screen);
}
update->setupflag = 1;
@@ -549,6 +560,10 @@ void Respa::run(int n)
bigint ntimestep;
for (int i = 0; i < n; i++) {
+ if (timer->check_timeout(i)) {
+ update->nsteps = i;
+ break;
+ }
ntimestep = ++update->ntimestep;
ev_set(ntimestep);
@@ -718,7 +733,7 @@ void Respa::recurse(int ilevel)
if (modify->n_pre_reverse) {
modify->pre_reverse(eflag,vflag);
timer->stamp(Timer::MODIFY);
- }
+ }
if (newton[ilevel]) {
comm->reverse_comm();
diff --git a/src/run.cpp b/src/run.cpp
index b120d2a14cb758f3e3a06846e300c01f1f5938cc..f238e7852b43ecc144fae4e23db81ecd3420c76d 100644
--- a/src/run.cpp
+++ b/src/run.cpp
@@ -42,6 +42,9 @@ void Run::command(int narg, char **arg)
if (domain->box_exist == 0)
error->all(FLERR,"Run command before simulation box is defined");
+ // ignore run command, if walltime limit was already reached
+ if (timer->is_timeout()) return;
+
bigint nsteps_input = force->bnumeric(FLERR,arg[0]);
// parse optional args
@@ -130,6 +133,9 @@ void Run::command(int narg, char **arg)
error->all(FLERR,"Run command stop value is before end of run");
}
+ if (!preflag && strstr(update->integrate_style,"respa"))
+ error->all(FLERR,"Run flag 'pre no' not compatible with r-RESPA");
+
// if nevery, make copies of arg strings that are commands
// required because re-parsing commands via input->one() will wipe out args
@@ -152,6 +158,7 @@ void Run::command(int narg, char **arg)
// if post, do full Finish, else just print time
update->whichflag = 1;
+ timer->init_timeout();
if (nevery == 0) {
update->nsteps = nsteps;
@@ -190,6 +197,9 @@ void Run::command(int narg, char **arg)
int iter = 0;
int nleft = nsteps;
while (nleft > 0 || iter == 0) {
+ if (timer->is_timeout()) break;
+ timer->init_timeout();
+
nsteps = MIN(nleft,nevery);
update->nsteps = nsteps;
diff --git a/src/timer.cpp b/src/timer.cpp
index a075ff3dec7a1c7a2605525003d682ab4e125d6b..165b6de64bd3b094e5dc96d027e8b4f352587255 100644
--- a/src/timer.cpp
+++ b/src/timer.cpp
@@ -13,9 +13,11 @@
#include <mpi.h>
#include <string.h>
+#include <stdlib.h>
#include "timer.h"
#include "comm.h"
#include "error.h"
+#include "force.h"
#include "memory.h"
#ifdef _WIN32
@@ -26,8 +28,37 @@
#include <sys/resource.h>
#endif
+#include <time.h>
+
using namespace LAMMPS_NS;
+// convert a timespec ([[HH:]MM:]SS) to seconds
+// the strings "off" and "unlimited" result in -1;
+
+static double timespec2seconds(char *timespec)
+{
+ double vals[3];
+ char *num;
+ int i = 0;
+
+ // first handle allowed textual inputs
+ if (strcmp(timespec,"off") == 0) return -1;
+ if (strcmp(timespec,"unlimited") == 0) return -1;
+
+ vals[0] = vals[1] = vals[2] = 0;
+
+ num = strtok(timespec,":");
+ while ((num != NULL) && (i < 3)) {
+ vals[i] = atoi(num);
+ ++i;
+ num = strtok(NULL,":");
+ }
+
+ if (i == 3) return (vals[0]*60 + vals[1])*60 + vals[2];
+ else if (i == 2) return vals[0]*60 + vals[1];
+ else return vals[0];
+}
+
// Return the CPU time for the current process in seconds very
// much in the same way as MPI_Wtime() returns the wall time.
@@ -64,6 +95,9 @@ Timer::Timer(LAMMPS *lmp) : Pointers(lmp)
{
_level = NORMAL;
_sync = OFF;
+ _timeout = -1.0;
+ _checkfreq = 10;
+ _nextcheck = -1;
}
/* ---------------------------------------------------------------------- */
@@ -176,12 +210,65 @@ void Timer::set_wall(enum ttype which, double newtime)
wall_array[which] = newtime;
}
+/* ---------------------------------------------------------------------- */
+
+void Timer::init_timeout()
+{
+ if (_timeout < 0)
+ _nextcheck = -1;
+ else
+ _nextcheck = _checkfreq;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void Timer::print_timeout(FILE *fp)
+{
+ if (!fp) return;
+
+ // format timeout setting
+ if (_timeout > 0) {
+ // time since init_timeout()
+ const double d = MPI_Wtime() - timeout_start;
+ // remaining timeout in seconds
+ int s = _timeout - d;
+ // remaining 1/100ths of seconds
+ const int hs = 100*((_timeout - d) - s);
+ // breaking s down into second/minutes/hours
+ const int seconds = s % 60;
+ s = (s - seconds) / 60;
+ const int minutes = s % 60;
+ const int hours = (s - minutes) / 60;
+ fprintf(fp," Walltime left : %d:%02d:%02d.%02d\n",
+ hours,minutes,seconds,hs);
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+bool Timer::_check_timeout()
+{
+ double walltime = MPI_Wtime() - timeout_start;
+ // broadcast time to insure all ranks act the same.
+ MPI_Bcast(&walltime,1,MPI_DOUBLE,0,world);
+
+ if (walltime < _timeout) {
+ _nextcheck += _checkfreq;
+ return false;
+ } else {
+ if (comm->me == 0)
+ error->warning(FLERR,"Wall time limit reached");
+ _timeout = 0.0;
+ return true;
+ }
+}
+
/* ----------------------------------------------------------------------
modify parameters of the Timer class
------------------------------------------------------------------------- */
static const char *timer_style[] = { "off", "loop", "normal", "full" };
static const char *timer_mode[] = { "nosync", "(dummy)", "sync" };
-static const char timer_fmt[] = "New timer settings: style=%s mode=%s\n";
+static const char timer_fmt[] = "New timer settings: style=%s mode=%s timeout=%s\n";
void Timer::modify_params(int narg, char **arg)
{
@@ -199,14 +286,37 @@ void Timer::modify_params(int narg, char **arg)
_sync = OFF;
} else if (strcmp(arg[iarg],timer_mode[NORMAL]) == 0) {
_sync = NORMAL;
+ } else if (strcmp(arg[iarg],"timeout") == 0) {
+ ++iarg;
+ if (iarg < narg) {
+ _timeout = timespec2seconds(arg[iarg]);
+ } else error->all(FLERR,"Illegal timers command");
+ } else if (strcmp(arg[iarg],"every") == 0) {
+ ++iarg;
+ if (iarg < narg) {
+ _checkfreq = force->inumeric(FLERR,arg[iarg]);
+ if (_checkfreq <= 0)
+ error->all(FLERR,"Illegal timers command");
+ } else error->all(FLERR,"Illegal timers command");
} else error->all(FLERR,"Illegal timers command");
++iarg;
}
+ timeout_start = MPI_Wtime();
if (comm->me == 0) {
+
+ // format timeout setting
+ char timebuf[32];
+ if (_timeout < 0) strcpy(timebuf,"off");
+ else {
+ time_t tv = _timeout;
+ struct tm *tm = gmtime(&tv);
+ strftime(timebuf,32,"%H:%M:%S",tm);
+ }
+
if (screen)
- fprintf(screen,timer_fmt,timer_style[_level],timer_mode[_sync]);
+ fprintf(screen,timer_fmt,timer_style[_level],timer_mode[_sync],timebuf);
if (logfile)
- fprintf(logfile,timer_fmt,timer_style[_level],timer_mode[_sync]);
+ fprintf(logfile,timer_fmt,timer_style[_level],timer_mode[_sync],timebuf);
}
}
diff --git a/src/timer.h b/src/timer.h
index ef5306cb2ec23dbce42cf40ec7bff8f04aa5ba7b..aec6dd7fa9e488b965e0aec0e86b11640f4293ed 100644
--- a/src/timer.h
+++ b/src/timer.h
@@ -47,6 +47,9 @@ class Timer : protected Pointers {
bool has_full() const { return (_level >= FULL); }
bool has_sync() const { return (_sync != OFF); }
+ // flag if wallclock time is expired
+ bool is_timeout() const { return (_timeout == 0.0); }
+
double elapsed(enum ttype);
double cpu(enum ttype);
@@ -57,6 +60,21 @@ class Timer : protected Pointers {
void set_wall(enum ttype, double);
+ // initialize timeout timer
+ void init_timeout();
+
+ // trigger enforced timeout
+ void force_timeout() { _timeout = 0.0; };
+
+ // print timeout message
+ void print_timeout(FILE *);
+
+ // check for timeout. inline wrapper around internal
+ // function to reduce overhead in case there is no check.
+ bool check_timeout(int step) {
+ if (_nextcheck != step) return false;
+ else return _check_timeout();
+ }
void modify_params(int, char **);
@@ -65,11 +83,18 @@ class Timer : protected Pointers {
double wall_array[NUM_TIMER];
double previous_cpu;
double previous_wall;
- int _level; // level of detail: off=0,loop=1,normal=2,full=3
- int _sync; // if nonzero, synchronize tasks before setting the timer
-
- // update requested timer array
+ double timeout_start;
+ int _level; // level of detail: off=0,loop=1,normal=2,full=3
+ int _sync; // if nonzero, synchronize tasks before setting the timer
+ int _timeout; // max allowed wall time in seconds. infinity if negative
+ int _checkfreq; // frequency of timeout checking
+ int _nextcheck; // timestep number of next timeout check
+
+ // update one specific timer array
void _stamp(enum ttype);
+
+ // check for timeout
+ bool _check_timeout();
};
}
diff --git a/src/verlet.cpp b/src/verlet.cpp
index 85b8201e8d7ec99effce6e2078845d4bbac18449..ce06a1c45688c03fa6281b195065061eab0d60ba 100644
--- a/src/verlet.cpp
+++ b/src/verlet.cpp
@@ -89,9 +89,10 @@ void Verlet::setup()
{
if (comm->me == 0 && screen) {
fprintf(screen,"Setting up Verlet run ...\n");
- fprintf(screen," Unit style : %s\n", update->unit_style);
- fprintf(screen," Current step: " BIGINT_FORMAT "\n", update->ntimestep);
- fprintf(screen," Time step : %g\n", update->dt);
+ fprintf(screen," Unit style : %s\n", update->unit_style);
+ fprintf(screen," Current step : " BIGINT_FORMAT "\n", update->ntimestep);
+ fprintf(screen," Time step : %g\n", update->dt);
+ timer->print_timeout(screen);
}
update->setupflag = 1;
@@ -229,6 +230,10 @@ void Verlet::run(int n)
else sortflag = 0;
for (int i = 0; i < n; i++) {
+ if (timer->check_timeout(i)) {
+ update->nsteps = i;
+ break;
+ }
ntimestep = ++update->ntimestep;
ev_set(ntimestep);