From e2707ca96f37c10fc982ea99fdde9db68b7c0784 Mon Sep 17 00:00:00 2001
From: sjplimp
Date: Thu, 28 Apr 2011 16:16:54 +0000
Subject: [PATCH 01/21] git-svn-id:
svn://svn.icms.temple.edu/lammps-ro/trunk@6026
f3b2605a-c512-4ea7-a41b-209d697bcdaa
---
src/fix_rigid.cpp | 35 +++++++++++------------------------
src/math_extra.cpp | 6 ++----
src/math_extra.h | 12 ++++++++++++
3 files changed, 25 insertions(+), 28 deletions(-)
diff --git a/src/fix_rigid.cpp b/src/fix_rigid.cpp
index 1099f09d7a..4d0686409d 100644
--- a/src/fix_rigid.cpp
+++ b/src/fix_rigid.cpp
@@ -648,8 +648,7 @@ void FixRigid::init()
sum[ibody][0] += 0.4 * massone * radius[i]*radius[i];
sum[ibody][1] += 0.4 * massone * radius[i]*radius[i];
sum[ibody][2] += 0.4 * massone * radius[i]*radius[i];
- }
- if (eflags[i] & INERTIA_ELLIPSOID) {
+ } else if (eflags[i] & INERTIA_ELLIPSOID) {
shape = ebonus[ellipsoid[i]].shape;
quatatom = ebonus[ellipsoid[i]].quat;
MathExtra::inertia_ellipsoid(shape,quatatom,massone,ivec);
@@ -665,11 +664,12 @@ void FixRigid::init()
MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world);
+ // diagonalize inertia tensor for each body via Jacobi rotations
// inertia = 3 eigenvalues = principal moments of inertia
- // ex_space,ey_space,ez_space = 3 eigenvectors = principal axes of rigid body
-
+ // evectors and exzy_space = 3 evectors = principal axes of rigid body
+
int ierror;
- double ez0,ez1,ez2;
+ double cross[3];
double tensor[3][3],evectors[3][3];
for (ibody = 0; ibody < nbody; ibody++) {
@@ -686,11 +686,9 @@ void FixRigid::init()
ex_space[ibody][0] = evectors[0][0];
ex_space[ibody][1] = evectors[1][0];
ex_space[ibody][2] = evectors[2][0];
-
ey_space[ibody][0] = evectors[0][1];
ey_space[ibody][1] = evectors[1][1];
ey_space[ibody][2] = evectors[2][1];
-
ez_space[ibody][0] = evectors[0][2];
ez_space[ibody][1] = evectors[1][2];
ez_space[ibody][2] = evectors[2][2];
@@ -706,21 +704,11 @@ void FixRigid::init()
if (inertia[ibody][2] < EPSILON*max) inertia[ibody][2] = 0.0;
// enforce 3 evectors as a right-handed coordinate system
- // flip 3rd evector if needed
-
- ez0 = ex_space[ibody][1]*ey_space[ibody][2] -
- ex_space[ibody][2]*ey_space[ibody][1];
- ez1 = ex_space[ibody][2]*ey_space[ibody][0] -
- ex_space[ibody][0]*ey_space[ibody][2];
- ez2 = ex_space[ibody][0]*ey_space[ibody][1] -
- ex_space[ibody][1]*ey_space[ibody][0];
-
- if (ez0*ez_space[ibody][0] + ez1*ez_space[ibody][1] +
- ez2*ez_space[ibody][2] < 0.0) {
- ez_space[ibody][0] = -ez_space[ibody][0];
- ez_space[ibody][1] = -ez_space[ibody][1];
- ez_space[ibody][2] = -ez_space[ibody][2];
- }
+ // flip 3rd vector if needed
+
+ MathExtra::cross3(ex_space[ibody],ey_space[ibody],cross);
+ if (MathExtra::dot3(cross,ez_space[ibody]) < 0.0)
+ MathExtra::negate3(ez_space[ibody]);
// create initial quaternion
@@ -823,8 +811,7 @@ void FixRigid::init()
sum[ibody][0] += 0.4 * massone * radius[i]*radius[i];
sum[ibody][1] += 0.4 * massone * radius[i]*radius[i];
sum[ibody][2] += 0.4 * massone * radius[i]*radius[i];
- }
- if (eflags[i] & INERTIA_ELLIPSOID) {
+ } else if (eflags[i] & INERTIA_ELLIPSOID) {
shape = ebonus[ellipsoid[i]].shape;
MathExtra::inertia_ellipsoid(shape,qorient[i],massone,ivec);
sum[ibody][0] += ivec[0];
diff --git a/src/math_extra.cpp b/src/math_extra.cpp
index c4318e8bbe..5160262aff 100644
--- a/src/math_extra.cpp
+++ b/src/math_extra.cpp
@@ -487,7 +487,7 @@ void inertia_line(double length, double theta, double mass, double *inertia)
/* ----------------------------------------------------------------------
compute space-frame inertia tensor of a triangle
v0,v1,v2 = 3 vertices of triangle
- from http://en.wikipedia.org/wiki/Inertia_tensor_of_triangle:
+ from http://en.wikipedia.org/wiki/Inertia_tensor_of_triangle
inertia tensor = a/24 (v0^2 + v1^2 + v2^2 + (v0+v1+v2)^2) I - a Vt S V
a = 2*area of tri = |(v1-v0) x (v2-v0)|
I = 3x3 identity matrix
@@ -523,9 +523,7 @@ void inertia_triangle(double *v0, double *v1, double *v2,
sub3(v2,v0,v2mv0);
cross3(v1mv0,v2mv0,normal);
double a = len3(normal);
- double inv24 = 1.0/24.0;
-
- // NOTE: use mass
+ double inv24 = mass/24.0;
inertia[0] = inv24*a*(sum-vtsv[0][0]);
inertia[1] = inv24*a*(sum-vtsv[1][1]);
diff --git a/src/math_extra.h b/src/math_extra.h
index 3ca98f8f12..44af2e9a8a 100755
--- a/src/math_extra.h
+++ b/src/math_extra.h
@@ -30,6 +30,7 @@ namespace MathExtra {
inline void norm3(double *v);
inline void normalize3(const double *v, double *ans);
inline void snormalize3(const double, const double *v, double *ans);
+ inline void negate3(double *v);
inline void add3(const double *v1, const double *v2, double *ans);
inline void sub3(const double *v1, const double *v2, double *ans);
inline double len3(const double *v);
@@ -156,6 +157,17 @@ void MathExtra::snormalize3(const double length, const double *v, double *ans)
ans[2] = v[2]*scale;
}
+/* ----------------------------------------------------------------------
+ negate vector v
+------------------------------------------------------------------------- */
+
+void MathExtra::negate3(double *v)
+{
+ v[0] = -v[0];
+ v[1] = -v[1];
+ v[2] = -v[2];
+}
+
/* ----------------------------------------------------------------------
ans = v1 + v2
------------------------------------------------------------------------- */
From bce11d1402ab028a34a9b63e10883c56e7be8da1 Mon Sep 17 00:00:00 2001
From: sjplimp
Date: Thu, 28 Apr 2011 18:52:25 +0000
Subject: [PATCH 02/21] git-svn-id:
svn://svn.icms.temple.edu/lammps-ro/trunk@6029
f3b2605a-c512-4ea7-a41b-209d697bcdaa
---
src/compute_cluster_atom.cpp | 2 +-
src/fix_langevin.cpp | 1 -
2 files changed, 1 insertion(+), 2 deletions(-)
diff --git a/src/compute_cluster_atom.cpp b/src/compute_cluster_atom.cpp
index 38a63d078c..0f3f3dd709 100644
--- a/src/compute_cluster_atom.cpp
+++ b/src/compute_cluster_atom.cpp
@@ -104,7 +104,7 @@ void ComputeClusterAtom::compute_peratom()
// grow clusterID array if necessary
- if (atom->nlocal > nmax) {
+ if (atom->nlocal+atom->nghost > nmax) {
memory->destroy(clusterID);
nmax = atom->nmax;
memory->create(clusterID,nmax,"cluster/atom:clusterID");
diff --git a/src/fix_langevin.cpp b/src/fix_langevin.cpp
index 47699bfac7..37da93f90d 100644
--- a/src/fix_langevin.cpp
+++ b/src/fix_langevin.cpp
@@ -338,7 +338,6 @@ void FixLangevin::post_force_no_tally()
}
}
}
-
}
/* ---------------------------------------------------------------------- */
From 995a92b9f3efa35a4a2574069583011d7e2772ca Mon Sep 17 00:00:00 2001
From: sjplimp
Date: Thu, 28 Apr 2011 18:52:32 +0000
Subject: [PATCH 03/21] git-svn-id:
svn://svn.icms.temple.edu/lammps-ro/trunk@6030
f3b2605a-c512-4ea7-a41b-209d697bcdaa
---
src/version.h | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/src/version.h b/src/version.h
index 6970bbeb4a..66a97ccc65 100644
--- a/src/version.h
+++ b/src/version.h
@@ -1 +1 @@
-#define LAMMPS_VERSION "20 Apr 2011"
+#define LAMMPS_VERSION "27 Apr 2011"
From 199c005d935e49b772befeedfc2150f1dbfe7f82 Mon Sep 17 00:00:00 2001
From: sjplimp
Date: Fri, 29 Apr 2011 15:52:26 +0000
Subject: [PATCH 04/21] git-svn-id:
svn://svn.icms.temple.edu/lammps-ro/trunk@6033
f3b2605a-c512-4ea7-a41b-209d697bcdaa
---
src/ASPHERE/compute_temp_asphere.cpp | 227 ++++++++++++++++++---------
src/ASPHERE/compute_temp_asphere.h | 2 +-
src/ASPHERE/fix_nve_asphere.cpp | 8 +-
src/USER-EFF/fix_langevin_eff.cpp | 16 +-
src/compute_temp_sphere.cpp | 143 +++++++++++------
src/compute_temp_sphere.h | 2 +-
src/fix_langevin.cpp | 176 +++++++++++++++++++--
src/fix_langevin.h | 6 +-
src/fix_rigid.cpp | 166 ++++++++++++++++++--
src/fix_rigid.h | 8 +
10 files changed, 601 insertions(+), 153 deletions(-)
diff --git a/src/ASPHERE/compute_temp_asphere.cpp b/src/ASPHERE/compute_temp_asphere.cpp
index e4e1177a7e..f2d34ac72a 100755
--- a/src/ASPHERE/compute_temp_asphere.cpp
+++ b/src/ASPHERE/compute_temp_asphere.cpp
@@ -32,13 +32,16 @@
using namespace LAMMPS_NS;
+enum{ROTATE,ALL};
+
+#define INERTIA 0.2 // moment of inertia for ellipsoid
+
/* ---------------------------------------------------------------------- */
ComputeTempAsphere::ComputeTempAsphere(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg)
{
- if (narg != 3 && narg != 4)
- error->all("Illegal compute temp/asphere command");
+ if (narg < 3) error->all("Illegal compute temp/asphere command");
scalar_flag = vector_flag = 1;
size_vector = 6;
@@ -48,11 +51,24 @@ ComputeTempAsphere::ComputeTempAsphere(LAMMPS *lmp, int narg, char **arg) :
tempbias = 0;
id_bias = NULL;
- if (narg == 4) {
- tempbias = 1;
- int n = strlen(arg[3]) + 1;
- id_bias = new char[n];
- strcpy(id_bias,arg[3]);
+ mode = ALL;
+
+ int iarg = 3;
+ while (iarg < narg) {
+ if (strcmp(arg[iarg],"bias") == 0) {
+ if (iarg+2 > narg) error->all("Illegal compute temp/asphere command");
+ tempbias = 1;
+ int n = strlen(arg[iarg+1]) + 1;
+ id_bias = new char[n];
+ strcpy(id_bias,arg[iarg+1]);
+ iarg += 2;
+ } else if (strcmp(arg[iarg],"dof") == 0) {
+ if (iarg+2 > narg) error->all("Illegal compute temp/asphere command");
+ if (strcmp(arg[iarg+1],"rotate") == 0) mode = ROTATE;
+ else if (strcmp(arg[iarg+1],"all") == 0) mode = ALL;
+ else error->all("Illegal compute temp/asphere command");
+ iarg += 2;
+ } else error->all("Illegal compute temp/asphere command");
}
vector = new double[6];
@@ -76,8 +92,7 @@ ComputeTempAsphere::~ComputeTempAsphere()
void ComputeTempAsphere::init()
{
- // check that all particles are finite-size
- // no point particles allowed, spherical is OK
+ // check that all particles are finite-size, no point particles allowed
int *ellipsoid = atom->ellipsoid;
int *mask = atom->mask;
@@ -114,18 +129,26 @@ void ComputeTempAsphere::init()
void ComputeTempAsphere::dof_compute()
{
// 6 dof for 3d, 3 dof for 2d
+ // which dof are included also depends on mode
// assume full rotation of extended particles
// user should correct this via compute_modify if needed
double natoms = group->count(igroup);
- int nper = 6;
- if (domain->dimension == 2) nper = 3;
+ int nper;
+ if (domain->dimension == 3) {
+ if (mode == ALL) nper = 6;
+ else nper = 3;
+ } else {
+ if (mode == ALL) nper = 3;
+ else nper = 1;
+ }
dof = nper*natoms;
// additional adjustments to dof
- if (tempbias == 1) dof -= tbias->dof_remove(-1) * natoms;
- else if (tempbias == 2) {
+ if (tempbias == 1) {
+ if (mode == ALL) dof -= tbias->dof_remove(-1) * natoms;
+ } else if (tempbias == 2) {
int *mask = atom->mask;
int nlocal = atom->nlocal;
int count = 0;
@@ -154,46 +177,73 @@ double ComputeTempAsphere::compute_scalar()
}
AtomVecEllipsoid::Bonus *bonus = avec->bonus;
- int *ellipsoid = atom->ellipsoid;
double **v = atom->v;
double **angmom = atom->angmom;
double *rmass = atom->rmass;
+ int *ellipsoid = atom->ellipsoid;
int *mask = atom->mask;
int nlocal = atom->nlocal;
double *shape,*quat;
double wbody[3],inertia[3];
double rot[3][3];
- double t = 0.0;
- // sum translationals and rotational energy for each particle
+ // sum translational and rotational energy for each particle
// no point particles since divide by inertia
- for (int i = 0; i < nlocal; i++)
- if (mask[i] & groupbit) {
+ double t = 0.0;
- shape = bonus[ellipsoid[i]].shape;
- quat = bonus[ellipsoid[i]].quat;
+ if (mode == ALL) {
+ for (int i = 0; i < nlocal; i++)
+ if (mask[i] & groupbit) {
+ t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * rmass[i];
- t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * rmass[i];
+ // principal moments of inertia
- // principal moments of inertia
+ shape = bonus[ellipsoid[i]].shape;
+ quat = bonus[ellipsoid[i]].quat;
- inertia[0] = rmass[i] * (shape[1]*shape[1]+shape[2]*shape[2]) / 5.0;
- inertia[1] = rmass[i] * (shape[0]*shape[0]+shape[2]*shape[2]) / 5.0;
- inertia[2] = rmass[i] * (shape[0]*shape[0]+shape[1]*shape[1]) / 5.0;
+ inertia[0] = INERTIA*rmass[i] * (shape[1]*shape[1]+shape[2]*shape[2]);
+ inertia[1] = INERTIA*rmass[i] * (shape[0]*shape[0]+shape[2]*shape[2]);
+ inertia[2] = INERTIA*rmass[i] * (shape[0]*shape[0]+shape[1]*shape[1]);
- // wbody = angular velocity in body frame
+ // wbody = angular velocity in body frame
- MathExtra::quat_to_mat(quat,rot);
- MathExtra::transpose_matvec(rot,angmom[i],wbody);
- wbody[0] /= inertia[0];
- wbody[1] /= inertia[1];
- wbody[2] /= inertia[2];
+ MathExtra::quat_to_mat(quat,rot);
+ MathExtra::transpose_matvec(rot,angmom[i],wbody);
+ wbody[0] /= inertia[0];
+ wbody[1] /= inertia[1];
+ wbody[2] /= inertia[2];
+
+ t += inertia[0]*wbody[0]*wbody[0] +
+ inertia[1]*wbody[1]*wbody[1] + inertia[2]*wbody[2]*wbody[2];
+ }
- t += inertia[0]*wbody[0]*wbody[0] +
- inertia[1]*wbody[1]*wbody[1] + inertia[2]*wbody[2]*wbody[2];
- }
+ } else {
+ for (int i = 0; i < nlocal; i++)
+ if (mask[i] & groupbit) {
+
+ // principal moments of inertia
+
+ shape = bonus[ellipsoid[i]].shape;
+ quat = bonus[ellipsoid[i]].quat;
+
+ inertia[0] = INERTIA*rmass[i] * (shape[1]*shape[1]+shape[2]*shape[2]);
+ inertia[1] = INERTIA*rmass[i] * (shape[0]*shape[0]+shape[2]*shape[2]);
+ inertia[2] = INERTIA*rmass[i] * (shape[0]*shape[0]+shape[1]*shape[1]);
+
+ // wbody = angular velocity in body frame
+
+ MathExtra::quat_to_mat(quat,rot);
+ MathExtra::transpose_matvec(rot,angmom[i],wbody);
+ wbody[0] /= inertia[0];
+ wbody[1] /= inertia[1];
+ wbody[2] /= inertia[2];
+
+ t += inertia[0]*wbody[0]*wbody[0] +
+ inertia[1]*wbody[1]*wbody[1] + inertia[2]*wbody[2]*wbody[2];
+ }
+ }
if (tempbias) tbias->restore_bias_all();
@@ -217,58 +267,93 @@ void ComputeTempAsphere::compute_vector()
}
AtomVecEllipsoid::Bonus *bonus = avec->bonus;
- int *ellipsoid = atom->ellipsoid;
double **v = atom->v;
double **angmom = atom->angmom;
double *rmass = atom->rmass;
+ int *ellipsoid = atom->ellipsoid;
int *mask = atom->mask;
int nlocal = atom->nlocal;
double *shape,*quat;
- double wbody[3],inertia[3];
+ double wbody[3],inertia[3],t[6];
double rot[3][3];
- double massone,t[6];
+ double massone;
+
+ // sum translational and rotational energy for each particle
+ // no point particles since divide by inertia
+
for (i = 0; i < 6; i++) t[i] = 0.0;
- for (i = 0; i < nlocal; i++)
- if (mask[i] & groupbit) {
+ if (mode == ALL) {
+ for (i = 0; i < nlocal; i++)
+ if (mask[i] & groupbit) {
+ massone = rmass[i];
+ t[0] += massone * v[i][0]*v[i][0];
+ t[1] += massone * v[i][1]*v[i][1];
+ t[2] += massone * v[i][2]*v[i][2];
+ t[3] += massone * v[i][0]*v[i][1];
+ t[4] += massone * v[i][0]*v[i][2];
+ t[5] += massone * v[i][1]*v[i][2];
+
+ // principal moments of inertia
- shape = bonus[ellipsoid[i]].shape;
- quat = bonus[ellipsoid[i]].quat;
+ shape = bonus[ellipsoid[i]].shape;
+ quat = bonus[ellipsoid[i]].quat;
- // translational kinetic energy
+ inertia[0] = INERTIA*massone * (shape[1]*shape[1]+shape[2]*shape[2]);
+ inertia[1] = INERTIA*massone * (shape[0]*shape[0]+shape[2]*shape[2]);
+ inertia[2] = INERTIA*massone * (shape[0]*shape[0]+shape[1]*shape[1]);
+
+ // wbody = angular velocity in body frame
+
+ MathExtra::quat_to_mat(quat,rot);
+ MathExtra::transpose_matvec(rot,angmom[i],wbody);
+ wbody[0] /= inertia[0];
+ wbody[1] /= inertia[1];
+ wbody[2] /= inertia[2];
+
+ // rotational kinetic energy
+
+ t[0] += inertia[0]*wbody[0]*wbody[0];
+ t[1] += inertia[1]*wbody[1]*wbody[1];
+ t[2] += inertia[2]*wbody[2]*wbody[2];
+ t[3] += inertia[0]*wbody[0]*wbody[1];
+ t[4] += inertia[1]*wbody[0]*wbody[2];
+ t[5] += inertia[2]*wbody[1]*wbody[2];
+ }
- massone = rmass[i];
- t[0] += massone * v[i][0]*v[i][0];
- t[1] += massone * v[i][1]*v[i][1];
- t[2] += massone * v[i][2]*v[i][2];
- t[3] += massone * v[i][0]*v[i][1];
- t[4] += massone * v[i][0]*v[i][2];
- t[5] += massone * v[i][1]*v[i][2];
+ } else {
+ for (i = 0; i < nlocal; i++)
+ if (mask[i] & groupbit) {
+
+ // principal moments of inertia
- // principal moments of inertia
+ shape = bonus[ellipsoid[i]].shape;
+ quat = bonus[ellipsoid[i]].quat;
+ massone = rmass[i];
- inertia[0] = rmass[i] * (shape[1]*shape[1]+shape[2]*shape[2]) / 5.0;
- inertia[1] = rmass[i] * (shape[0]*shape[0]+shape[2]*shape[2]) / 5.0;
- inertia[2] = rmass[i] * (shape[0]*shape[0]+shape[1]*shape[1]) / 5.0;
-
- // wbody = angular velocity in body frame
-
- MathExtra::quat_to_mat(quat,rot);
- MathExtra::transpose_matvec(rot,angmom[i],wbody);
- wbody[0] /= inertia[0];
- wbody[1] /= inertia[1];
- wbody[2] /= inertia[2];
-
- // rotational kinetic energy
-
- t[0] += inertia[0]*wbody[0]*wbody[0];
- t[1] += inertia[1]*wbody[1]*wbody[1];
- t[2] += inertia[2]*wbody[2]*wbody[2];
- t[3] += inertia[0]*wbody[0]*wbody[1];
- t[4] += inertia[1]*wbody[0]*wbody[2];
- t[5] += inertia[2]*wbody[1]*wbody[2];
- }
+ inertia[0] = INERTIA*massone * (shape[1]*shape[1]+shape[2]*shape[2]);
+ inertia[1] = INERTIA*massone * (shape[0]*shape[0]+shape[2]*shape[2]);
+ inertia[2] = INERTIA*massone * (shape[0]*shape[0]+shape[1]*shape[1]);
+
+ // wbody = angular velocity in body frame
+
+ MathExtra::quat_to_mat(quat,rot);
+ MathExtra::transpose_matvec(rot,angmom[i],wbody);
+ wbody[0] /= inertia[0];
+ wbody[1] /= inertia[1];
+ wbody[2] /= inertia[2];
+
+ // rotational kinetic energy
+
+ t[0] += inertia[0]*wbody[0]*wbody[0];
+ t[1] += inertia[1]*wbody[1]*wbody[1];
+ t[2] += inertia[2]*wbody[2]*wbody[2];
+ t[3] += inertia[0]*wbody[0]*wbody[1];
+ t[4] += inertia[1]*wbody[0]*wbody[2];
+ t[5] += inertia[2]*wbody[1]*wbody[2];
+ }
+ }
if (tempbias) tbias->restore_bias_all();
diff --git a/src/ASPHERE/compute_temp_asphere.h b/src/ASPHERE/compute_temp_asphere.h
index 19e29ebf1b..dde67a1bd5 100755
--- a/src/ASPHERE/compute_temp_asphere.h
+++ b/src/ASPHERE/compute_temp_asphere.h
@@ -36,7 +36,7 @@ class ComputeTempAsphere : public Compute {
void restore_bias(int, double *);
private:
- int fix_dof;
+ int fix_dof,mode;
double tfactor;
char *id_bias;
class Compute *tbias; // ptr to additional bias compute
diff --git a/src/ASPHERE/fix_nve_asphere.cpp b/src/ASPHERE/fix_nve_asphere.cpp
index 4f8c94e208..9e4155581f 100755
--- a/src/ASPHERE/fix_nve_asphere.cpp
+++ b/src/ASPHERE/fix_nve_asphere.cpp
@@ -29,6 +29,8 @@
using namespace LAMMPS_NS;
+#define INERTIA 0.2 // moment of inertia for ellipsoid
+
/* ---------------------------------------------------------------------- */
FixNVEAsphere::FixNVEAsphere(LAMMPS *lmp, int narg, char **arg) :
@@ -103,9 +105,9 @@ void FixNVEAsphere::initial_integrate(int vflag)
shape = bonus[ellipsoid[i]].shape;
quat = bonus[ellipsoid[i]].quat;
- inertia[0] = rmass[i] * (shape[1]*shape[1]+shape[2]*shape[2]) / 5.0;
- inertia[1] = rmass[i] * (shape[0]*shape[0]+shape[2]*shape[2]) / 5.0;
- inertia[2] = rmass[i] * (shape[0]*shape[0]+shape[1]*shape[1]) / 5.0;
+ inertia[0] = INERTIA*rmass[i] * (shape[1]*shape[1]+shape[2]*shape[2]);
+ inertia[1] = INERTIA*rmass[i] * (shape[0]*shape[0]+shape[2]*shape[2]);
+ inertia[2] = INERTIA*rmass[i] * (shape[0]*shape[0]+shape[1]*shape[1]);
// compute omega at 1/2 step from angmom at 1/2 step and current q
// update quaternion a full step via Richardson iteration
diff --git a/src/USER-EFF/fix_langevin_eff.cpp b/src/USER-EFF/fix_langevin_eff.cpp
index 2758ebea07..7328fdb23b 100644
--- a/src/USER-EFF/fix_langevin_eff.cpp
+++ b/src/USER-EFF/fix_langevin_eff.cpp
@@ -88,7 +88,9 @@ void FixLangevinEff::post_force_no_tally()
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);
- if (abs(spin[i])==1) erforce[i] += 0.75*gamma1*ervel[i] + 0.866025404*gamma2*(random->uniform()-0.5);
+ if (abs(spin[i])==1)
+ erforce[i] += 0.75*gamma1*ervel[i] +
+ 0.866025404*gamma2*(random->uniform()-0.5);
}
}
} else if (which == BIAS) {
@@ -105,7 +107,8 @@ void FixLangevinEff::post_force_no_tally()
if (v[i][2] != 0.0)
f[i][2] += gamma1*v[i][2] + gamma2*(random->uniform()-0.5);
if (abs(spin[i])==1 && ervel[i] != 0.0)
- erforce[i] += 0.75*gamma1*ervel[i] + 0.866025404*gamma2*(random->uniform()-0.5);
+ erforce[i] += 0.75*gamma1*ervel[i] +
+ 0.866025404*gamma2*(random->uniform()-0.5);
temperature->restore_bias(i,v[i]);
}
}
@@ -158,7 +161,8 @@ void FixLangevinEff::post_force_tally()
flangevin[i][0] = gamma1*v[i][0] + gamma2*(random->uniform()-0.5);
flangevin[i][1] = gamma1*v[i][1] + gamma2*(random->uniform()-0.5);
flangevin[i][2] = gamma1*v[i][2] + gamma2*(random->uniform()-0.5);
- erforcelangevin[i] = 0.75*gamma1*ervel[i]+0.866025404*gamma2*(random->uniform()-0.5);
+ erforcelangevin[i] = 0.75*gamma1*ervel[i] +
+ 0.866025404*gamma2*(random->uniform()-0.5);
f[i][0] += flangevin[i][0];
f[i][1] += flangevin[i][1];
f[i][2] += flangevin[i][2];
@@ -175,14 +179,16 @@ void FixLangevinEff::post_force_tally()
flangevin[i][0] = gamma1*v[i][0] + gamma2*(random->uniform()-0.5);
flangevin[i][1] = gamma1*v[i][1] + gamma2*(random->uniform()-0.5);
flangevin[i][2] = gamma1*v[i][2] + gamma2*(random->uniform()-0.5);
- erforcelangevin[i] = 0.75*gamma1*ervel[i]+0.866025404*gamma2*(random->uniform()-0.5);
+ erforcelangevin[i] = 0.75*gamma1*ervel[i] +
+ 0.866025404*gamma2*(random->uniform()-0.5);
if (v[i][0] != 0.0) f[i][0] += flangevin[i][0];
else flangevin[i][0] = 0.0;
if (v[i][1] != 0.0) f[i][1] += flangevin[i][1];
else flangevin[i][1] = 0.0;
if (v[i][2] != 0.0) f[i][2] += flangevin[i][2];
else flangevin[i][2] = 0.0;
- if (abs(spin[i])==1 && ervel[i] != 0.0) erforce[i] += erforcelangevin[i];
+ if (abs(spin[i])==1 && ervel[i] != 0.0)
+ erforce[i] += erforcelangevin[i];
temperature->restore_bias(i,v[i]);
}
}
diff --git a/src/compute_temp_sphere.cpp b/src/compute_temp_sphere.cpp
index bad55efdb8..93c9ec74aa 100644
--- a/src/compute_temp_sphere.cpp
+++ b/src/compute_temp_sphere.cpp
@@ -26,6 +26,8 @@
using namespace LAMMPS_NS;
+enum{ROTATE,ALL};
+
#define INERTIA 0.4 // moment of inertia for sphere
/* ---------------------------------------------------------------------- */
@@ -33,8 +35,7 @@ using namespace LAMMPS_NS;
ComputeTempSphere::ComputeTempSphere(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg)
{
- if (narg != 3 && narg != 4)
- error->all("Illegal compute temp/sphere command");
+ if (narg < 3) error->all("Illegal compute temp/sphere command");
scalar_flag = vector_flag = 1;
size_vector = 6;
@@ -44,11 +45,24 @@ ComputeTempSphere::ComputeTempSphere(LAMMPS *lmp, int narg, char **arg) :
tempbias = 0;
id_bias = NULL;
- if (narg == 4) {
- tempbias = 1;
- int n = strlen(arg[3]) + 1;
- id_bias = new char[n];
- strcpy(id_bias,arg[3]);
+ mode = ALL;
+
+ int iarg = 3;
+ while (iarg < narg) {
+ if (strcmp(arg[iarg],"bias") == 0) {
+ if (iarg+2 > narg) error->all("Illegal compute temp/sphere command");
+ tempbias = 1;
+ int n = strlen(arg[iarg+1]) + 1;
+ id_bias = new char[n];
+ strcpy(id_bias,arg[iarg+1]);
+ iarg += 2;
+ } else if (strcmp(arg[iarg],"dof") == 0) {
+ if (iarg+2 > narg) error->all("Illegal compute temp/sphere command");
+ if (strcmp(arg[iarg+1],"rotate") == 0) mode = ROTATE;
+ else if (strcmp(arg[iarg+1],"all") == 0) mode = ALL;
+ else error->all("Illegal compute temp/sphere command");
+ iarg += 2;
+ } else error->all("Illegal compute temp/sphere command");
}
vector = new double[6];
@@ -100,27 +114,34 @@ void ComputeTempSphere::dof_compute()
// 6 or 3 dof for extended/point particles for 3d
// 3 or 2 dof for extended/point particles for 2d
+ // which dof are included also depends on mode
// assume full rotation of extended particles
// user should correct this via compute_modify if needed
- int dimension = domain->dimension;
-
double *radius = atom->radius;
int *mask = atom->mask;
int nlocal = atom->nlocal;
count = 0;
- if (dimension == 3) {
+ if (domain->dimension == 3) {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
- if (radius[i] == 0.0) count += 3;
- else count += 6;
+ if (radius[i] == 0.0) {
+ if (mode == ALL) count += 3;
+ } else {
+ if (mode == ALL) count += 6;
+ else count += 3;
+ }
}
} else {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
- if (radius[i] == 0.0) count += 2;
- else count += 3;
+ if (radius[i] == 0.0) {
+ if (mode == ALL) count += 2;
+ } else {
+ if (mode == ALL) count += 3;
+ else count += 1;
+ }
}
}
@@ -130,28 +151,38 @@ void ComputeTempSphere::dof_compute()
// additional adjustments to dof
if (tempbias == 1) {
- double natoms = group->count(igroup);
- dof -= tbias->dof_remove(-1) * natoms;
+ if (mode == ALL) {
+ double natoms = group->count(igroup);
+ dof -= tbias->dof_remove(-1) * natoms;
+ }
} else if (tempbias == 2) {
int *mask = atom->mask;
int nlocal = atom->nlocal;
count = 0;
- if (dimension == 3) {
+ if (domain->dimension == 3) {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
if (tbias->dof_remove(i)) {
- if (radius[i] == 0.0) count += 3;
- else count += 6;
+ if (radius[i] == 0.0) {
+ if (mode == ALL) count += 3;
+ } else {
+ if (mode == ALL) count += 6;
+ else count += 3;
+ }
}
}
} else {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
if (tbias->dof_remove(i)) {
- if (radius[i] == 0.0) count += 2;
- else count += 3;
+ if (radius[i] == 0.0) {
+ if (mode == ALL) count += 2;
+ } else {
+ if (mode == ALL) count += 3;
+ else count += 1;
+ }
}
}
}
@@ -187,12 +218,19 @@ double ComputeTempSphere::compute_scalar()
double t = 0.0;
- for (int i = 0; i < nlocal; i++)
- if (mask[i] & groupbit) {
- t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * rmass[i];
- t += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] +
- omega[i][2]*omega[i][2]) * INERTIA*radius[i]*radius[i]*rmass[i];
- }
+ if (mode == ALL) {
+ for (int i = 0; i < nlocal; i++)
+ if (mask[i] & groupbit) {
+ t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * rmass[i];
+ t += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] +
+ omega[i][2]*omega[i][2]) * INERTIA*rmass[i]*radius[i]*radius[i];
+ }
+ } else {
+ for (int i = 0; i < nlocal; i++)
+ if (mask[i] & groupbit)
+ t += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] +
+ omega[i][2]*omega[i][2]) * INERTIA*rmass[i]*radius[i]*radius[i];
+ }
if (tempbias) tbias->restore_bias_all();
@@ -225,25 +263,38 @@ void ComputeTempSphere::compute_vector()
double massone,inertiaone,t[6];
for (int i = 0; i < 6; i++) t[i] = 0.0;
- for (int i = 0; i < nlocal; i++)
- if (mask[i] & groupbit) {
- massone = rmass[i];
- t[0] += massone * v[i][0]*v[i][0];
- t[1] += massone * v[i][1]*v[i][1];
- t[2] += massone * v[i][2]*v[i][2];
- t[3] += massone * v[i][0]*v[i][1];
- t[4] += massone * v[i][0]*v[i][2];
- t[5] += massone * v[i][1]*v[i][2];
-
- inertiaone = INERTIA*radius[i]*radius[i]*rmass[i];
- t[0] += inertiaone * omega[i][0]*omega[i][0];
- t[1] += inertiaone * omega[i][1]*omega[i][1];
- t[2] += inertiaone * omega[i][2]*omega[i][2];
- t[3] += inertiaone * omega[i][0]*omega[i][1];
- t[4] += inertiaone * omega[i][0]*omega[i][2];
- t[5] += inertiaone * omega[i][1]*omega[i][2];
- }
-
+ if (mode == ALL) {
+ for (int i = 0; i < nlocal; i++)
+ if (mask[i] & groupbit) {
+ massone = rmass[i];
+ t[0] += massone * v[i][0]*v[i][0];
+ t[1] += massone * v[i][1]*v[i][1];
+ t[2] += massone * v[i][2]*v[i][2];
+ t[3] += massone * v[i][0]*v[i][1];
+ t[4] += massone * v[i][0]*v[i][2];
+ t[5] += massone * v[i][1]*v[i][2];
+
+ inertiaone = INERTIA*rmass[i]*radius[i]*radius[i];
+ t[0] += inertiaone * omega[i][0]*omega[i][0];
+ t[1] += inertiaone * omega[i][1]*omega[i][1];
+ t[2] += inertiaone * omega[i][2]*omega[i][2];
+ t[3] += inertiaone * omega[i][0]*omega[i][1];
+ t[4] += inertiaone * omega[i][0]*omega[i][2];
+ t[5] += inertiaone * omega[i][1]*omega[i][2];
+ }
+ } else {
+ for (int i = 0; i < nlocal; i++)
+ if (mask[i] & groupbit) {
+ inertiaone = INERTIA*rmass[i]*radius[i]*radius[i];
+ t[0] += inertiaone * omega[i][0]*omega[i][0];
+ t[1] += inertiaone * omega[i][1]*omega[i][1];
+ t[2] += inertiaone * omega[i][2]*omega[i][2];
+ t[3] += inertiaone * omega[i][0]*omega[i][1];
+ t[4] += inertiaone * omega[i][0]*omega[i][2];
+ t[5] += inertiaone * omega[i][1]*omega[i][2];
+ }
+ }
+
if (tempbias) tbias->restore_bias_all();
MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world);
diff --git a/src/compute_temp_sphere.h b/src/compute_temp_sphere.h
index 86285061bd..c0b29dce59 100644
--- a/src/compute_temp_sphere.h
+++ b/src/compute_temp_sphere.h
@@ -36,7 +36,7 @@ class ComputeTempSphere : public Compute {
void restore_bias(int, double *);
private:
- int fix_dof;
+ int fix_dof,mode;
double tfactor;
double *inertia;
char *id_bias;
diff --git a/src/fix_langevin.cpp b/src/fix_langevin.cpp
index 37da93f90d..0d233126ef 100644
--- a/src/fix_langevin.cpp
+++ b/src/fix_langevin.cpp
@@ -20,7 +20,9 @@
#include "string.h"
#include "stdlib.h"
#include "fix_langevin.h"
+#include "math_extra.h"
#include "atom.h"
+#include "atom_vec_ellipsoid.h"
#include "force.h"
#include "update.h"
#include "modify.h"
@@ -38,6 +40,9 @@ using namespace LAMMPS_NS;
enum{NOBIAS,BIAS};
+#define SINERTIA 0.4 // moment of inertia for sphere
+#define EINERTIA 0.2 // moment of inertia for ellipsoid
+
/* ---------------------------------------------------------------------- */
FixLangevin::FixLangevin(LAMMPS *lmp, int narg, char **arg) :
@@ -71,6 +76,7 @@ FixLangevin::FixLangevin(LAMMPS *lmp, int narg, char **arg) :
// optional args
for (int i = 1; i <= atom->ntypes; i++) ratio[i] = 1.0;
+ oflag = aflag = 0;
tally = 0;
zeroflag = 0;
@@ -96,9 +102,29 @@ FixLangevin::FixLangevin(LAMMPS *lmp, int narg, char **arg) :
else if (strcmp(arg[iarg+1],"yes") == 0) zeroflag = 1;
else error->all("Illegal fix langevin command");
iarg += 2;
+ } else if (strcmp(arg[iarg],"omega") == 0) {
+ if (iarg+2 > narg) error->all("Illegal fix langevin command");
+ if (strcmp(arg[iarg+1],"no") == 0) oflag = 0;
+ else if (strcmp(arg[iarg+1],"yes") == 0) oflag = 1;
+ else error->all("Illegal fix langevin command");
+ iarg += 2;
+ } else if (strcmp(arg[iarg],"angmom") == 0) {
+ if (iarg+2 > narg) error->all("Illegal fix langevin command");
+ if (strcmp(arg[iarg+1],"no") == 0) aflag = 0;
+ else if (strcmp(arg[iarg+1],"yes") == 0) aflag = 1;
+ else error->all("Illegal fix langevin command");
+ iarg += 2;
} else error->all("Illegal fix langevin command");
}
+ // error check
+
+ if (aflag) {
+ avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
+ if (!avec)
+ error->all("Fix langevin angmom requires atom style ellipsoid");
+ }
+
// set temperature = NULL, user can override via fix_modify if wants bias
id_temp = NULL;
@@ -140,6 +166,35 @@ int FixLangevin::setmask()
void FixLangevin::init()
{
+ if (oflag && !atom->sphere_flag)
+ error->all("Fix langevin omega require atom style sphere");
+ if (aflag && !atom->ellipsoid_flag)
+ error->all("Fix langevin angmom require atom style ellipsoid");
+
+ // if oflag or aflag set, check that all group particles are finite-size
+
+ if (oflag) {
+ double *radius = atom->radius;
+ int *mask = atom->mask;
+ int nlocal = atom->nlocal;
+
+ for (int i = 0; i < nlocal; i++)
+ if (mask[i] & groupbit)
+ if (radius[i] == 0.0)
+ error->one("Fix langevin omega requires extended particles");
+ }
+
+ if (aflag) {
+ int *ellipsoid = atom->ellipsoid;
+ int *mask = atom->mask;
+ int nlocal = atom->nlocal;
+
+ for (int i = 0; i < nlocal; i++)
+ if (mask[i] & groupbit)
+ if (ellipsoid[i] < 0)
+ error->one("Fix langevin angmom requires extended particles");
+ }
+
// set force prefactors
if (!atom->rmass) {
@@ -219,6 +274,11 @@ void FixLangevin::post_force_no_tally()
double fran[3],fsum[3],fsumall[3];
fsum[0] = fsum[1] = fsum[2] = 0.0;
bigint count;
+
+ double boltz = force->boltz;
+ double dt = update->dt;
+ double mvv2e = force->mvv2e;
+ double ftm2v = force->ftm2v;
if (zeroflag) {
count = group->count(igroup);
@@ -227,11 +287,6 @@ void FixLangevin::post_force_no_tally()
}
if (rmass) {
- double boltz = force->boltz;
- double dt = update->dt;
- double mvv2e = force->mvv2e;
- double ftm2v = force->ftm2v;
-
if (which == NOBIAS) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
@@ -280,7 +335,6 @@ void FixLangevin::post_force_no_tally()
} else {
if (which == NOBIAS) {
-
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
gamma1 = gfactor1[type[i]];
@@ -295,7 +349,6 @@ void FixLangevin::post_force_no_tally()
fsum[1] += fran[1];
fsum[2] += fran[2];
}
-
}
} else if (which == BIAS) {
@@ -338,6 +391,11 @@ void FixLangevin::post_force_no_tally()
}
}
}
+
+ // thermostat omega and angmom
+
+ if (oflag) omega_thermostat(tsqrt);
+ if (aflag) angmom_thermostat(tsqrt);
}
/* ---------------------------------------------------------------------- */
@@ -373,12 +431,12 @@ void FixLangevin::post_force_tally()
// test v = 0 since some computes mask non-participating atoms via v = 0
// and added force has extra term not multiplied by v = 0
- if (rmass) {
- double boltz = force->boltz;
- double dt = update->dt;
- double mvv2e = force->mvv2e;
- double ftm2v = force->ftm2v;
+ double boltz = force->boltz;
+ double dt = update->dt;
+ double mvv2e = force->mvv2e;
+ double ftm2v = force->ftm2v;
+ if (rmass) {
if (which == NOBIAS) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
@@ -454,6 +512,100 @@ void FixLangevin::post_force_tally()
}
}
}
+
+ // thermostat omega and angmom
+
+ if (oflag) omega_thermostat(tsqrt);
+ if (aflag) angmom_thermostat(tsqrt);
+}
+
+/* ----------------------------------------------------------------------
+ thermostat rotational dof via omega
+------------------------------------------------------------------------- */
+
+void FixLangevin::omega_thermostat(double tsqrt)
+{
+ double gamma1,gamma2;
+
+ double boltz = force->boltz;
+ double dt = update->dt;
+ double mvv2e = force->mvv2e;
+ double ftm2v = force->ftm2v;
+
+ double **torque = atom->torque;
+ double **omega = atom->omega;
+ double *radius = atom->radius;
+ double *rmass = atom->rmass;
+ int *mask = atom->mask;
+ int *type = atom->type;
+ int nlocal = atom->nlocal;
+
+ double tran[3];
+ double inertiaone;
+
+ for (int i = 0; i < nlocal; i++) {
+ if (mask[i] & groupbit) {
+ inertiaone = SINERTIA*radius[i]*radius[i]*rmass[i];
+ gamma1 = -inertiaone / t_period / ftm2v;
+ gamma2 = sqrt(inertiaone) * sqrt(24.0*boltz/t_period/dt/mvv2e) / ftm2v;
+ gamma1 *= 1.0/ratio[type[i]];
+ gamma2 *= 1.0/sqrt(ratio[type[i]]) * tsqrt;
+ tran[0] = gamma2*(random->uniform()-0.5);
+ tran[1] = gamma2*(random->uniform()-0.5);
+ tran[2] = gamma2*(random->uniform()-0.5);
+ torque[i][0] += gamma1*omega[i][0] + tran[0];
+ torque[i][1] += gamma1*omega[i][1] + tran[1];
+ torque[i][2] += gamma1*omega[i][2] + tran[2];
+ }
+ }
+}
+
+/* ----------------------------------------------------------------------
+ thermostat rotational dof via angmom
+------------------------------------------------------------------------- */
+
+void FixLangevin::angmom_thermostat(double tsqrt)
+{
+ double gamma1,gamma2;
+
+ double boltz = force->boltz;
+ double dt = update->dt;
+ double mvv2e = force->mvv2e;
+ double ftm2v = force->ftm2v;
+
+ AtomVecEllipsoid::Bonus *bonus = avec->bonus;
+ double **torque = atom->torque;
+ double **angmom = atom->angmom;
+ double *rmass = atom->rmass;
+ int *ellipsoid = atom->ellipsoid;
+ int *mask = atom->mask;
+ int *type = atom->type;
+ int nlocal = atom->nlocal;
+
+ double inertia[3],wbody[3],omega[3],tran[3],rot[3][3];
+ double *shape,*quat;
+
+ for (int i = 0; i < nlocal; i++) {
+ if (mask[i] & groupbit) {
+ shape = bonus[ellipsoid[i]].shape;
+ inertia[0] = EINERTIA*rmass[i] * (shape[1]*shape[1]+shape[2]*shape[2]);
+ inertia[1] = EINERTIA*rmass[i] * (shape[0]*shape[0]+shape[2]*shape[2]);
+ inertia[2] = EINERTIA*rmass[i] * (shape[0]*shape[0]+shape[1]*shape[1]);
+ quat = bonus[ellipsoid[i]].quat;
+ MathExtra::mq_to_omega(angmom[i],quat,inertia,omega);
+
+ gamma1 = -1.0 / t_period / ftm2v;
+ gamma2 = sqrt(24.0*boltz/t_period/dt/mvv2e) / ftm2v;
+ gamma1 *= 1.0/ratio[type[i]];
+ gamma2 *= 1.0/sqrt(ratio[type[i]]) * tsqrt;
+ tran[0] = sqrt(inertia[0])*gamma2*(random->uniform()-0.5);
+ tran[1] = sqrt(inertia[1])*gamma2*(random->uniform()-0.5);
+ tran[2] = sqrt(inertia[2])*gamma2*(random->uniform()-0.5);
+ torque[i][0] += inertia[0]*gamma1*omega[0] + tran[0];
+ torque[i][1] += inertia[1]*gamma1*omega[1] + tran[1];
+ torque[i][2] += inertia[2]*gamma1*omega[2] + tran[2];
+ }
+ }
}
/* ----------------------------------------------------------------------
diff --git a/src/fix_langevin.h b/src/fix_langevin.h
index f325befd46..735f6cdcf0 100644
--- a/src/fix_langevin.h
+++ b/src/fix_langevin.h
@@ -41,11 +41,13 @@ class FixLangevin : public Fix {
double memory_usage();
protected:
- int which,tally,zeroflag;
+ int which,tally,zeroflag,oflag,aflag;
double t_start,t_stop,t_period;
double *gfactor1,*gfactor2,*ratio;
double energy,energy_onestep;
+ class AtomVecEllipsoid *avec;
+
int nmax;
double **flangevin;
@@ -57,6 +59,8 @@ class FixLangevin : public Fix {
virtual void post_force_no_tally();
virtual void post_force_tally();
+ void omega_thermostat(double);
+ void angmom_thermostat(double);
};
}
diff --git a/src/fix_rigid.cpp b/src/fix_rigid.cpp
index 4d0686409d..6d43b9949e 100644
--- a/src/fix_rigid.cpp
+++ b/src/fix_rigid.cpp
@@ -25,6 +25,7 @@
#include "modify.h"
#include "group.h"
#include "comm.h"
+#include "random_mars.h"
#include "force.h"
#include "output.h"
#include "memory.h"
@@ -45,11 +46,16 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
{
int i,ibody;
+ scalar_flag = 1;
+ extscalar = 0;
time_integrate = 1;
rigid_flag = 1;
virial_flag = 1;
create_attribute = 1;
+ MPI_Comm_rank(world,&me);
+ MPI_Comm_size(world,&nprocs);
+
// perform initial allocation of atom-based arrays
// register with Atom class
@@ -193,12 +199,14 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
memory->create(imagebody,nbody,"rigid:imagebody");
memory->create(fflag,nbody,3,"rigid:fflag");
memory->create(tflag,nbody,3,"rigid:tflag");
+ memory->create(langextra,nbody,6,"rigid:langextra");
memory->create(sum,nbody,6,"rigid:sum");
memory->create(all,nbody,6,"rigid:all");
memory->create(remapflag,nbody,4,"rigid:remapflag");
// initialize force/torque flags to default = 1.0
+ // for 2d: fz, tx, ty = 0.0
array_flag = 1;
size_array_rows = nbody;
@@ -209,10 +217,13 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
for (i = 0; i < nbody; i++) {
fflag[i][0] = fflag[i][1] = fflag[i][2] = 1.0;
tflag[i][0] = tflag[i][1] = tflag[i][2] = 1.0;
+ if (domain->dimension == 2) fflag[i][2] = tflag[i][0] = tflag[i][1] = 0.0;
}
// parse optional args
+ int seed;
+ langflag = 0;
tempflag = 0;
pressflag = 0;
t_chain = 10;
@@ -238,6 +249,9 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
else if (strcmp(arg[iarg+4],"on") == 0) zflag = 1.0;
else error->all("Illegal fix rigid command");
+ if (domain->dimension == 2 && zflag == 1.0)
+ error->all("Fix rigid z force cannot be on for 2d simulation");
+
int count = 0;
for (int m = mlo; m <= mhi; m++) {
fflag[m-1][0] = xflag;
@@ -266,6 +280,9 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
else if (strcmp(arg[iarg+4],"on") == 0) zflag = 1.0;
else error->all("Illegal fix rigid command");
+ if (domain->dimension == 2 && (xflag == 1.0 || yflag == 1.0))
+ error->all("Fix rigid xy torque cannot be on for 2d simulation");
+
int count = 0;
for (int m = mlo; m <= mhi; m++) {
tflag[m-1][0] = xflag;
@@ -277,10 +294,24 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
iarg += 5;
+ } else if (strcmp(arg[iarg],"langevin") == 0) {
+ if (iarg+5 > narg) error->all("Illegal fix rigid command");
+ if (strcmp(style,"rigid") != 0 && strcmp(style,"rigid/nve") != 0)
+ error->all("Illegal fix rigid command");
+ langflag = 1;
+ t_start = atof(arg[iarg+1]);
+ t_stop = atof(arg[iarg+2]);
+ t_period = atof(arg[iarg+3]);
+ seed = atoi(arg[iarg+4]);
+ if (t_period <= 0.0)
+ error->all("Fix rigid langevin period must be > 0.0");
+ if (seed <= 0) error->all("Illegal fix rigid command");
+ iarg += 5;
+
} else if (strcmp(arg[iarg],"temp") == 0) {
if (iarg+4 > narg) error->all("Illegal fix rigid command");
if (strcmp(style,"rigid/nvt") != 0 && strcmp(style,"rigid/npt") != 0)
- error->all("Illegal fix/rigid command");
+ error->all("Illegal fix rigid command");
tempflag = 1;
t_start = atof(arg[iarg+1]);
t_stop = atof(arg[iarg+2]);
@@ -290,7 +321,7 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
} else if (strcmp(arg[iarg],"press") == 0) {
if (iarg+4 > narg) error->all("Illegal fix rigid command");
if (strcmp(style,"rigid/npt") != 0)
- error->all("Illegal fix/rigid command");
+ error->all("Illegal fix rigid command");
pressflag = 1;
p_start = atof(arg[iarg+1]);
p_stop = atof(arg[iarg+2]);
@@ -300,7 +331,7 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
} else if (strcmp(arg[iarg],"tparam") == 0) {
if (iarg+4 > narg) error->all("Illegal fix rigid command");
if (strcmp(style,"rigid/nvt") != 0)
- error->all("Illegal fix/rigid command");
+ error->all("Illegal fix rigid command");
t_chain = atoi(arg[iarg+1]);
t_iter = atoi(arg[iarg+2]);
t_order = atoi(arg[iarg+3]);
@@ -309,13 +340,18 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
} else if (strcmp(arg[iarg],"pparam") == 0) {
if (iarg+2 > narg) error->all("Illegal fix rigid command");
if (strcmp(style,"rigid/npt") != 0)
- error->all("Illegal fix/rigid command");
+ error->all("Illegal fix rigid command");
p_chain = atoi(arg[iarg+1]);
iarg += 2;
} else error->all("Illegal fix rigid command");
}
+ // initialize Marsaglia RNG with processor-unique seed
+
+ if (langflag) random = new RanMars(lmp,seed + me);
+ else random = NULL;
+
// initialize vector output quantities in case accessed before run
for (i = 0; i < nbody; i++) {
@@ -369,7 +405,7 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
int nsum = 0;
for (ibody = 0; ibody < nbody; ibody++) nsum += nrigid[ibody];
- if (comm->me == 0) {
+ if (me == 0) {
if (screen) fprintf(screen,"%d rigid bodies with %d atoms\n",nbody,nsum);
if (logfile) fprintf(logfile,"%d rigid bodies with %d atoms\n",nbody,nsum);
}
@@ -383,6 +419,8 @@ FixRigid::~FixRigid()
atom->delete_callback(id,0);
+ delete random;
+
// delete locally stored arrays
memory->destroy(body);
@@ -409,6 +447,7 @@ FixRigid::~FixRigid()
memory->destroy(imagebody);
memory->destroy(fflag);
memory->destroy(tflag);
+ memory->destroy(langextra);
memory->destroy(sum);
memory->destroy(all);
@@ -422,6 +461,7 @@ int FixRigid::setmask()
int mask = 0;
mask |= INITIAL_INTEGRATE;
mask |= FINAL_INTEGRATE;
+ if (langflag) mask |= POST_FORCE;
mask |= PRE_NEIGHBOR;
mask |= INITIAL_INTEGRATE_RESPA;
mask |= FINAL_INTEGRATE_RESPA;
@@ -441,7 +481,7 @@ void FixRigid::init()
int count = 0;
for (int i = 0; i < modify->nfix; i++)
if (strcmp(modify->fix[i]->style,"rigid") == 0) count++;
- if (count > 1 && comm->me == 0) error->warning("More than one fix rigid");
+ if (count > 1 && me == 0) error->warning("More than one fix rigid");
// error if npt,nph fix comes before rigid fix
@@ -855,6 +895,15 @@ void FixRigid::init()
fabs(all[ibody][5]/norm) > TOLERANCE)
error->all("Fix rigid: Bad principal moments");
}
+
+ // temperature scale factor
+
+ double ndof = 0.0;
+ for (ibody = 0; ibody < nbody; ibody++) {
+ ndof += fflag[ibody][0] + fflag[ibody][1] + fflag[ibody][2];
+ ndof += tflag[ibody][0] + tflag[ibody][1] + tflag[ibody][2];
+ }
+ tfactor = force->mvv2e / (ndof * force->boltz);
}
/* ---------------------------------------------------------------------- */
@@ -998,6 +1047,13 @@ void FixRigid::setup(int vflag)
torque[ibody][2] = all[ibody][5];
}
+ // zero langextra in case Langevin thermostat not used
+ // no point to calling post_force() here since langextra
+ // is only added to fcm/torque in final_integrate()
+
+ for (ibody = 0; ibody < nbody; ibody++)
+ for (i = 0; i < 6; i++) langextra[ibody][i] = 0.0;
+
// virial setup before call to set_v
if (vflag) v_setup(vflag);
@@ -1072,6 +1128,50 @@ void FixRigid::initial_integrate(int vflag)
set_xv();
}
+/* ----------------------------------------------------------------------
+ apply Langevin thermostat to all 6 DOF of rigid bodies
+ computed by proc 0, broadcast to other procs
+ unlike fix langevin, this stores extra force in extra arrays,
+ which are added in when final_integrate() calculates a new fcm/torque
+------------------------------------------------------------------------- */
+
+void FixRigid::post_force(int vflag)
+{
+ if (me == 0) {
+ double gamma1,gamma2;
+
+ double delta = update->ntimestep - update->beginstep;
+ delta /= update->endstep - update->beginstep;
+ double t_target = t_start + delta * (t_stop-t_start);
+ double tsqrt = sqrt(t_target);
+
+ double boltz = force->boltz;
+ double dt = update->dt;
+ double mvv2e = force->mvv2e;
+ double ftm2v = force->ftm2v;
+
+ for (int i = 0; i < nbody; i++) {
+ gamma1 = -masstotal[i] / t_period / ftm2v;
+ gamma2 = sqrt(masstotal[i]) * tsqrt *
+ sqrt(24.0*boltz/t_period/dt/mvv2e) / ftm2v;
+ langextra[i][0] = gamma1*vcm[i][0] + gamma2*(random->uniform()-0.5);
+ langextra[i][1] = gamma1*vcm[i][1] + gamma2*(random->uniform()-0.5);
+ langextra[i][2] = gamma1*vcm[i][2] + gamma2*(random->uniform()-0.5);
+
+ gamma1 = -1.0 / t_period / ftm2v;
+ gamma2 = tsqrt * sqrt(24.0*boltz/t_period/dt/mvv2e) / ftm2v;
+ langextra[i][3] = inertia[i][0]*gamma1*omega[i][0] +
+ sqrt(inertia[i][0])*gamma2*(random->uniform()-0.5);
+ langextra[i][4] = inertia[i][1]*gamma1*omega[i][1] +
+ sqrt(inertia[i][1])*gamma2*(random->uniform()-0.5);
+ langextra[i][5] = inertia[i][2]*gamma1*omega[i][2] +
+ sqrt(inertia[i][2])*gamma2*(random->uniform()-0.5);
+ }
+ }
+
+ MPI_Bcast(&langextra[0][0],6*nbody,MPI_DOUBLE,0,world);
+}
+
/* ---------------------------------------------------------------------- */
void FixRigid::final_integrate()
@@ -1150,13 +1250,17 @@ void FixRigid::final_integrate()
MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world);
+ // update vcm and angmom
+ // include Langevin thermostat forces
+ // fflag,tflag = 0 for some dimensions in 2d
+
for (ibody = 0; ibody < nbody; ibody++) {
- fcm[ibody][0] = all[ibody][0];
- fcm[ibody][1] = all[ibody][1];
- fcm[ibody][2] = all[ibody][2];
- torque[ibody][0] = all[ibody][3];
- torque[ibody][1] = all[ibody][4];
- torque[ibody][2] = all[ibody][5];
+ fcm[ibody][0] = all[ibody][0] + langextra[ibody][0];
+ fcm[ibody][1] = all[ibody][1] + langextra[ibody][1];
+ fcm[ibody][2] = all[ibody][2] + langextra[ibody][2];
+ torque[ibody][0] = all[ibody][3] + langextra[ibody][3];
+ torque[ibody][1] = all[ibody][4] + langextra[ibody][4];
+ torque[ibody][2] = all[ibody][5] + langextra[ibody][5];
// update vcm by 1/2 step
@@ -1360,7 +1464,7 @@ int FixRigid::dof(int igroup)
if (nall[ibody]+mall[ibody] > 0 &&
nall[ibody]+mall[ibody] != nrigid[ibody]) flag = 1;
}
- if (flag && comm->me == 0)
+ if (flag && me == 0)
error->warning("Computing temperature of portions of rigid bodies");
// remove appropriate DOFs for each rigid body wholly in temperature group
@@ -1834,6 +1938,42 @@ void FixRigid::reset_dt()
dtq = 0.5 * update->dt;
}
+/* ----------------------------------------------------------------------
+ return temperature of collection of rigid bodies
+ non-active DOF are removed by fflag/tflag and in tfactor
+------------------------------------------------------------------------- */
+
+double FixRigid::compute_scalar()
+{
+ double wbody[3],rot[3][3];
+
+ double t = 0.0;
+
+ for (int i = 0; i < nbody; i++) {
+ t += masstotal[i] * (fflag[i][0]*vcm[i][0]*vcm[i][0] +
+ fflag[i][1]*vcm[i][1]*vcm[i][1] + \
+ fflag[i][2]*vcm[i][2]*vcm[i][2]);
+
+ // wbody = angular velocity in body frame
+
+ MathExtra::quat_to_mat(quat[i],rot);
+ MathExtra::transpose_matvec(rot,angmom[i],wbody);
+ if (inertia[i][0] == 0.0) wbody[0] = 0.0;
+ else wbody[0] /= inertia[i][0];
+ if (inertia[i][1] == 0.0) wbody[1] = 0.0;
+ else wbody[1] /= inertia[i][1];
+ if (inertia[i][2] == 0.0) wbody[2] = 0.0;
+ else wbody[2] /= inertia[i][2];
+
+ t += tflag[i][0]*inertia[i][0]*wbody[0]*wbody[0] +
+ tflag[i][1]*inertia[i][1]*wbody[1]*wbody[1] +
+ tflag[i][2]*inertia[i][2]*wbody[2]*wbody[2];
+ }
+
+ t *= tfactor;
+ return t;
+}
+
/* ----------------------------------------------------------------------
return attributes of a rigid body
15 values per body
diff --git a/src/fix_rigid.h b/src/fix_rigid.h
index 3aa343015a..06121ad47a 100644
--- a/src/fix_rigid.h
+++ b/src/fix_rigid.h
@@ -32,9 +32,11 @@ class FixRigid : public Fix {
virtual void init();
virtual void setup(int);
virtual void initial_integrate(int);
+ void post_force(int);
virtual void final_integrate();
void initial_integrate_respa(int, int, int);
void final_integrate_respa(int, int);
+ virtual double compute_scalar();
double memory_usage();
void grow_arrays(int);
@@ -50,6 +52,7 @@ class FixRigid : public Fix {
double compute_array(int, int);
protected:
+ int me,nprocs;
double dtv,dtf,dtq;
double *step_respa;
int triclinic;
@@ -70,6 +73,7 @@ class FixRigid : public Fix {
int *imagebody; // image flags of xcm of each rigid body
double **fflag; // flag for on/off of center-of-mass force
double **tflag; // flag for on/off of center-of-mass torque
+ double **langextra; // Langevin thermostat forces and torques
int *body; // which body each atom is part of (-1 if none)
double **displace; // displacement of each atom in body coords
@@ -85,6 +89,9 @@ class FixRigid : public Fix {
double **qorient; // rotation state of ext particle wrt rigid body
double **dorient; // orientation of dipole mu wrt rigid body
+ double tfactor; // scale factor on temperature of rigid bodies
+ int langflag; // 0/1 = no/yes Langevin thermostat
+
int tempflag; // NVT settings
double t_start,t_stop;
double t_period,t_freq;
@@ -95,6 +102,7 @@ class FixRigid : public Fix {
double p_period,p_freq;
int p_chain;
+ class RanMars *random;
class AtomVecEllipsoid *avec_ellipsoid;
// bitmasks for eflags
From 72644dcb8d26290fcc21da6992876b810be9b737 Mon Sep 17 00:00:00 2001
From: sjplimp
Date: Fri, 29 Apr 2011 16:27:56 +0000
Subject: [PATCH 05/21] git-svn-id:
svn://svn.icms.temple.edu/lammps-ro/trunk@6034
f3b2605a-c512-4ea7-a41b-209d697bcdaa
---
doc/compute_temp_asphere.html | 52 ++++++++++++++++------
doc/compute_temp_asphere.txt | 47 ++++++++++++++------
doc/compute_temp_sphere.html | 56 ++++++++++++++++-------
doc/compute_temp_sphere.txt | 51 ++++++++++++++-------
doc/fix_langevin.html | 38 +++++++++++++---
doc/fix_langevin.txt | 37 ++++++++++++---
doc/fix_rigid.html | 84 +++++++++++++++++++++++++----------
doc/fix_rigid.txt | 84 +++++++++++++++++++++++++----------
8 files changed, 333 insertions(+), 116 deletions(-)
diff --git a/doc/compute_temp_asphere.html b/doc/compute_temp_asphere.html
index 3b29b68e74..daaad528a9 100644
--- a/doc/compute_temp_asphere.html
+++ b/doc/compute_temp_asphere.html
@@ -13,16 +13,29 @@
Syntax:
-compute ID group-ID temp/asphere bias-ID
+compute ID group-ID temp/asphere keyword value ...
-
-If a bias-ID is specified it must be the ID of a temperature compute
-that removes a "bias" velocity from each atom. This allows compute
-temp/sphere to compute its thermal temperature after the translational
-kinetic energy components have been altered in a prescribed way,
-e.g. to remove a velocity profile. Thermostats that use this compute
-will work with this bias term. See the doc pages for individual
-computes that calculate a temperature and the doc pages for fixes that
-perform thermostatting for more details.
-
This compute subtracts out translational degrees-of-freedom due to
fixes that constrain molecular motion, such as fix
shake and fix rigid. This means the
@@ -96,6 +100,26 @@ be altered using the extra option of the
discussion of different ways to compute temperature and perform
thermostatting.
+
+
+The keyword/value option pairs are used in the following ways.
+
+For the bias keyword, bias-ID refers to the ID of a temperature
+compute that removes a "bias" velocity from each atom. This allows
+compute temp/sphere to compute its thermal temperature after the
+translational kinetic energy components have been altered in a
+prescribed way, e.g. to remove a velocity profile. Thermostats that
+use this compute will work with this bias term. See the doc pages for
+individual computes that calculate a temperature and the doc pages for
+fixes that perform thermostatting for more details.
+
+For the dof keyword, a setting of all calculates a temperature
+that includes both translational and rotational degrees of freedom. A
+setting of rotate calculates a temperature that includes only
+rotational degrees of freedom.
+
+
+
Output info:
This compute calculates a global scalar (the temperature) and a global
diff --git a/doc/compute_temp_asphere.txt b/doc/compute_temp_asphere.txt
index b22256fbe8..cdd8870981 100755
--- a/doc/compute_temp_asphere.txt
+++ b/doc/compute_temp_asphere.txt
@@ -10,16 +10,24 @@ compute temp/asphere command :h3
[Syntax:]
-compute ID group-ID temp/asphere bias-ID :pre
+compute ID group-ID temp/asphere keyword value ... :pre
-ID, group-ID are documented in "compute"_compute.html command
-temp/asphere = style name of this compute command
-bias-ID = ID of a temperature compute that removes a velocity bias (optional) :ul
+ID, group-ID are documented in "compute"_compute.html command :ulb,l
+temp/asphere = style name of this compute command :l
+zero or more keyword/value pairs may be appended :l
+keyword = {bias} or {dof} :l
+ {bias} value = bias-ID{uniform} or {gaussian}
+ bias-ID = ID of a temperature compute that removes a velocity bias
+ {dof} value = {all} or {rotate}
+ all = compute temperature of translational and rotational degrees of freedom
+ rotate = compute temperature of just rotational degrees of freedom :pre
+:ule
[Examples:]
compute 1 all temp/asphere
-compute myTemp mobile temp/asphere tempCOM :pre
+compute myTemp mobile temp/asphere bias tempCOM
+compute myTemp mobile temp/asphere dof rotate :pre
[Description:]
@@ -72,15 +80,6 @@ The number of atoms contributing to the temperature is assumed to be
constant for the duration of the run; use the {dynamic} option of the
"compute_modify"_compute_modify.html command if this is not the case.
-If a {bias-ID} is specified it must be the ID of a temperature compute
-that removes a "bias" velocity from each atom. This allows compute
-temp/sphere to compute its thermal temperature after the translational
-kinetic energy components have been altered in a prescribed way,
-e.g. to remove a velocity profile. Thermostats that use this compute
-will work with this bias term. See the doc pages for individual
-computes that calculate a temperature and the doc pages for fixes that
-perform thermostatting for more details.
-
This compute subtracts out translational degrees-of-freedom due to
fixes that constrain molecular motion, such as "fix
shake"_fix_shake.html and "fix rigid"_fix_rigid.html. This means the
@@ -93,6 +92,26 @@ See "this howto section"_Section_howto.html#4_16 of the manual for a
discussion of different ways to compute temperature and perform
thermostatting.
+:line
+
+The keyword/value option pairs are used in the following ways.
+
+For the {bias} keyword, {bias-ID} refers to the ID of a temperature
+compute that removes a "bias" velocity from each atom. This allows
+compute temp/sphere to compute its thermal temperature after the
+translational kinetic energy components have been altered in a
+prescribed way, e.g. to remove a velocity profile. Thermostats that
+use this compute will work with this bias term. See the doc pages for
+individual computes that calculate a temperature and the doc pages for
+fixes that perform thermostatting for more details.
+
+For the {dof} keyword, a setting of {all} calculates a temperature
+that includes both translational and rotational degrees of freedom. A
+setting of {rotate} calculates a temperature that includes only
+rotational degrees of freedom.
+
+:line
+
[Output info:]
This compute calculates a global scalar (the temperature) and a global
diff --git a/doc/compute_temp_sphere.html b/doc/compute_temp_sphere.html
index 31e73a05f5..23e18d16b5 100644
--- a/doc/compute_temp_sphere.html
+++ b/doc/compute_temp_sphere.html
@@ -13,16 +13,29 @@
Syntax:
-compute ID group-ID temp/sphere bias-ID
+compute ID group-ID temp/sphere keyword value ...
-
- ID, group-ID are documented in compute command
-
- temp/sphere = style name of this compute command
-
- bias-ID = ID of a temperature compute that removes a velocity bias (optional)
+
Examples:
compute 1 all temp/sphere
-compute myTemp mobile temp/sphere tempCOM
+compute myTemp mobile temp/sphere bias tempCOM
+compute myTemp mobile temp/sphere dof rotate
Description:
@@ -66,15 +79,6 @@ the vector are ordered xx, yy, zz, xy, xz, yz.
constant for the duration of the run; use the dynamic option of the
compute_modify command if this is not the case.
-If a bias-ID is specified it must be the ID of a temperature compute
-that removes a "bias" velocity from each atom. This allows compute
-temp/sphere to compute its thermal temperature after the translational
-kinetic energy components have been altered in a prescribed way,
-e.g. to remove a velocity profile. Thermostats that use this compute
-will work with this bias term. See the doc pages for individual
-computes that calculate a temperature and the doc pages for fixes that
-perform thermostatting for more details.
-
This compute subtracts out translational degrees-of-freedom due to
fixes that constrain molecular motion, such as fix
shake and fix rigid. This means the
@@ -87,6 +91,26 @@ be altered using the extra option of the
discussion of different ways to compute temperature and perform
thermostatting.
+
+
+The keyword/value option pairs are used in the following ways.
+
+For the bias keyword, bias-ID refers to the ID of a temperature
+compute that removes a "bias" velocity from each atom. This allows
+compute temp/sphere to compute its thermal temperature after the
+translational kinetic energy components have been altered in a
+prescribed way, e.g. to remove a velocity profile. Thermostats that
+use this compute will work with this bias term. See the doc pages for
+individual computes that calculate a temperature and the doc pages for
+fixes that perform thermostatting for more details.
+
+For the dof keyword, a setting of all calculates a temperature
+that includes both translational and rotational degrees of freedom. A
+setting of rotate calculates a temperature that includes only
+rotational degrees of freedom.
+
+
+
Output info:
This compute calculates a global scalar (the temperature) and a global
@@ -116,6 +140,8 @@ particles with radius = 0.0.
compute temp, compute
temp/asphere
-Default: none
+
Default:
+
+The option defaults are no bias and dof = all.