From 199c005d935e49b772befeedfc2150f1dbfe7f82 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Fri, 29 Apr 2011 15:52:26 +0000 Subject: [PATCH] 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