From d55f968432a462193c4625af7c88f3d605b2fc39 Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Fri, 17 Jun 2016 23:48:15 +0000
Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@15203
f3b2605a-c512-4ea7-a41b-209d697bcdaa
---
src/KOKKOS/verlet_kokkos.cpp | 14 +-
src/MANYBODY/fix_qeq_comb.cpp | 16 +-
src/MANYBODY/fix_qeq_comb.h | 2 +-
src/MISC/fix_efield.cpp | 16 +-
src/MISC/fix_efield.h | 2 +-
src/MISC/fix_orient_bcc.cpp | 16 +-
src/MISC/fix_orient_bcc.h | 2 +-
src/MISC/fix_orient_fcc.cpp | 16 +-
src/MISC/fix_orient_fcc.h | 2 +-
src/RIGID/fix_ehex.cpp | 8 +-
src/RIGID/fix_ehex.h | 2 +-
src/RIGID/fix_rattle.cpp | 11 +-
src/RIGID/fix_rattle.h | 1 -
src/USER-DPD/atom_vec_dpd.cpp | 111 ++--
src/USER-DPD/atom_vec_dpd.h | 7 +-
src/USER-DPD/compute_dpd.cpp | 8 +-
src/USER-DPD/compute_dpd.h | 4 +-
src/USER-DPD/compute_dpd_atom.h | 2 +-
src/USER-DPD/fix_eos_cv.cpp | 20 +-
src/USER-DPD/fix_eos_cv.h | 2 +-
src/USER-DPD/fix_eos_table.cpp | 26 +-
src/USER-DPD/fix_eos_table.h | 8 +-
src/USER-DPD/fix_eos_table_rx.cpp | 88 +--
src/USER-DPD/fix_eos_table_rx.h | 6 +-
src/USER-DPD/fix_rx.cpp | 211 ++++----
src/USER-DPD/fix_rx.h | 4 +-
src/USER-DPD/fix_shardlow.cpp | 505 ++++++++---------
src/USER-DPD/fix_shardlow.h | 16 +-
src/USER-DPD/pair_dpd_fdt.cpp | 10 +-
src/USER-DPD/pair_dpd_fdt.h | 2 +-
src/USER-DPD/pair_dpd_fdt_energy.cpp | 10 +-
src/USER-DPD/pair_dpd_fdt_energy.h | 2 +-
src/USER-DPD/pair_exp6_rx.cpp | 600 ++++++++++-----------
src/USER-DPD/pair_exp6_rx.h | 4 +-
src/USER-DPD/pair_multi_lucy.cpp | 50 +-
src/USER-DPD/pair_multi_lucy.h | 2 +-
src/USER-DPD/pair_multi_lucy_rx.cpp | 140 +++--
src/USER-DPD/pair_multi_lucy_rx.h | 2 +-
src/USER-DPD/pair_table_rx.cpp | 100 ++--
src/USER-DPD/pair_table_rx.h | 4 +-
src/USER-INTEL/improper_cvff_intel.cpp | 4 +-
src/USER-INTEL/improper_cvff_intel.h | 2 +-
src/USER-INTEL/improper_harmonic_intel.cpp | 4 +-
src/USER-INTEL/improper_harmonic_intel.h | 2 +-
src/USER-MISC/fix_addtorque.cpp | 16 +-
src/USER-MISC/fix_addtorque.h | 2 +-
src/USER-MISC/fix_smd.cpp | 35 +-
src/USER-MISC/fix_smd.h | 2 +-
src/USER-MISC/fix_ttm_mod.cpp | 69 +--
src/USER-MISC/fix_ttm_mod.h | 4 +-
src/USER-OMP/fix_gravity_omp.cpp | 2 +-
src/USER-OMP/fix_qeq_comb_omp.cpp | 6 +-
src/USER-OMP/respa_omp.cpp | 13 +-
src/USER-SMD/fix_smd_setvel.cpp | 10 +-
src/atom.cpp | 6 +-
src/atom.h | 3 +-
src/balance.cpp | 8 +-
src/dump_image.cpp | 1 -
src/fix.cpp | 10 +
src/fix.h | 2 +
src/fix_addforce.cpp | 17 +-
src/fix_addforce.h | 2 +-
src/fix_aveforce.cpp | 13 +-
src/fix_aveforce.h | 2 +-
src/fix_drag.cpp | 16 +-
src/fix_drag.h | 2 +-
src/fix_gravity.cpp | 16 +-
src/fix_gravity.h | 2 +-
src/fix_indent.cpp | 16 +-
src/fix_indent.h | 2 +-
src/fix_restrain.cpp | 16 +-
src/fix_restrain.h | 2 +-
src/fix_setforce.cpp | 12 +-
src/fix_setforce.h | 3 +-
src/fix_spring.cpp | 16 +-
src/fix_spring.h | 2 +-
src/fix_spring_chunk.cpp | 16 +-
src/fix_spring_chunk.h | 2 +-
src/fix_spring_rg.cpp | 16 +-
src/fix_spring_rg.h | 2 +-
src/fix_spring_self.cpp | 15 +-
src/fix_spring_self.h | 2 +-
src/fix_viscous.cpp | 19 +-
src/fix_viscous.h | 2 +-
src/fix_wall.cpp | 18 +-
src/fix_wall.h | 3 +-
src/fix_wall_region.cpp | 16 +-
src/fix_wall_region.h | 2 +-
src/improper.cpp | 2 +
src/improper.h | 1 +
src/improper_hybrid.cpp | 9 +
src/improper_hybrid.h | 1 +
src/min.cpp | 8 +-
src/min.h | 4 +
src/min_cg.cpp | 7 +-
src/min_fire.cpp | 7 +-
src/min_hftn.cpp | 16 +-
src/min_linesearch.cpp | 4 -
src/min_quickmin.cpp | 7 +-
src/min_sd.cpp | 7 +-
src/minimize.cpp | 4 +
src/modify.h | 2 +
src/neigh_shardlow.cpp | 162 ++----
src/neighbor.cpp | 19 +-
src/neighbor.h | 11 -
src/random_mars.cpp | 6 +-
src/random_park.cpp | 6 +-
src/respa.cpp | 19 +-
src/run.cpp | 10 +
src/timer.cpp | 116 +++-
src/timer.h | 33 +-
src/verlet.cpp | 11 +-
112 files changed, 1576 insertions(+), 1401 deletions(-)
diff --git a/src/KOKKOS/verlet_kokkos.cpp b/src/KOKKOS/verlet_kokkos.cpp
index d063294023..f3c4e3e07d 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 c2d9cd9281..3c5954677b 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 d627ef669b..384ea4bbcc 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 2151a05092..85af46ed3b 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 93272ba9bf..4a360b6464 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 88db82bd82..24c587645d 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 aecbaca793..2035e85b15 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 df41981f3d..c0c836dca7 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 3d5d1207aa..8bdc098289 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 3d39bba93d..c968887374 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 1e7a0fe482..3ae974fc0d 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 3a5fef69cb..359be6b885 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 7b94d2cf0d..e45dfb1a4a 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 92028b60e6..586ddca02b 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 71ef48d053..077d18485f 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 59716362c7..08aa19bcdf 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 c2e8133cd4..695048eef3 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 e591e4c48e..c77490be76 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 4a8d213362..df764f7082 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 043b8b82d1..f53c3177af 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 b9579cfaf7..98276d0f01 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 7c33cf6059..1623fa21f8 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 f7acada9d5..5f382436e8 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 bb2e507c44..b513e292be 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 10b297753e..ce1cb46183 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 912b75ed5f..97b9bd55cc 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 919689e618..7248df2e67 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 327e6357dc..6421737b3d 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 de313b820e..21c5e420fa 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 a11d963ef9..f7721853be 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 2fd3c0bf55..34e21918d1 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 ef1afed4e8..c296e8f1a3 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 e6aeb9ac7d..d4710e1a8d 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 ece29dcda7..45abfab77a 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 37f235b8e3..501ef2d04c 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 ee6f8c4431..f3c67e4fa4 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 6a349101d5..b90252f66f 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 8600a3d2c6..a85bc6cd75 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 71b48b2776..6b824b415a 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 00b83bafb7..42e27390de 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 6c0735a755..0fb02420b9 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 154bb5bd18..95ccd8f9d2 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 7693793f02..071ff548ea 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 d5b0e4ee0f..4e38383863 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 0e1ffaf0e1..ed5d8875f7 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 513b658c87..49631ec916 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 3f6ca3f58a..b527e4b788 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 687c1900b2..4f86456f99 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 5a58dd8489..f7b81c7206 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 11a68e4485..3415baf5b1 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 b853d844f0..1e65d1f764 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 f8de18498d..0c69876b86 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 d925088eda..738538a209 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 568ffb7977..92c4a00581 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 ac428b9a95..acb3cbce7c 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 4168df32cf..5e6f257789 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 c615db6112..db7361703c 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 b9ff534c66..5e6701ff34 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 165536d834..1a3b6a42c9 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 4bf2609eff..a2aa3782cc 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 2b5c682266..ac8384f2c0 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 d590463a5e..d19c0db0d1 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 443ae6ea57..09a685cb1a 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 2ec3110343..a6228be7df 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 6a28e89806..3fda219126 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 dc0075143d..fa444a4b35 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 a08213e3d0..95f899ae86 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 f37a5bb411..35be23a40a 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 4b85dd2788..a6693e98a2 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 98f8d081ea..c3a8882734 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 97e898494e..a92352e563 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 0a5669a7e8..a6f23f86c6 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 d73ab0db63..08f177ae8c 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 a445f342ed..54d17049f7 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 c90feaaafe..5b96e1e95d 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 f784ad9d34..1bb551068b 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 d6dc36a73c..16b3d7939e 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 c95d9063ca..1453a8e621 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 4145e01047..b7fd33ce77 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 c3a0129843..17337be6fd 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 830af42162..9517cf0c31 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 3eb5945abd..0e5ddda127 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 bc5059ae69..9c82b34d8d 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 8045e29fa4..7d42a66cba 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 fcebd77a9a..83d34c3069 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 65163dbe07..82634c349f 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 ca17304f36..824f92eb89 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 96dd6c17dc..cecd7faa33 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 3b5453de59..0b75ea3a22 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 5e0b525b6d..e029233ae9 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 2b7d0cde8d..d819365b74 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 cab5184891..dbba347cf9 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 5cdd5428de..cce07e1b8b 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 417cef1c79..639f87ed66 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 9bf41d1687..eda59bd3d4 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 186c8232a8..698b00b7a6 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 78dbe2836b..d08df65057 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 3b5566e897..ddc0c64321 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 48aee919e3..491387406d 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 5f97a117e3..73a3867d8a 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 89d63f5e72..e4e89804bf 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 17d409a46e..5ca7042810 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 0163e6eb4d..3a1246684b 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 f1766dd81e..f82a20acd9 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 9496db647b..b44e0fde00 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 96398b170c..fbac44501d 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 28580c5cf5..80460a39fd 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 b96e5e79cb..6c1be33921 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 b120d2a14c..f238e7852b 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 a075ff3dec..165b6de64b 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 ef5306cb2e..aec6dd7fa9 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 85b8201e8d..ce06a1c456 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);
--
GitLab