From fae9ea3b7914baaa01b3d9b7498d0dc446e716a1 Mon Sep 17 00:00:00 2001
From: athomps <athomps@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Sun, 27 Mar 2011 23:34:04 +0000
Subject: [PATCH] Added the option zero yes/no
git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@5857 f3b2605a-c512-4ea7-a41b-209d697bcdaa
---
src/fix_langevin.cpp | 89 ++++++++++++++++++++++++++++++++++++++------
src/fix_langevin.h | 2 +-
2 files changed, 78 insertions(+), 13 deletions(-)
diff --git a/src/fix_langevin.cpp b/src/fix_langevin.cpp
index ae2cb4f2e5..e05d830881 100644
--- a/src/fix_langevin.cpp
+++ b/src/fix_langevin.cpp
@@ -32,6 +32,7 @@
#include "random_mars.h"
#include "memory.h"
#include "error.h"
+#include "group.h"
using namespace LAMMPS_NS;
@@ -71,6 +72,7 @@ FixLangevin::FixLangevin(LAMMPS *lmp, int narg, char **arg) :
for (int i = 1; i <= atom->ntypes; i++) ratio[i] = 1.0;
tally = 0;
+ zerosum = 0;
int iarg = 7;
while (iarg < narg) {
@@ -88,6 +90,12 @@ FixLangevin::FixLangevin(LAMMPS *lmp, int narg, char **arg) :
else if (strcmp(arg[iarg+1],"yes") == 0) tally = 1;
else error->all("Illegal fix langevin command");
iarg += 2;
+ } else if (strcmp(arg[iarg],"zerosum") == 0) {
+ if (iarg+2 > narg) error->all("Illegal fix langevin command");
+ if (strcmp(arg[iarg+1],"no") == 0) zerosum = 0;
+ else if (strcmp(arg[iarg+1],"yes") == 0) zerosum = 1;
+ else error->all("Illegal fix langevin command");
+ iarg += 2;
} else error->all("Illegal fix langevin command");
}
@@ -204,7 +212,20 @@ void FixLangevin::post_force_no_tally()
// computed on current nlocal atoms to remove bias
// test v = 0 since some computes mask non-participating atoms via v = 0
// and added force has extra term not multiplied by v = 0
+ // for ZEROSUM:
+ // Sum random force over all atoms in group
+ // Subtract sum/count from each atom in group
+ double fran[3],fsum[3],fsumall[3];
+ fsum[0] = fsum[1] = fsum[2] = 0.0;
+ bigint count;
+
+ if (zerosum) {
+ count = group->count(igroup);
+ if (count == 0)
+ error->all("Cannot zero Langevin force of 0 atoms");
+ }
+
if (rmass) {
double boltz = force->boltz;
double dt = update->dt;
@@ -218,9 +239,15 @@ void FixLangevin::post_force_no_tally()
gamma2 = sqrt(rmass[i]) * sqrt(24.0*boltz/t_period/dt/mvv2e) / ftm2v;
gamma1 *= 1.0/ratio[type[i]];
gamma2 *= 1.0/sqrt(ratio[type[i]]) * tsqrt;
- f[i][0] += gamma1*v[i][0] + gamma2*(random->uniform()-0.5);
- f[i][1] += gamma1*v[i][1] + gamma2*(random->uniform()-0.5);
- f[i][2] += gamma1*v[i][2] + gamma2*(random->uniform()-0.5);
+ fran[0] = gamma2*(random->uniform()-0.5);
+ fran[1] = gamma2*(random->uniform()-0.5);
+ fran[2] = gamma2*(random->uniform()-0.5);
+ f[i][0] += gamma1*v[i][0] + fran[0];
+ f[i][1] += gamma1*v[i][1] + fran[1];
+ f[i][2] += gamma1*v[i][2] + fran[2];
+ fsum[0] += fran[0];
+ fsum[1] += fran[1];
+ fsum[2] += fran[2];
}
}
@@ -233,27 +260,42 @@ void FixLangevin::post_force_no_tally()
gamma1 *= 1.0/ratio[type[i]];
gamma2 *= 1.0/sqrt(ratio[type[i]]) * tsqrt;
temperature->remove_bias(i,v[i]);
+ fran[0] = gamma2*(random->uniform()-0.5);
+ fran[1] = gamma2*(random->uniform()-0.5);
+ fran[2] = gamma2*(random->uniform()-0.5);
if (v[i][0] != 0.0)
- f[i][0] += gamma1*v[i][0] + gamma2*(random->uniform()-0.5);
+ f[i][0] += gamma1*v[i][0] + fran[0];
if (v[i][1] != 0.0)
- f[i][1] += gamma1*v[i][1] + gamma2*(random->uniform()-0.5);
+ f[i][1] += gamma1*v[i][1] + fran[1];
if (v[i][2] != 0.0)
- f[i][2] += gamma1*v[i][2] + gamma2*(random->uniform()-0.5);
+ f[i][2] += gamma1*v[i][2] + fran[2];
+ fsum[0] += fran[0];
+ fsum[1] += fran[1];
+ fsum[2] += fran[2];
temperature->restore_bias(i,v[i]);
}
}
}
} else {
+
if (which == NOBIAS) {
+
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
gamma1 = gfactor1[type[i]];
gamma2 = gfactor2[type[i]] * tsqrt;
- f[i][0] += gamma1*v[i][0] + gamma2*(random->uniform()-0.5);
- f[i][1] += gamma1*v[i][1] + gamma2*(random->uniform()-0.5);
- f[i][2] += gamma1*v[i][2] + gamma2*(random->uniform()-0.5);
+ fran[0] = gamma2*(random->uniform()-0.5);
+ fran[1] = gamma2*(random->uniform()-0.5);
+ fran[2] = gamma2*(random->uniform()-0.5);
+ f[i][0] += gamma1*v[i][0] + fran[0];
+ f[i][1] += gamma1*v[i][1] + fran[1];
+ f[i][2] += gamma1*v[i][2] + fran[2];
+ fsum[0] += fran[0];
+ fsum[1] += fran[1];
+ fsum[2] += fran[2];
}
+
}
} else if (which == BIAS) {
@@ -263,17 +305,40 @@ void FixLangevin::post_force_no_tally()
gamma1 = gfactor1[type[i]];
gamma2 = gfactor2[type[i]] * tsqrt;
temperature->remove_bias(i,v[i]);
+ fran[0] = gamma2*(random->uniform()-0.5);
+ fran[1] = gamma2*(random->uniform()-0.5);
+ fran[2] = gamma2*(random->uniform()-0.5);
if (v[i][0] != 0.0)
- f[i][0] += gamma1*v[i][0] + gamma2*(random->uniform()-0.5);
+ f[i][0] += gamma1*v[i][0] + fran[0];
if (v[i][1] != 0.0)
- f[i][1] += gamma1*v[i][1] + gamma2*(random->uniform()-0.5);
+ f[i][1] += gamma1*v[i][1] + fran[1];
if (v[i][2] != 0.0)
- f[i][2] += gamma1*v[i][2] + gamma2*(random->uniform()-0.5);
+ f[i][2] += gamma1*v[i][2] + fran[2];
+ fsum[0] += fran[0];
+ fsum[1] += fran[1];
+ fsum[2] += fran[2];
temperature->restore_bias(i,v[i]);
}
}
}
}
+
+ // Set total force to zero
+
+ if (zerosum) {
+ MPI_Allreduce(fsum,fsumall,3,MPI_DOUBLE,MPI_SUM,world);
+ fsumall[0] /= count;
+ fsumall[1] /= count;
+ fsumall[2] /= count;
+ for (int i = 0; i < nlocal; i++) {
+ if (mask[i] & groupbit) {
+ f[i][0] -= fsumall[0];
+ f[i][1] -= fsumall[1];
+ f[i][2] -= fsumall[2];
+ }
+ }
+ }
+
}
/* ---------------------------------------------------------------------- */
diff --git a/src/fix_langevin.h b/src/fix_langevin.h
index 5953079362..1aa3b15678 100644
--- a/src/fix_langevin.h
+++ b/src/fix_langevin.h
@@ -41,7 +41,7 @@ class FixLangevin : public Fix {
double memory_usage();
protected:
- int which,tally;
+ int which,tally,zerosum;
double t_start,t_stop,t_period;
double *gfactor1,*gfactor2,*ratio;
double energy,energy_onestep;
--
GitLab