From 838fbc7cdca531169bc90d0d11da3d80558fb5cb Mon Sep 17 00:00:00 2001 From: sjplimp Date: Tue, 29 Apr 2014 20:21:23 +0000 Subject: [PATCH 01/10] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@11835 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/GRANULAR/fix_pour.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/GRANULAR/fix_pour.cpp b/src/GRANULAR/fix_pour.cpp index e135206598..f3fac80314 100644 --- a/src/GRANULAR/fix_pour.cpp +++ b/src/GRANULAR/fix_pour.cpp @@ -209,8 +209,8 @@ FixPour::FixPour(LAMMPS *lmp, int narg, char **arg) : // volume_one = volume of inserted particle (with max possible radius) // in 3d, insure dy >= 1, for quasi-2d simulations - double volume,volume_one=0.0; - dstyle = -1; + double volume,volume_one; + if (domain->dimension == 3) { if (region_style == 1) { double dy = yhi - ylo; From 08ff12cc4986452d813c52b4c4362f27ee802895 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Tue, 29 Apr 2014 21:30:31 +0000 Subject: [PATCH 02/10] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@11836 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/MISC/fix_deposit.cpp | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/MISC/fix_deposit.cpp b/src/MISC/fix_deposit.cpp index 753fa3b888..4f5d204362 100644 --- a/src/MISC/fix_deposit.cpp +++ b/src/MISC/fix_deposit.cpp @@ -466,17 +466,17 @@ void FixDeposit::pre_exchange() for (j = 0; j < nfix; j++) if (fix[j]->create_attribute) fix[j]->set_arrays(n); } - - // FixRigidSmall::set_molecule stores rigid body attributes - // coord is new position of geometric center of mol, not COM - // FixShake::set_molecule stores shake info for molecule - - if (rigidflag) - fixrigid->set_molecule(nlocalprev,maxtag_all,coord,vnew,quat); - else if (shakeflag) - fixshake->set_molecule(nlocalprev,maxtag_all,coord,vnew,quat); } + // FixRigidSmall::set_molecule stores rigid body attributes + // coord is new position of geometric center of mol, not COM + // FixShake::set_molecule stores shake info for molecule + + if (rigidflag) + fixrigid->set_molecule(nlocalprev,maxtag_all,coord,vnew,quat); + else if (shakeflag) + fixshake->set_molecule(nlocalprev,maxtag_all,coord,vnew,quat); + // old code: unsuccessful if no proc performed insertion of an atom // don't think that check is necessary // if get this far, should always be succesful From 9c889c3f34c96f421c4cb83f06479d6fc7f663b8 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Tue, 29 Apr 2014 22:21:39 +0000 Subject: [PATCH 03/10] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@11837 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/RIGID/fix_rigid_nh_small.cpp | 1534 +++++++++++++++++++++++++++++ src/RIGID/fix_rigid_nh_small.h | 195 ++++ src/RIGID/fix_rigid_nph_small.cpp | 92 ++ src/RIGID/fix_rigid_nph_small.h | 53 + src/RIGID/fix_rigid_npt_small.cpp | 104 ++ src/RIGID/fix_rigid_npt_small.h | 65 ++ src/RIGID/fix_rigid_nve_small.cpp | 28 + src/RIGID/fix_rigid_nve_small.h | 36 + src/RIGID/fix_rigid_nvt_small.cpp | 50 + src/RIGID/fix_rigid_nvt_small.h | 60 ++ src/RIGID/fix_rigid_small.cpp | 135 +++ src/RIGID/fix_rigid_small.h | 19 +- 12 files changed, 2370 insertions(+), 1 deletion(-) create mode 100644 src/RIGID/fix_rigid_nh_small.cpp create mode 100644 src/RIGID/fix_rigid_nh_small.h create mode 100644 src/RIGID/fix_rigid_nph_small.cpp create mode 100644 src/RIGID/fix_rigid_nph_small.h create mode 100644 src/RIGID/fix_rigid_npt_small.cpp create mode 100644 src/RIGID/fix_rigid_npt_small.h create mode 100644 src/RIGID/fix_rigid_nve_small.cpp create mode 100644 src/RIGID/fix_rigid_nve_small.h create mode 100644 src/RIGID/fix_rigid_nvt_small.cpp create mode 100644 src/RIGID/fix_rigid_nvt_small.h diff --git a/src/RIGID/fix_rigid_nh_small.cpp b/src/RIGID/fix_rigid_nh_small.cpp new file mode 100644 index 0000000000..29223847c6 --- /dev/null +++ b/src/RIGID/fix_rigid_nh_small.cpp @@ -0,0 +1,1534 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Trung Dac Nguyen (ORNL) + references: Kamberaj et al., J. Chem. Phys. 122, 224114 (2005) + Miller et al., J Chem Phys. 116, 8649-8659 (2002) +------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdio.h" +#include "string.h" +#include "fix_rigid_nh_small.h" +#include "math_extra.h" +#include "atom.h" +#include "compute.h" +#include "domain.h" +#include "update.h" +#include "modify.h" +#include "fix_deform.h" +#include "group.h" +#include "comm.h" +#include "force.h" +#include "kspace.h" +#include "output.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; +using namespace FixConst; + +enum{NONE,XYZ,XY,YZ,XZ}; // same as in FixRigid +enum{ISO,ANISO,TRICLINIC}; // same as in FixRigid + +#define EPSILON 1.0e-7 + +enum{FULL_BODY,INITIAL,FINAL,FORCE_TORQUE,VCM_ANGMOM,XCM_MASS,ITENSOR,DOF}; + +/* ---------------------------------------------------------------------- */ + +FixRigidNHSmall::FixRigidNHSmall(LAMMPS *lmp, int narg, char **arg) : + FixRigidSmall(lmp, narg, arg) +{ + // error checks + + if ((p_flag[0] == 1 && p_period[0] <= 0.0) || + (p_flag[1] == 1 && p_period[1] <= 0.0) || + (p_flag[2] == 1 && p_period[2] <= 0.0)) + error->all(FLERR,"Fix rigid/small npt/nph period must be > 0.0"); + + if (domain->dimension == 2 && p_flag[2]) + error->all(FLERR,"Invalid fix rigid/small npt/nph command for a 2d simulation"); + if (domain->dimension == 2 && (pcouple == YZ || pcouple == XZ)) + error->all(FLERR,"Invalid fix rigid/small npt/nph command for a 2d simulation"); + + if (pcouple == XYZ && (p_flag[0] == 0 || p_flag[1] == 0)) + error->all(FLERR,"Invalid fix rigid/small npt/nph command pressure settings"); + if (pcouple == XYZ && domain->dimension == 3 && p_flag[2] == 0) + error->all(FLERR,"Invalid fix rigid/small npt/nph command pressure settings"); + if (pcouple == XY && (p_flag[0] == 0 || p_flag[1] == 0)) + error->all(FLERR,"Invalid fix rigid/small npt/nph command pressure settings"); + if (pcouple == YZ && (p_flag[1] == 0 || p_flag[2] == 0)) + error->all(FLERR,"Invalid fix rigid/small npt/nph command pressure settings"); + if (pcouple == XZ && (p_flag[0] == 0 || p_flag[2] == 0)) + error->all(FLERR,"Invalid fix rigid/small npt/nph command pressure settings"); + + // require periodicity in tensile dimension + + if (p_flag[0] && domain->xperiodic == 0) + error->all(FLERR, + "Cannot use fix rigid/small npt/nph on a non-periodic dimension"); + if (p_flag[1] && domain->yperiodic == 0) + error->all(FLERR, + "Cannot use fix rigid/small npt/nph on a non-periodic dimension"); + if (p_flag[2] && domain->zperiodic == 0) + error->all(FLERR, + "Cannot use fix rigid/small npt/nph on a non-periodic dimension"); + + if (pcouple == XYZ && domain->dimension == 3 && + (p_start[0] != p_start[1] || p_start[0] != p_start[2] || + p_stop[0] != p_stop[1] || p_stop[0] != p_stop[2] || + p_period[0] != p_period[1] || p_period[0] != p_period[2])) + error->all(FLERR,"Invalid fix rigid/small npt/nph pressure settings"); + if (pcouple == XYZ && domain->dimension == 2 && + (p_start[0] != p_start[1] || p_stop[0] != p_stop[1] || + p_period[0] != p_period[1])) + error->all(FLERR,"Invalid fix rigid/small npt/nph pressure settings"); + if (pcouple == XY && + (p_start[0] != p_start[1] || p_stop[0] != p_stop[1] || + p_period[0] != p_period[1])) + error->all(FLERR,"Invalid fix rigid/small npt/nph pressure settings"); + if (pcouple == YZ && + (p_start[1] != p_start[2] || p_stop[1] != p_stop[2] || + p_period[1] != p_period[2])) + error->all(FLERR,"Invalid fix rigid/small npt/nph pressure settings"); + if (pcouple == XZ && + (p_start[0] != p_start[2] || p_stop[0] != p_stop[2] || + p_period[0] != p_period[2])) + error->all(FLERR,"Invalid fix rigid/small npt/nph pressure settings"); + + if ((tstat_flag && t_period <= 0.0) || + (p_flag[0] && p_period[0] <= 0.0) || + (p_flag[1] && p_period[1] <= 0.0) || + (p_flag[2] && p_period[2] <= 0.0)) + error->all(FLERR,"Fix rigid/small nvt/npt/nph damping parameters must be > 0.0"); + + // memory allocation and initialization + + if (tstat_flag || pstat_flag) { + allocate_chain(); + allocate_order(); + } + + if (tstat_flag) { + eta_t[0] = eta_r[0] = 0.0; + eta_dot_t[0] = eta_dot_r[0] = 0.0; + f_eta_t[0] = f_eta_r[0] = 0.0; + + for (int i = 1; i < t_chain; i++) { + eta_t[i] = eta_r[i] = 0.0; + eta_dot_t[i] = eta_dot_r[i] = 0.0; + } + } + + if (pstat_flag) { + epsilon_dot[0] = epsilon_dot[1] = epsilon_dot[2] = 0.0; + eta_b[0] = eta_dot_b[0] = f_eta_b[0] = 0.0; + for (int i = 1; i < p_chain; i++) + eta_b[i] = eta_dot_b[i] = 0.0; + } + + // rigid body pointers + + nrigidfix = 0; + rfix = NULL; + + vol0 = 0.0; + t0 = 1.0; + + tcomputeflag = 0; + pcomputeflag = 0; +} + +/* ---------------------------------------------------------------------- */ + +FixRigidNHSmall::~FixRigidNHSmall() +{ + if (tstat_flag || pstat_flag) { + deallocate_chain(); + deallocate_order(); + } + + if (rfix) delete [] rfix; + + if (tcomputeflag) { + modify->delete_compute(id_temp); + delete [] id_temp; + } + + // delete pressure if fix created it + + if (pstat_flag) { + if (pcomputeflag) modify->delete_compute(id_press); + delete [] id_press; + } +} + +/* ---------------------------------------------------------------------- */ + +int FixRigidNHSmall::setmask() +{ + int mask = 0; + mask = FixRigidSmall::setmask(); + if (tstat_flag || pstat_flag) mask |= THERMO_ENERGY; + + return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixRigidNHSmall::init() +{ + FixRigidSmall::init(); + + // recheck that dilate group has not been deleted + + if (allremap == 0) { + int idilate = group->find(id_dilate); + if (idilate == -1) + error->all(FLERR,"Fix rigid npt/nph dilate group ID does not exist"); + dilate_group_bit = group->bitmask[idilate]; + } + + // initialize thermostats + // set timesteps, constants + // store Yoshida-Suzuki integrator parameters + + dtv = update->dt; + dtf = 0.5 * update->dt * force->ftm2v; + dtq = 0.5 * update->dt; + + boltz = force->boltz; + nktv2p = force->nktv2p; + mvv2e = force->mvv2e; + dimension = domain->dimension; + + if (force->kspace) kspace_flag = 1; + else kspace_flag = 0; + + // see Table 1 in Kamberaj et al + + if (tstat_flag || pstat_flag) { + if (t_order == 3) { + w[0] = 1.0 / (2.0 - pow(2.0, 1.0/3.0)); + w[1] = 1.0 - 2.0*w[0]; + w[2] = w[0]; + } else if (t_order == 5) { + w[0] = 1.0 / (4.0 - pow(4.0, 1.0/3.0)); + w[1] = w[0]; + w[2] = 1.0 - 4.0 * w[0]; + w[3] = w[0]; + w[4] = w[0]; + } + } + + int icompute; + if (tcomputeflag) { + icompute = modify->find_compute(id_temp); + if (icompute < 0) + error->all(FLERR,"Temp ID for fix rigid npt/nph does not exist"); + temperature = modify->compute[icompute]; + } + + if (pstat_flag) { + if (domain->triclinic) + error->all(FLERR,"fix rigid npt/nph does not yet allow triclinic box"); + + // ensure no conflict with fix deform + + for (int i = 0; i < modify->nfix; i++) + if (strcmp(modify->fix[i]->style,"deform") == 0) { + int *dimflag = ((FixDeform *) modify->fix[i])->dimflag; + if ((p_flag[0] && dimflag[0]) || (p_flag[1] && dimflag[1]) || + (p_flag[2] && dimflag[2])) + error->all(FLERR,"Cannot use fix rigid npt/nph and fix deform on " + "same component of stress tensor"); + } + + // set frequency + + p_freq_max = 0.0; + p_freq_max = MAX(p_freq[0],p_freq[1]); + p_freq_max = MAX(p_freq_max,p_freq[2]); + + // tally the number of dimensions that are barostatted + // set initial volume and reference cell, if not already done + + pdim = p_flag[0] + p_flag[1] + p_flag[2]; + if (vol0 == 0.0) { + if (dimension == 2) vol0 = domain->xprd * domain->yprd; + else vol0 = domain->xprd * domain->yprd * domain->zprd; + } + + // set pressure compute ptr + + icompute = modify->find_compute(id_press); + if (icompute < 0) + error->all(FLERR,"Press ID for fix rigid npt/nph does not exist"); + pressure = modify->compute[icompute]; + + // detect if any rigid fixes exist so rigid bodies move on remap + // rfix[] = indices to each fix rigid + // this will include self + + if (rfix) delete [] rfix; + nrigidfix = 0; + rfix = NULL; + + for (int i = 0; i < modify->nfix; i++) + if (modify->fix[i]->rigid_flag) nrigidfix++; + if (nrigidfix) { + rfix = new int[nrigidfix]; + nrigidfix = 0; + for (int i = 0; i < modify->nfix; i++) + if (modify->fix[i]->rigid_flag) rfix[nrigidfix++] = i; + } + } +} + +/* ---------------------------------------------------------------------- */ + +void FixRigidNHSmall::setup(int vflag) +{ + FixRigidSmall::setup(vflag); + + // total translational and rotational degrees of freedom + + int k,ibody; + double *inertia; + + nf_t = nf_r = dimension * nlocal_body; + for (ibody = 0; ibody < nlocal_body; ibody++) { + inertia = body[ibody].inertia; + for (k = 0; k < domain->dimension; k++) + if (fabs(body[ibody].inertia[k]) < EPSILON) nf_r--; + } + + double nf[2], nfall[2]; + nf[0] = nf_t; + nf[1] = nf_r; + MPI_Allreduce(nf,nfall,2,MPI_DOUBLE,MPI_SUM,world); + nf_t = nfall[0]; + nf_r = nfall[1]; + + g_f = nf_t + nf_r; + onednft = 1.0 + (double)(dimension) / (double)g_f; + onednfr = (double) (dimension) / (double)g_f; + + double mbody[3]; + akin_t = akin_r = 0.0; + for (int ibody = 0; ibody < nlocal_body; ibody++) { + Body *b = &body[ibody]; + MathExtra::transpose_matvec(b->ex_space,b->ey_space,b->ez_space, + b->angmom,mbody); + MathExtra::quatvec(b->quat,mbody,b->conjqm); + b->conjqm[0] *= 2.0; + b->conjqm[1] *= 2.0; + b->conjqm[2] *= 2.0; + b->conjqm[3] *= 2.0; + + if (tstat_flag || pstat_flag) { + akin_t += b->mass*(b->vcm[0]*b->vcm[0] + b->vcm[1]*b->vcm[1] + + b->vcm[2]*b->vcm[2]); + akin_r += b->angmom[0]*b->omega[0] + b->angmom[1]*b->omega[1] + + b->angmom[2]*b->omega[2]; + } + } + + // accumulate translational and rotational kinetic energies + + if (tstat_flag || pstat_flag) { + double ke[2],keall[2]; + ke[0] = akin_t; + ke[1] = akin_r; + MPI_Allreduce(ke,keall,2,MPI_DOUBLE,MPI_SUM,world); + akin_t = keall[0]; + akin_r = keall[1]; + } + + // compute target temperature + + if (tstat_flag) compute_temp_target(); + else if (pstat_flag) { + t0 = temperature->compute_scalar(); + if (t0 == 0.0) { + if (strcmp(update->unit_style,"lj") == 0) t0 = 1.0; + else t0 = 300.0; + } + t_target = t0; + } + + // compute target pressure + // compute current pressure + // trigger virial computation on next timestep + + if (pstat_flag) { + compute_press_target(); + + temperature->compute_scalar(); + if (pstyle == ISO) pressure->compute_scalar(); + else pressure->compute_vector(); + couple(); + pressure->addstep(update->ntimestep+1); + } + + // initialize thermostat/barostat settings + + double kt, t_mass, tb_mass; + kt = boltz * t_target; + + if (tstat_flag) { + t_mass = kt / (t_freq*t_freq); + q_t[0] = nf_t * t_mass; + q_r[0] = nf_r * t_mass; + for (int i = 1; i < t_chain; i++) + q_t[i] = q_r[i] = t_mass; + + for (int i = 1; i < t_chain; i++) { + f_eta_t[i] = (q_t[i-1] * eta_dot_t[i-1] * eta_dot_t[i-1] - kt)/q_t[i]; + f_eta_r[i] = (q_r[i-1] * eta_dot_r[i-1] * eta_dot_r[i-1] - kt)/q_r[i]; + } + } + + // initial forces on barostat thermostat variables + + if (pstat_flag) { + for (int i = 0; i < 3; i++) + if (p_flag[i]) { + epsilon_mass[i] = (g_f + dimension) * kt / (p_freq[i]*p_freq[i]); + epsilon[i] = log(vol0)/dimension; + } + + tb_mass = kt / (p_freq_max * p_freq_max); + q_b[0] = dimension * dimension * tb_mass; + for (int i = 1; i < p_chain; i++) { + q_b[i] = tb_mass; + f_eta_b[i] = (q_b[i] * eta_dot_b[i-1] * eta_dot_b[i-1] - kt)/q_b[i]; + } + } + + // update order/timestep dependent coefficients + + if (tstat_flag || pstat_flag) { + for (int i = 0; i < t_order; i++) { + wdti1[i] = w[i] * dtv / t_iter; + wdti2[i] = wdti1[i] / 2.0; + wdti4[i] = wdti1[i] / 4.0; + } + } +} + +/* ---------------------------------------------------------------------- + perform preforce velocity Verlet integration + see Kamberaj paper for step references +------------------------------------------------------------------------- */ + +void FixRigidNHSmall::initial_integrate(int vflag) +{ + double tmp,scale_r,scale_t[3],scale_v[3]; + double dtfm,mbody[3],tbody[3],fquat[4]; + double dtf2 = dtf * 2.0; + + // compute target temperature + // update thermostat chains coupled to particles + + if (tstat_flag) { + compute_temp_target(); + nhc_temp_integrate(); + } + + // compute target pressure + // update epsilon dot + // update thermostat coupled to barostat + + if (pstat_flag) { + nhc_press_integrate(); + + if (pstyle == ISO) { + temperature->compute_scalar(); + pressure->compute_scalar(); + } else { + temperature->compute_vector(); + pressure->compute_vector(); + } + couple(); + pressure->addstep(update->ntimestep+1); + + compute_press_target(); + nh_epsilon_dot(); + } + + // compute scale variables + + scale_t[0] = scale_t[1] = scale_t[2] = 1.0; + scale_v[0] = scale_v[1] = scale_v[2] = 1.0; + scale_r = 1.0; + + if (tstat_flag) { + tmp = exp(-dtq * eta_dot_t[0]); + scale_t[0] = scale_t[1] = scale_t[2] = tmp; + tmp = exp(-dtq * eta_dot_r[0]); + scale_r = tmp; + } + + if (pstat_flag) { + scale_t[0] *= exp(-dtq * (epsilon_dot[0] + mtk_term2)); + scale_t[1] *= exp(-dtq * (epsilon_dot[1] + mtk_term2)); + scale_t[2] *= exp(-dtq * (epsilon_dot[2] + mtk_term2)); + scale_r *= exp(-dtq * (pdim * mtk_term2)); + + tmp = dtq * epsilon_dot[0]; + scale_v[0] = dtv * exp(tmp) * maclaurin_series(tmp); + tmp = dtq * epsilon_dot[1]; + scale_v[1] = dtv * exp(tmp) * maclaurin_series(tmp); + tmp = dtq * epsilon_dot[2]; + scale_v[2] = dtv * exp(tmp) * maclaurin_series(tmp); + } + + // update xcm, vcm, quat, conjqm and angmom + + for (int ibody = 0; ibody < nlocal_body; ibody++) { + Body *b = &body[ibody]; + + // step 1.1 - update vcm by 1/2 step + + dtfm = dtf / b->mass; + b->vcm[0] += dtfm * b->fcm[0]; + b->vcm[1] += dtfm * b->fcm[1]; + b->vcm[2] += dtfm * b->fcm[2]; + + if (tstat_flag || pstat_flag) { + b->vcm[0] *= scale_t[0]; + b->vcm[1] *= scale_t[1]; + b->vcm[2] *= scale_t[2]; + } + + // step 1.2 - update xcm by full step + + if (!pstat_flag) { + b->xcm[0] += dtv * b->vcm[0]; + b->xcm[1] += dtv * b->vcm[1]; + b->xcm[2] += dtv * b->vcm[2]; + } else { + b->xcm[0] += scale_v[0] * b->vcm[0]; + b->xcm[1] += scale_v[1] * b->vcm[1]; + b->xcm[2] += scale_v[2] * b->vcm[2]; + } + + // step 1.3 - apply torque (body coords) to quaternion momentum + + MathExtra::transpose_matvec(b->ex_space,b->ey_space,b->ez_space, + b->torque,tbody); + MathExtra::quatvec(b->quat,tbody,fquat); + + b->conjqm[0] += dtf2 * fquat[0]; + b->conjqm[1] += dtf2 * fquat[1]; + b->conjqm[2] += dtf2 * fquat[2]; + b->conjqm[3] += dtf2 * fquat[3]; + + if (tstat_flag || pstat_flag) { + b->conjqm[0] *= scale_r; + b->conjqm[1] *= scale_r; + b->conjqm[2] *= scale_r; + b->conjqm[3] *= scale_r; + } + + // step 1.4 to 1.13 - use no_squish rotate to update p and q + + no_squish_rotate(3,b->conjqm,b->quat,b->inertia,dtq); + no_squish_rotate(2,b->conjqm,b->quat,b->inertia,dtq); + no_squish_rotate(1,b->conjqm,b->quat,b->inertia,dtv); + no_squish_rotate(2,b->conjqm,b->quat,b->inertia,dtq); + no_squish_rotate(3,b->conjqm,b->quat,b->inertia,dtq); + + // update exyz_space + // transform p back to angmom + // update angular velocity + + MathExtra::q_to_exyz(b->quat,b->ex_space,b->ey_space, + b->ez_space); + MathExtra::invquatvec(b->quat,b->conjqm,mbody); + MathExtra::matvec(b->ex_space,b->ey_space,b->ez_space, + mbody,b->angmom); + + b->angmom[0] *= 0.5; + b->angmom[1] *= 0.5; + b->angmom[2] *= 0.5; + + MathExtra::angmom_to_omega(b->angmom,b->ex_space,b->ey_space, + b->ez_space,b->inertia,b->omega); + } + + // virial setup before call to set_xv + + if (vflag) v_setup(vflag); + else evflag = 0; + + // forward communicate updated info of all bodies + + commflag = INITIAL; + comm->forward_comm_variable_fix(this); + + // accumulate translational and rotational kinetic energies + + if (tstat_flag || pstat_flag) { + + akin_t = akin_r = 0.0; + for (int ibody = 0; ibody < nlocal_body; ibody++) { + Body *b = &body[ibody]; + akin_t += b->mass*(b->vcm[0]*b->vcm[0] + b->vcm[1]*b->vcm[1] + + b->vcm[2]*b->vcm[2]); + akin_r += b->angmom[0]*b->omega[0] + b->angmom[1]*b->omega[1] + + b->angmom[2]*b->omega[2]; + } + + double ke[2],keall[2]; + ke[0] = akin_t; + ke[1] = akin_r; + MPI_Allreduce(ke,keall,2,MPI_DOUBLE,MPI_SUM,world); + akin_t = keall[0]; + akin_r = keall[1]; + } + + // remap simulation box by 1/2 step + + if (pstat_flag) remap(); + + // set coords/orient and velocity/rotation of atoms in rigid bodies + // from quarternion and omega + + set_xv(); + + // remap simulation box by full step + // redo KSpace coeffs since volume has changed + + if (pstat_flag) { + remap(); + if (kspace_flag) force->kspace->setup(); + } +} + +/* ---------------------------------------------------------------------- */ + +void FixRigidNHSmall::final_integrate() +{ + int i,ibody; + double tmp,scale_t[3],scale_r; + double dtfm,xy,xz,yz; + double mbody[3],tbody[3],fquat[4]; + double dtf2 = dtf * 2.0; + + // compute scale variables + + scale_t[0] = scale_t[1] = scale_t[2] = 1.0; + scale_r = 1.0; + + if (tstat_flag) { + tmp = exp(-1.0 * dtq * eta_dot_t[0]); + scale_t[0] = scale_t[1] = scale_t[2] = tmp; + scale_r = exp(-1.0 * dtq * eta_dot_r[0]); + } + + if (pstat_flag) { + scale_t[0] *= exp(-dtq * (epsilon_dot[0] + mtk_term2)); + scale_t[1] *= exp(-dtq * (epsilon_dot[1] + mtk_term2)); + scale_t[2] *= exp(-dtq * (epsilon_dot[2] + mtk_term2)); + scale_r *= exp(-dtq * (pdim * mtk_term2)); + } + + // sum over atoms to get force and torque on rigid body + + imageint *image = atom->image; + double **x = atom->x; + double **f = atom->f; + int nlocal = atom->nlocal; + + double dx,dy,dz; + double unwrap[3]; + double *xcm,*fcm,*tcm; + + for (ibody = 0; ibody < nlocal_body+nghost_body; ibody++) { + fcm = body[ibody].fcm; + fcm[0] = fcm[1] = fcm[2] = 0.0; + tcm = body[ibody].torque; + tcm[0] = tcm[1] = tcm[2] = 0.0; + } + + for (i = 0; i < nlocal; i++) { + if (atom2body[i] < 0) continue; + Body *b = &body[atom2body[i]]; + + fcm = b->fcm; + fcm[0] += f[i][0]; + fcm[1] += f[i][1]; + fcm[2] += f[i][2]; + + domain->unmap(x[i],image[i],unwrap); + xcm = b->xcm; + dx = unwrap[0] - xcm[0]; + dy = unwrap[1] - xcm[1]; + dz = unwrap[2] - xcm[2]; + + tcm = b->torque; + tcm[0] += dy*f[i][2] - dz*f[i][1]; + tcm[1] += dz*f[i][0] - dx*f[i][2]; + tcm[2] += dx*f[i][1] - dy*f[i][0]; + } + + // extended particles add their torque to torque of body + + if (extended) { + double **torque = atom->torque; + + for (i = 0; i < nlocal; i++) { + if (atom2body[i] < 0) continue; + + if (eflags[i] & TORQUE) { + tcm = body[atom2body[i]].torque; + tcm[0] += torque[i][0]; + tcm[1] += torque[i][1]; + tcm[2] += torque[i][2]; + } + } + } + + // reverse communicate fcm, torque of all bodies + + commflag = FORCE_TORQUE; + comm->reverse_comm_variable_fix(this); + + // include Langevin thermostat forces and torques + + if (langflag) { + for (int ibody = 0; ibody < nlocal_body; ibody++) { + fcm = body[ibody].fcm; + fcm[0] += langextra[ibody][0]; + fcm[1] += langextra[ibody][1]; + fcm[2] += langextra[ibody][2]; + tcm = body[ibody].torque; + tcm[0] += langextra[ibody][3]; + tcm[1] += langextra[ibody][4]; + tcm[2] += langextra[ibody][5]; + } + } + + // update vcm and angmom + // include Langevin thermostat forces + // fflag,tflag = 0 for some dimensions in 2d + + for (ibody = 0; ibody < nbody; ibody++) { + Body *b = &body[ibody]; + + // update vcm by 1/2 step + + dtfm = dtf / b->mass; + if (tstat_flag || pstat_flag) { + b->vcm[0] *= scale_t[0]; + b->vcm[1] *= scale_t[1]; + b->vcm[2] *= scale_t[2]; + } + + b->vcm[0] += dtfm * b->fcm[0]; + b->vcm[1] += dtfm * b->fcm[1]; + b->vcm[2] += dtfm * b->fcm[2]; + + // update conjqm, then transform to angmom, set velocity again + // virial is already setup from initial_integrate + + MathExtra::transpose_matvec(b->ex_space,b->ey_space, + b->ez_space,b->torque,tbody); + MathExtra::quatvec(b->quat,tbody,fquat); + + if (tstat_flag || pstat_flag) { + b->conjqm[0] = scale_r * b->conjqm[0] + dtf2 * fquat[0]; + b->conjqm[1] = scale_r * b->conjqm[1] + dtf2 * fquat[1]; + b->conjqm[2] = scale_r * b->conjqm[2] + dtf2 * fquat[2]; + b->conjqm[3] = scale_r * b->conjqm[3] + dtf2 * fquat[3]; + } else { + b->conjqm[0] += dtf2 * fquat[0]; + b->conjqm[1] += dtf2 * fquat[1]; + b->conjqm[2] += dtf2 * fquat[2]; + b->conjqm[3] += dtf2 * fquat[3]; + } + + MathExtra::invquatvec(b->quat,b->conjqm,mbody); + MathExtra::matvec(b->ex_space,b->ey_space,b->ez_space,mbody,b->angmom); + + b->angmom[0] *= 0.5; + b->angmom[1] *= 0.5; + b->angmom[2] *= 0.5; + + MathExtra::angmom_to_omega(b->angmom,b->ex_space,b->ey_space, + b->ez_space,b->inertia,b->omega); + } + + // forward communicate updated info of all bodies + + commflag = FINAL; + comm->forward_comm_variable_fix(this); + + // accumulate translational and rotational kinetic energies + + if (pstat_flag) { + + akin_t = akin_r = 0.0; + for (int ibody = 0; ibody < nlocal_body; ibody++) { + Body *b = &body[ibody]; + akin_t += b->mass*(b->vcm[0]*b->vcm[0] + b->vcm[1]*b->vcm[1] + + b->vcm[2]*b->vcm[2]); + akin_r += b->angmom[0]*b->omega[0] + b->angmom[1]*b->omega[1] + + b->angmom[2]*b->omega[2]; + } + + double ke[2],keall[2]; + ke[0] = akin_t; + ke[1] = akin_r; + MPI_Allreduce(ke,keall,2,MPI_DOUBLE,MPI_SUM,world); + akin_t = keall[0]; + akin_r = keall[1]; + } + + // set velocity/rotation of atoms in rigid bodies + // virial is already setup from initial_integrate + + set_v(); + + // compute temperature and pressure tensor + // couple to compute current pressure components + // trigger virial computation on next timestep + + if (tcomputeflag) t_current = temperature->compute_scalar(); + if (pstat_flag) { + if (pstyle == ISO) pressure->compute_scalar(); + else pressure->compute_vector(); + couple(); + pressure->addstep(update->ntimestep+1); + } + + if (pstat_flag) nh_epsilon_dot(); + + // update eta_dot_t and eta_dot_r + // update eta_dot_b + + if (tstat_flag) nhc_temp_integrate(); + if (pstat_flag) nhc_press_integrate(); +} + +/* ---------------------------------------------------------------------- */ + +void FixRigidNHSmall::nhc_temp_integrate() +{ + int i,j,k; + double kt,gfkt_t,gfkt_r,tmp,ms,s,s2; + + kt = boltz * t_target; + gfkt_t = nf_t * kt; + gfkt_r = nf_r * kt; + + // update thermostat masses + + double t_mass = boltz * t_target / (t_freq * t_freq); + q_t[0] = nf_t * t_mass; + q_r[0] = nf_r * t_mass; + for (i = 1; i < t_chain; i++) + q_t[i] = q_r[i] = t_mass; + + // update force of thermostats coupled to particles + + f_eta_t[0] = (akin_t * mvv2e - gfkt_t) / q_t[0]; + f_eta_r[0] = (akin_r * mvv2e - gfkt_r) / q_r[0]; + + // multiple timestep iteration + + for (i = 0; i < t_iter; i++) { + for (j = 0; j < t_order; j++) { + + // update thermostat velocities half step + + eta_dot_t[t_chain-1] += wdti2[j] * f_eta_t[t_chain-1]; + eta_dot_r[t_chain-1] += wdti2[j] * f_eta_r[t_chain-1]; + + for (k = 1; k < t_chain; k++) { + tmp = wdti4[j] * eta_dot_t[t_chain-k]; + ms = maclaurin_series(tmp); + s = exp(-1.0 * tmp); + s2 = s * s; + eta_dot_t[t_chain-k-1] = eta_dot_t[t_chain-k-1] * s2 + + wdti2[j] * f_eta_t[t_chain-k-1] * s * ms; + + tmp = wdti4[j] * eta_dot_r[t_chain-k]; + ms = maclaurin_series(tmp); + s = exp(-1.0 * tmp); + s2 = s * s; + eta_dot_r[t_chain-k-1] = eta_dot_r[t_chain-k-1] * s2 + + wdti2[j] * f_eta_r[t_chain-k-1] * s * ms; + } + + // update thermostat positions a full step + + for (k = 0; k < t_chain; k++) { + eta_t[k] += wdti1[j] * eta_dot_t[k]; + eta_r[k] += wdti1[j] * eta_dot_r[k]; + } + + // update thermostat forces + + for (k = 1; k < t_chain; k++) { + f_eta_t[k] = q_t[k-1] * eta_dot_t[k-1] * eta_dot_t[k-1] - kt; + f_eta_t[k] /= q_t[k]; + f_eta_r[k] = q_r[k-1] * eta_dot_r[k-1] * eta_dot_r[k-1] - kt; + f_eta_r[k] /= q_r[k]; + } + + // update thermostat velocities a full step + + for (k = 0; k < t_chain-1; k++) { + tmp = wdti4[j] * eta_dot_t[k+1]; + ms = maclaurin_series(tmp); + s = exp(-1.0 * tmp); + s2 = s * s; + eta_dot_t[k] = eta_dot_t[k] * s2 + wdti2[j] * f_eta_t[k] * s * ms; + tmp = q_t[k] * eta_dot_t[k] * eta_dot_t[k] - kt; + f_eta_t[k+1] = tmp / q_t[k+1]; + + tmp = wdti4[j] * eta_dot_r[k+1]; + ms = maclaurin_series(tmp); + s = exp(-1.0 * tmp); + s2 = s * s; + eta_dot_r[k] = eta_dot_r[k] * s2 + wdti2[j] * f_eta_r[k] * s * ms; + tmp = q_r[k] * eta_dot_r[k] * eta_dot_r[k] - kt; + f_eta_r[k+1] = tmp / q_r[k+1]; + } + + eta_dot_t[t_chain-1] += wdti2[j] * f_eta_t[t_chain-1]; + eta_dot_r[t_chain-1] += wdti2[j] * f_eta_r[t_chain-1]; + } + } +} + +/* ---------------------------------------------------------------------- */ + +void FixRigidNHSmall::nhc_press_integrate() +{ + int i,k; + double tmp,s,s2,ms,kecurrent; + double kt = boltz * t_target; + double lkt_press = kt; + + // update thermostat masses + + double tb_mass = kt / (p_freq_max * p_freq_max); + q_b[0] = tb_mass; + for (int i = 1; i < p_chain; i++) { + q_b[i] = tb_mass; + f_eta_b[i] = q_b[i-1] * eta_dot_b[i-1] * eta_dot_b[i-1] - kt; + f_eta_b[i] /= q_b[i]; + } + + // update forces acting on thermostat + + kecurrent = 0.0; + for (i = 0; i < 3; i++) + if (p_flag[i]) { + epsilon_mass[i] = (g_f + dimension) * kt / (p_freq[i] * p_freq[i]); + kecurrent += epsilon_mass[i] * epsilon_dot[i] * epsilon_dot[i]; + } + + f_eta_b[0] = (kecurrent - lkt_press) / q_b[0]; + + // update thermostat velocities a half step + + eta_dot_b[p_chain-1] += 0.5 * dtq * f_eta_b[p_chain-1]; + + for (k = 0; k < p_chain-1; k++) { + tmp = 0.5 * dtq * eta_dot_b[p_chain-k-1]; + ms = maclaurin_series(tmp); + s = exp(-0.5 * tmp); + s2 = s * s; + eta_dot_b[p_chain-k-2] = eta_dot_b[p_chain-k-2] * s2 + + dtq * f_eta_b[p_chain-k-2] * s * ms; + } + + // update thermostat positions + + for (k = 0; k < p_chain; k++) + eta_b[k] += dtv * eta_dot_b[k]; + + // update epsilon dot + + s = exp(-1.0 * dtq * eta_dot_b[0]); + for (i = 0; i < 3; i++) + if (p_flag[i]) epsilon_dot[i] *= s; + + kecurrent = 0.0; + for (i = 0; i < 3; i++) + if (p_flag[i]) + kecurrent += epsilon_mass[i] * epsilon_dot[i] * epsilon_dot[i]; + + f_eta_b[0] = (kecurrent - lkt_press) / q_b[0]; + + // update thermostat velocites a full step + + for (k = 0; k < p_chain-1; k++) { + tmp = 0.5 * dtq * eta_dot_b[k+1]; + ms = maclaurin_series(tmp); + s = exp(-0.5 * tmp); + s2 = s * s; + eta_dot_b[k] = eta_dot_b[k] * s2 + dtq * f_eta_b[k] * s * ms; + tmp = q_b[k] * eta_dot_b[k] * eta_dot_b[k] - kt; + f_eta_b[k+1] = tmp / q_b[k+1]; + } + + eta_dot_b[p_chain-1] += 0.5 * dtq * f_eta_b[p_chain-1]; + +} + +/* ---------------------------------------------------------------------- + compute kinetic energy in the extended Hamiltonian + conserved quantity = sum of returned energy and potential energy +-----------------------------------------------------------------------*/ + +double FixRigidNHSmall::compute_scalar() +{ + int i,k,ibody; + double kt = boltz * t_target; + double energy,ke_t,ke_q,tmp,Pkq[4]; + + double *vcm,*inertia,*quat; + + // compute the kinetic parts of H_NVE in Kameraj et al (JCP 2005, pp 224114) + + // translational and rotational kinetic energies + + ke_t = 0.0; + ke_q = 0.0; + + for (int i = 0; i < nlocal_body; i++) { + vcm = body[i].vcm; + quat = body[i].quat; + ke_t += body[i].mass * (vcm[0]*vcm[0] + vcm[1]*vcm[1] + + vcm[2]*vcm[2]); + + for (k = 1; k < 4; k++) { + if (k == 1) { + Pkq[0] = -quat[1]; + Pkq[1] = quat[0]; + Pkq[2] = quat[3]; + Pkq[3] = -quat[2]; + } else if (k == 2) { + Pkq[0] = -quat[2]; + Pkq[1] = -quat[3]; + Pkq[2] = quat[0]; + Pkq[3] = quat[1]; + } else if (k == 3) { + Pkq[0] = -quat[3]; + Pkq[1] = quat[2]; + Pkq[2] = -quat[1]; + Pkq[3] = quat[0]; + } + + tmp = body[i].conjqm[0]*Pkq[0] + body[i].conjqm[1]*Pkq[1] + + body[i].conjqm[2]*Pkq[2] + body[i].conjqm[3]*Pkq[3]; + tmp *= tmp; + + if (fabs(body[i].inertia[k-1]) < 1e-6) tmp = 0.0; + else tmp /= (8.0 * body[i].inertia[k-1]); + ke_q += tmp; + } + } + + double ke[2],keall[2]; + ke[0] = ke_t; + ke[1] = ke_q; + MPI_Allreduce(ke,keall,2,MPI_DOUBLE,MPI_SUM,world); + ke_t = keall[0]; + ke_q = keall[1]; + + energy = (ke_t + ke_q) * mvv2e; + + if (tstat_flag) { + + // thermostat chain energy: from equation 12 in Kameraj et al (JCP 2005) + + energy += kt * (nf_t * eta_t[0] + nf_r * eta_r[0]); + + for (i = 1; i < t_chain; i++) + energy += kt * (eta_t[i] + eta_r[i]); + + for (i = 0; i < t_chain; i++) { + energy += 0.5 * q_t[i] * (eta_dot_t[i] * eta_dot_t[i]); + energy += 0.5 * q_r[i] * (eta_dot_r[i] * eta_dot_r[i]); + } + } + + if (pstat_flag) { + + // using equation 22 in Kameraj et al for H_NPT + + for (i = 0; i < 3; i++) + energy += 0.5 * epsilon_mass[i] * epsilon_dot[i] * epsilon_dot[i]; + + double vol; + if (dimension == 2) vol = domain->xprd * domain->yprd; + else vol = domain->xprd * domain->yprd * domain->zprd; + + double p0 = (p_target[0] + p_target[1] + p_target[2]) / 3.0; + energy += p0 * vol / nktv2p; + + for (i = 0; i < p_chain; i++) { + energy += kt * eta_b[i]; + energy += 0.5 * q_b[i] * (eta_dot_b[i] * eta_dot_b[i]); + } + } + + return energy; +} + +/* ---------------------------------------------------------------------- */ + +void FixRigidNHSmall::couple() +{ + double *tensor = pressure->vector; + + if (pstyle == ISO) { + p_current[0] = p_current[1] = p_current[2] = pressure->scalar; + } else if (pcouple == XYZ) { + double ave = 1.0/3.0 * (tensor[0] + tensor[1] + tensor[2]); + p_current[0] = p_current[1] = p_current[2] = ave; + } else if (pcouple == XY) { + double ave = 0.5 * (tensor[0] + tensor[1]); + p_current[0] = p_current[1] = ave; + p_current[2] = tensor[2]; + } else if (pcouple == YZ) { + double ave = 0.5 * (tensor[1] + tensor[2]); + p_current[1] = p_current[2] = ave; + p_current[0] = tensor[0]; + } else if (pcouple == XZ) { + double ave = 0.5 * (tensor[0] + tensor[2]); + p_current[0] = p_current[2] = ave; + p_current[1] = tensor[1]; + } else { + p_current[0] = tensor[0]; + p_current[1] = tensor[1]; + p_current[2] = tensor[2]; + } +} + +/* ---------------------------------------------------------------------- */ + +void FixRigidNHSmall::remap() +{ + int i; + double oldlo,oldhi,ctr,expfac; + + double **x = atom->x; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + // epsilon is not used, except for book-keeping + + for (i = 0; i < 3; i++) epsilon[i] += dtq * epsilon_dot[i]; + + // convert pertinent atoms and rigid bodies to lamda coords + + if (allremap) domain->x2lamda(nlocal); + else { + for (i = 0; i < nlocal; i++) + if (mask[i] & dilate_group_bit) + domain->x2lamda(x[i],x[i]); + } + + if (nrigidfix) + for (i = 0; i < nrigidfix; i++) + modify->fix[rfix[i]]->deform(0); + + // reset global and local box to new size/shape + + for (i = 0; i < 3; i++) { + if (p_flag[i]) { + oldlo = domain->boxlo[i]; + oldhi = domain->boxhi[i]; + ctr = 0.5 * (oldlo + oldhi); + expfac = exp(dtq * epsilon_dot[i]); + domain->boxlo[i] = (oldlo-ctr)*expfac + ctr; + domain->boxhi[i] = (oldhi-ctr)*expfac + ctr; + } + } + + domain->set_global_box(); + domain->set_local_box(); + + // convert pertinent atoms and rigid bodies back to box coords + + if (allremap) domain->lamda2x(nlocal); + else { + for (i = 0; i < nlocal; i++) + if (mask[i] & dilate_group_bit) + domain->lamda2x(x[i],x[i]); + } + + if (nrigidfix) + for (i = 0; i< nrigidfix; i++) + modify->fix[rfix[i]]->deform(1); +} + +/* ---------------------------------------------------------------------- + compute target temperature and kinetic energy +-----------------------------------------------------------------------*/ + +void FixRigidNHSmall::compute_temp_target() +{ + double delta = update->ntimestep - update->beginstep; + if (delta != 0.0) delta /= update->endstep - update->beginstep; + + t_target = t_start + delta * (t_stop-t_start); +} + +/* ---------------------------------------------------------------------- + compute hydrostatic target pressure +-----------------------------------------------------------------------*/ + +void FixRigidNHSmall::compute_press_target() +{ + double delta = update->ntimestep - update->beginstep; + if (delta != 0.0) delta /= update->endstep - update->beginstep; + + p_hydro = 0.0; + for (int i = 0; i < 3; i++) + if (p_flag[i]) { + p_target[i] = p_start[i] + delta * (p_stop[i]-p_start[i]); + p_hydro += p_target[i]; + } + p_hydro /= pdim; +} + +/* ---------------------------------------------------------------------- + apply evolution operators to quat, quat momentum + see Miller paper cited in fix rigid/nvt and fix rigid/npt +------------------------------------------------------------------------- */ + +void FixRigidNHSmall::no_squish_rotate(int k, double *p, double *q, + double *inertia, double dt) +{ + double phi,c_phi,s_phi,kp[4],kq[4]; + + // apply permuation operator on p and q, get kp and kq + + if (k == 1) { + kq[0] = -q[1]; kp[0] = -p[1]; + kq[1] = q[0]; kp[1] = p[0]; + kq[2] = q[3]; kp[2] = p[3]; + kq[3] = -q[2]; kp[3] = -p[2]; + } else if (k == 2) { + kq[0] = -q[2]; kp[0] = -p[2]; + kq[1] = -q[3]; kp[1] = -p[3]; + kq[2] = q[0]; kp[2] = p[0]; + kq[3] = q[1]; kp[3] = p[1]; + } else if (k == 3) { + kq[0] = -q[3]; kp[0] = -p[3]; + kq[1] = q[2]; kp[1] = p[2]; + kq[2] = -q[1]; kp[2] = -p[1]; + kq[3] = q[0]; kp[3] = p[0]; + } + + // obtain phi, cosines and sines + + phi = p[0]*kq[0] + p[1]*kq[1] + p[2]*kq[2] + p[3]*kq[3]; + if (fabs(inertia[k-1]) < 1e-6) phi *= 0.0; + else phi /= 4.0 * inertia[k-1]; + c_phi = cos(dt * phi); + s_phi = sin(dt * phi); + + // advance p and q + + p[0] = c_phi*p[0] + s_phi*kp[0]; + p[1] = c_phi*p[1] + s_phi*kp[1]; + p[2] = c_phi*p[2] + s_phi*kp[2]; + p[3] = c_phi*p[3] + s_phi*kp[3]; + + q[0] = c_phi*q[0] + s_phi*kq[0]; + q[1] = c_phi*q[1] + s_phi*kq[1]; + q[2] = c_phi*q[2] + s_phi*kq[2]; + q[3] = c_phi*q[3] + s_phi*kq[3]; +} + +/* ---------------------------------------------------------------------- + update epsilon_dot +-----------------------------------------------------------------------*/ + +void FixRigidNHSmall::nh_epsilon_dot() +{ + int i; + double volume,scale,f_epsilon; + + if (dimension == 2) volume = domain->xprd*domain->yprd; + else volume = domain->xprd*domain->yprd*domain->zprd; + + // MTK terms + + mtk_term1 = (akin_t + akin_r) * mvv2e / g_f; + + scale = exp(-1.0 * dtq * eta_dot_b[0]); + + for (i = 0; i < 3; i++) + if (p_flag[i]) { + f_epsilon = (p_current[i]-p_hydro)*volume / nktv2p + mtk_term1; + f_epsilon /= epsilon_mass[i]; + epsilon_dot[i] += dtq * f_epsilon; + epsilon_dot[i] *= scale; + } + + mtk_term2 = 0.0; + for (i = 0; i < 3; i++) + if (p_flag[i]) mtk_term2 += epsilon_dot[i]; + mtk_term2 /= g_f; +} + +/* ---------------------------------------------------------------------- + pack entire state of Fix into one write +------------------------------------------------------------------------- */ + +void FixRigidNHSmall::write_restart(FILE *fp) +{ + if (tstat_flag == 0 && pstat_flag == 0) return; + + int nsize = 2; // tstat_flag and pstat_flag + + if (tstat_flag) { + nsize += 1; // t_chain + nsize += 4*t_chain; // eta_t, eta_r, eta_dot_t, eta_dot_r + } + + if (pstat_flag) { + nsize += 7; // p_chain, epsilon(3) and epsilon_dot(3) + nsize += 2*p_chain; + } + + double *list; + memory->create(list,nsize,"rigid_nh:list"); + + int n = 0; + + list[n++] = tstat_flag; + if (tstat_flag) { + list[n++] = t_chain; + for (int i = 0; i < t_chain; i++) { + list[n++] = eta_t[i]; + list[n++] = eta_r[i]; + list[n++] = eta_dot_t[i]; + list[n++] = eta_dot_r[i]; + } + } + + list[n++] = pstat_flag; + if (pstat_flag) { + list[n++] = epsilon[0]; + list[n++] = epsilon[1]; + list[n++] = epsilon[2]; + list[n++] = epsilon_dot[0]; + list[n++] = epsilon_dot[1]; + list[n++] = epsilon_dot[2]; + + list[n++] = p_chain; + for (int i = 0; i < p_chain; i++) { + list[n++] = eta_b[i]; + list[n++] = eta_dot_b[i]; + } + } + + if (comm->me == 0) { + int size = (nsize)*sizeof(double); + fwrite(&size,sizeof(int),1,fp); + fwrite(list,sizeof(double),nsize,fp); + } + + memory->destroy(list); +} + +/* ---------------------------------------------------------------------- + use state info from restart file to restart the Fix +------------------------------------------------------------------------- */ + +void FixRigidNHSmall::restart(char *buf) +{ + int n = 0; + double *list = (double *) buf; + int flag = static_cast (list[n++]); + + if (flag) { + int m = static_cast (list[n++]); + if (tstat_flag && m == t_chain) { + for (int i = 0; i < t_chain; i++) { + eta_t[i] = list[n++]; + eta_r[i] = list[n++]; + eta_dot_t[i] = list[n++]; + eta_dot_r[i] = list[n++]; + } + } else n += 4*m; + } + + flag = static_cast (list[n++]); + if (flag) { + epsilon[0] = list[n++]; + epsilon[1] = list[n++]; + epsilon[2] = list[n++]; + epsilon_dot[0] = list[n++]; + epsilon_dot[1] = list[n++]; + epsilon_dot[2] = list[n++]; + + int m = static_cast (list[n++]); + if (pstat_flag && m == p_chain) { + for (int i = 0; i < p_chain; i++) { + eta_b[i] = list[n++]; + eta_dot_b[i] = list[n++]; + } + } else n += 2*m; + } +} + +/* ---------------------------------------------------------------------- */ + +int FixRigidNHSmall::modify_param(int narg, char **arg) +{ + if (strcmp(arg[0],"temp") == 0) { + if (narg < 2) error->all(FLERR,"Illegal fix_modify command"); + if (!pstat_flag) error->all(FLERR,"Illegal fix_modify command"); + if (tcomputeflag) { + modify->delete_compute(id_temp); + tcomputeflag = 0; + } + delete [] id_temp; + int n = strlen(arg[1]) + 1; + id_temp = new char[n]; + strcpy(id_temp,arg[1]); + + int icompute = modify->find_compute(arg[1]); + if (icompute < 0) + error->all(FLERR,"Could not find fix_modify temperature ID"); + temperature = modify->compute[icompute]; + + if (temperature->tempflag == 0) + error->all(FLERR, + "Fix_modify temperature ID does not compute temperature"); + if (temperature->igroup != 0 && comm->me == 0) + error->warning(FLERR,"Temperature for fix modify is not for group all"); + + // reset id_temp of pressure to new temperature ID + + if (pstat_flag) { + icompute = modify->find_compute(id_press); + if (icompute < 0) + error->all(FLERR,"Pressure ID for fix modify does not exist"); + modify->compute[icompute]->reset_extra_compute_fix(id_temp); + } + + return 2; + + } else if (strcmp(arg[0],"press") == 0) { + if (narg < 2) error->all(FLERR,"Illegal fix_modify command"); + if (!pstat_flag) error->all(FLERR,"Illegal fix_modify command"); + if (pcomputeflag) { + modify->delete_compute(id_press); + pcomputeflag = 0; + } + delete [] id_press; + int n = strlen(arg[1]) + 1; + id_press = new char[n]; + strcpy(id_press,arg[1]); + + int icompute = modify->find_compute(arg[1]); + if (icompute < 0) error->all(FLERR,"Could not find fix_modify pressure ID"); + pressure = modify->compute[icompute]; + + if (pressure->pressflag == 0) + error->all(FLERR,"Fix_modify pressure ID does not compute pressure"); + return 2; + } + + return 0; +} + +/* ---------------------------------------------------------------------- */ + +void FixRigidNHSmall::allocate_chain() +{ + if (tstat_flag) { + q_t = new double[t_chain]; + q_r = new double[t_chain]; + eta_t = new double[t_chain]; + eta_r = new double[t_chain]; + eta_dot_t = new double[t_chain]; + eta_dot_r = new double[t_chain]; + f_eta_t = new double[t_chain]; + f_eta_r = new double[t_chain]; + } + + if (pstat_flag) { + q_b = new double[p_chain]; + eta_b = new double[p_chain]; + eta_dot_b = new double[p_chain]; + f_eta_b = new double[p_chain]; + } +} + +/* ---------------------------------------------------------------------- */ + +void FixRigidNHSmall::reset_target(double t_new) +{ + t_start = t_stop = t_new; +} + +/* ---------------------------------------------------------------------- */ + +void FixRigidNHSmall::allocate_order() +{ + w = new double[t_order]; + wdti1 = new double[t_order]; + wdti2 = new double[t_order]; + wdti4 = new double[t_order]; +} + +/* ---------------------------------------------------------------------- */ + +void FixRigidNHSmall::deallocate_chain() +{ + if (tstat_flag) { + delete [] q_t; + delete [] q_r; + delete [] eta_t; + delete [] eta_r; + delete [] eta_dot_t; + delete [] eta_dot_r; + delete [] f_eta_t; + delete [] f_eta_r; + } + + if (pstat_flag) { + delete [] q_b; + delete [] eta_b; + delete [] eta_dot_b; + delete [] f_eta_b; + } +} + +/* ---------------------------------------------------------------------- */ + +void FixRigidNHSmall::deallocate_order() +{ + delete [] w; + delete [] wdti1; + delete [] wdti2; + delete [] wdti4; +} + diff --git a/src/RIGID/fix_rigid_nh_small.h b/src/RIGID/fix_rigid_nh_small.h new file mode 100644 index 0000000000..3dd8d80994 --- /dev/null +++ b/src/RIGID/fix_rigid_nh_small.h @@ -0,0 +1,195 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifndef LMP_FIX_RIGID_NH_SMALL_H +#define LMP_FIX_RIGID_NH_SMALL_H + +#include "fix_rigid_small.h" + +namespace LAMMPS_NS { + +class FixRigidNHSmall : public FixRigidSmall { + public: + FixRigidNHSmall(class LAMMPS *, int, char **); + virtual ~FixRigidNHSmall(); + virtual int setmask(); + virtual void init(); + virtual void setup(int); + virtual void initial_integrate(int); + virtual void final_integrate(); + virtual double compute_scalar(); + int modify_param(int, char **); + void write_restart(FILE *); + void restart(char *buf); + void reset_target(double); + + protected: + double boltz,nktv2p,mvv2e; // boltzman constant, conversion factors + + int dimension; // # of dimensions + int nf_t,nf_r; // trans/rot degrees of freedom + double onednft,onednfr; // factors 1 + dimension/trans(rot) degrees of freedom + double *w,*wdti1,*wdti2,*wdti4; // Yoshida-Suzuki coefficients + double *q_t,*q_r; // trans/rot thermostat masses + double *eta_t,*eta_r; // trans/rot thermostat positions + double *eta_dot_t,*eta_dot_r; // trans/rot thermostat velocities + double *f_eta_t,*f_eta_r; // trans/rot thermostat forces + + double epsilon_mass[3], *q_b; // baro/thermo masses + double epsilon[3],*eta_b; // baro/thermo positions + double epsilon_dot[3],*eta_dot_b; // baro/thermo velocities + double *f_eta_b; // thermo forces + double akin_t,akin_r; // translational/rotational kinetic energies + + int kspace_flag; // 1 if KSpace invoked, 0 if not + int nrigidfix; // number of rigid fixes + int *rfix; // indicies of rigid fixes + + double vol0; // reference volume + double t0; // reference temperature + int pdim,g_f; // number of barostatted dims, total DoFs + double p_hydro; // hydrostatic target pressure + double p_freq_max; // maximum barostat frequency + + double mtk_term1,mtk_term2; // Martyna-Tobias-Klein corrections + + double t_target,t_current; + double t_freq; + + char *id_temp,*id_press; + class Compute *temperature,*pressure; + int tcomputeflag,pcomputeflag; + + void couple(); + void remap(); + void nhc_temp_integrate(); + void nhc_press_integrate(); + + virtual void compute_temp_target(); + void compute_press_target(); + void nh_epsilon_dot(); + + void allocate_chain(); + void allocate_order(); + void deallocate_chain(); + void deallocate_order(); + + void no_squish_rotate(int, double *, double *, double *, double); + inline double maclaurin_series(double); +}; + +inline double FixRigidNHSmall::maclaurin_series(double x) +{ + double x2,x4; + x2 = x * x; + x4 = x2 * x2; + return (1.0 + (1.0/6.0) * x2 + (1.0/120.0) * x4 + (1.0/5040.0) * x2 * x4 + + (1.0/362880.0) * x4 * x4); +} + +} + +#endif + +/* ERROR/WARNING messages: + +E: Fix rigid npt/nph period must be > 0.0 + +Self-explanatory. + +E: Invalid fix rigid npt/nph command for a 2d simulation + +Cannot control z dimension in a 2d model. + +E: Invalid fix rigid npt/nph command pressure settings + +If multiple dimensions are coupled, those dimensions must be +specified. + +E: Cannot use fix rigid npt/nph on a non-periodic dimension + +When specifying a diagonal pressure component, the dimension must be +periodic. + +E: Invalid fix rigid npt/nph pressure settings + +Settings for coupled dimensions must be the same. + +E: Fix rigid nvt/npt/nph damping parameters must be > 0.0 + +Self-explanatory. + +E: Fix rigid npt/nph dilate group ID does not exist + +Self-explanatory. + +E: Temp ID for fix rigid npt/nph does not exist + +Specified compute temperature must be valid. + +E: fix rigid npt/nph does not yet allow triclinic box + +Self-explanatory. + +E: Cannot use fix rigid npt/nph and fix deform on same component of stress tensor + +This would be changing the same box dimension twice. + +E: Press ID for fix rigid npt/nph does not exist + +Specified compute pressure must be valid. + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +E: Could not find fix_modify temperature ID + +The compute ID for computing temperature does not exist. + +E: Fix_modify temperature ID does not compute temperature + +The compute ID assigned to the fix must compute temperature. + +W: Temperature for fix modify is not for group all + +The temperature compute is being used with a pressure calculation +which does operate on group all, so this may be inconsistent. + +E: Pressure ID for fix modify does not exist + +Self-explanatory. + +E: Could not find fix_modify pressure ID + +The compute ID for computing pressure does not exist. + +E: Fix_modify pressure ID does not compute pressure + +The compute ID assigned to the fix must compute pressure. + +U: Target temperature for fix rigid nvt/npt cannot be 0.0 + +Self-explanatory. + +U: Temperature ID for fix rigid npt/nph does not exist + +Self-explanatory. + +U: Pressure ID for fix rigid npt/nph does not exist + +Self-explanatory. + +*/ diff --git a/src/RIGID/fix_rigid_nph_small.cpp b/src/RIGID/fix_rigid_nph_small.cpp new file mode 100644 index 0000000000..714360dc61 --- /dev/null +++ b/src/RIGID/fix_rigid_nph_small.cpp @@ -0,0 +1,92 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Trung Dac Nguyen (ORNL) + references: Kamberaj et al., J. Chem. Phys. 122, 224114 (2005) + Miller et al., J Chem Phys. 116, 8649-8659 (2002) +------------------------------------------------------------------------- */ + +#include "string.h" +#include "fix_rigid_nph_small.h" +#include "domain.h" +#include "modify.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +FixRigidNPHSmall::FixRigidNPHSmall(LAMMPS *lmp, int narg, char **arg) : + FixRigidNHSmall(lmp, narg, arg) +{ + // other setting are made by parent + + scalar_flag = 1; + restart_global = 1; + box_change_size = 1; + extscalar = 1; + + // error checks + + if (pstat_flag == 0) + error->all(FLERR,"Pressure control must be used with fix nph"); + if (tstat_flag == 1) + error->all(FLERR,"Temperature control must not be used with fix nph"); + if (p_start[0] < 0.0 || p_start[1] < 0.0 || p_start[2] < 0.0 || + p_stop[0] < 0.0 || p_stop[1] < 0.0 || p_stop[2] < 0.0) + error->all(FLERR,"Target pressure for fix rigid/nph cannot be 0.0"); + + // convert input periods to frequency + p_freq[0] = p_freq[1] = p_freq[2] = 0.0; + + if (p_flag[0]) p_freq[0] = 1.0 / p_period[0]; + if (p_flag[1]) p_freq[1] = 1.0 / p_period[1]; + if (p_flag[2]) p_freq[2] = 1.0 / p_period[2]; + + // create a new compute temp style + // id = fix-ID + temp + // compute group = all since pressure is always global (group all) + // and thus its KE/temperature contribution should use group all + + int n = strlen(id) + 6; + id_temp = new char[n]; + strcpy(id_temp,id); + strcat(id_temp,"_temp"); + + char **newarg = new char*[3]; + newarg[0] = id_temp; + newarg[1] = (char *) "all"; + newarg[2] = (char *) "temp"; + modify->add_compute(3,newarg); + delete [] newarg; + tcomputeflag = 1; + + // create a new compute pressure style + // id = fix-ID + press, compute group = all + // pass id_temp as 4th arg to pressure constructor + + n = strlen(id) + 7; + id_press = new char[n]; + strcpy(id_press,id); + strcat(id_press,"_press"); + + newarg = new char*[4]; + newarg[0] = id_press; + newarg[1] = (char *) "all"; + newarg[2] = (char *) "pressure"; + newarg[3] = id_temp; + modify->add_compute(4,newarg); + delete [] newarg; + pcomputeflag = 1; +} diff --git a/src/RIGID/fix_rigid_nph_small.h b/src/RIGID/fix_rigid_nph_small.h new file mode 100644 index 0000000000..ff772f8e4d --- /dev/null +++ b/src/RIGID/fix_rigid_nph_small.h @@ -0,0 +1,53 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef FIX_CLASS + +FixStyle(rigid/nph/small,FixRigidNPHSmall) + +#else + +#ifndef LMP_FIX_RIGID_NPH_SMALL_H +#define LMP_FIX_RIGID_NPH_SMALL_H + +#include "fix_rigid_nh_small.h" + +namespace LAMMPS_NS { + +class FixRigidNPHSmall : public FixRigidNHSmall { + public: + FixRigidNPHSmall(class LAMMPS *, int, char **); + ~FixRigidNPHSmall() {} +}; + + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Pressure control must be used with fix rigid nph/small + +UNDOCUMENTED + +E: Temperature control must not be used with fix rigid/nph/small + +UNDOCUMENTED + +E: Target pressure for fix rigid/nph/small cannot be 0.0 + +UNDOCUMENTED + +*/ diff --git a/src/RIGID/fix_rigid_npt_small.cpp b/src/RIGID/fix_rigid_npt_small.cpp new file mode 100644 index 0000000000..c627e568ba --- /dev/null +++ b/src/RIGID/fix_rigid_npt_small.cpp @@ -0,0 +1,104 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Trung Dac Nguyen (ORNL) + references: Kamberaj et al., J. Chem. Phys. 122, 224114 (2005) + Miller et al., J Chem Phys. 116, 8649-8659 (2002) +------------------------------------------------------------------------- */ + +#include "string.h" +#include "fix_rigid_npt_small.h" +#include "domain.h" +#include "modify.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +FixRigidNPTSmall::FixRigidNPTSmall(LAMMPS *lmp, int narg, char **arg) : + FixRigidNHSmall(lmp, narg, arg) +{ + // other setting are made by parent + + scalar_flag = 1; + restart_global = 1; + box_change_size = 1; + extscalar = 1; + + // error checks + + if (tstat_flag == 0 || pstat_flag == 0) + error->all(FLERR,"Did not set temp or press for fix rigid/npt/small"); + if (t_start <= 0.0 || t_stop <= 0.0) + error->all(FLERR,"Target temperature for fix rigid/npt/small cannot be 0.0"); + if (p_start[0] < 0.0 || p_start[1] < 0.0 || p_start[2] < 0.0 || + p_stop[0] < 0.0 || p_stop[1] < 0.0 || p_stop[2] < 0.0) + error->all(FLERR,"Target pressure for fix rigid/npt/small cannot be 0.0"); + + if (t_period <= 0.0) error->all(FLERR,"Fix rigid/npt/small period must be > 0.0"); + + // thermostat chain parameters + + if (t_chain < 1) error->all(FLERR,"Illegal fix_modify command"); + if (t_iter < 1) error->all(FLERR,"Illegal fix_modify command"); + if (t_order != 3 && t_order != 5) + error->all(FLERR,"Fix_modify order must be 3 or 5"); + + // convert input periods to frequency + + t_freq = 0.0; + p_freq[0] = p_freq[1] = p_freq[2] = 0.0; + + t_freq = 1.0 / t_period; + if (p_flag[0]) p_freq[0] = 1.0 / p_period[0]; + if (p_flag[1]) p_freq[1] = 1.0 / p_period[1]; + if (p_flag[2]) p_freq[2] = 1.0 / p_period[2]; + + // create a new compute temp style + // id = fix-ID + temp + // compute group = all since pressure is always global (group all) + // and thus its KE/temperature contribution should use group all + + int n = strlen(id) + 6; + id_temp = new char[n]; + strcpy(id_temp,id); + strcat(id_temp,"_temp"); + + char **newarg = new char*[3]; + newarg[0] = id_temp; + newarg[1] = (char *) "all"; + newarg[2] = (char *) "temp"; + modify->add_compute(3,newarg); + delete [] newarg; + tcomputeflag = 1; + + // create a new compute pressure style + // id = fix-ID + press, compute group = all + // pass id_temp as 4th arg to pressure constructor + + n = strlen(id) + 7; + id_press = new char[n]; + strcpy(id_press,id); + strcat(id_press,"_press"); + + newarg = new char*[4]; + newarg[0] = id_press; + newarg[1] = (char *) "all"; + newarg[2] = (char *) "pressure"; + newarg[3] = id_temp; + modify->add_compute(4,newarg); + delete [] newarg; + pcomputeflag = 1; +} diff --git a/src/RIGID/fix_rigid_npt_small.h b/src/RIGID/fix_rigid_npt_small.h new file mode 100644 index 0000000000..fa697beb44 --- /dev/null +++ b/src/RIGID/fix_rigid_npt_small.h @@ -0,0 +1,65 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef FIX_CLASS + +FixStyle(rigid/npt/small,FixRigidNPTSmall) + +#else + +#ifndef LMP_FIX_RIGID_NPT_SMALL_H +#define LMP_FIX_RIGID_NPT_SMALL_H + +#include "fix_rigid_nh_small.h" + +namespace LAMMPS_NS { + +class FixRigidNPTSmall : public FixRigidNHSmall { + public: + FixRigidNPTSmall(class LAMMPS *, int, char **); + ~FixRigidNPTSmall() {} +}; + + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Did not set temp or press for fix rigid/npt/small + +UNDOCUMENTED + +E: Target temperature for fix rigid/npt/small cannot be 0.0 + +UNDOCUMENTED + +E: Target pressure for fix rigid/npt/small cannot be 0.0 + +UNDOCUMENTED + +E: Fix rigid/npt/small period must be > 0.0 + +UNDOCUMENTED + +E: Illegal ... command + +UNDOCUMENTED + +E: Fix_modify order must be 3 or 5 + +UNDOCUMENTED + +*/ diff --git a/src/RIGID/fix_rigid_nve_small.cpp b/src/RIGID/fix_rigid_nve_small.cpp new file mode 100644 index 0000000000..fdea05db09 --- /dev/null +++ b/src/RIGID/fix_rigid_nve_small.cpp @@ -0,0 +1,28 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Trung Dac Nguyen (ORNL) + references: Kamberaj et al., J. Chem. Phys. 122, 224114 (2005) + Miller et al., J Chem Phys. 116, 8649-8659 (2002) +------------------------------------------------------------------------- */ + +#include "fix_rigid_nve_small.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +FixRigidNVESmall::FixRigidNVESmall(LAMMPS *lmp, int narg, char **arg) : + FixRigidNHSmall(lmp, narg, arg) {} + diff --git a/src/RIGID/fix_rigid_nve_small.h b/src/RIGID/fix_rigid_nve_small.h new file mode 100644 index 0000000000..13d510f6fd --- /dev/null +++ b/src/RIGID/fix_rigid_nve_small.h @@ -0,0 +1,36 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef FIX_CLASS + +FixStyle(rigid/nve/small,FixRigidNVESmall) + +#else + +#ifndef LMP_FIX_RIGID_NVE_SMALL_H +#define LMP_FIX_RIGID_NVE_SMALL_H + +#include "fix_rigid_nh_small.h" + +namespace LAMMPS_NS { + +class FixRigidNVESmall : public FixRigidNHSmall { + public: + FixRigidNVESmall(class LAMMPS *, int, char **); + ~FixRigidNVESmall() {} +}; + +} + +#endif +#endif diff --git a/src/RIGID/fix_rigid_nvt_small.cpp b/src/RIGID/fix_rigid_nvt_small.cpp new file mode 100644 index 0000000000..f21f7168ac --- /dev/null +++ b/src/RIGID/fix_rigid_nvt_small.cpp @@ -0,0 +1,50 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Trung Dac Nguyen (ORNL) + references: Kamberaj et al., J. Chem. Phys. 122, 224114 (2005) + Miller et al., J Chem Phys. 116, 8649-8659 (2002) +------------------------------------------------------------------------- */ + +#include "fix_rigid_nvt_small.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +FixRigidNVTSmall::FixRigidNVTSmall(LAMMPS *lmp, int narg, char **arg) : + FixRigidNHSmall(lmp, narg, arg) +{ + // other settings are made by parent + + scalar_flag = 1; + restart_global = 1; + extscalar = 1; + + // error checking + // convert input period to frequency + + if (tstat_flag == 0) + error->all(FLERR,"Did not set temp for fix rigid/nvt/small"); + if (t_start < 0.0 || t_stop <= 0.0) + error->all(FLERR,"Target temperature for fix rigid/nvt/small cannot be 0.0"); + if (t_period <= 0.0) error->all(FLERR,"Fix rigid/nvt/small period must be > 0.0"); + t_freq = 1.0 / t_period; + + if (t_chain < 1) error->all(FLERR,"Fix rigid nvt/small t_chain should not be less than 1"); + if (t_iter < 1) error->all(FLERR,"Fix rigid nvt/small t_iter should not be less than 1"); + if (t_order != 3 && t_order != 5) + error->all(FLERR,"Fix rigid nvt/small t_order must be 3 or 5"); +} diff --git a/src/RIGID/fix_rigid_nvt_small.h b/src/RIGID/fix_rigid_nvt_small.h new file mode 100644 index 0000000000..39120a3d21 --- /dev/null +++ b/src/RIGID/fix_rigid_nvt_small.h @@ -0,0 +1,60 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef FIX_CLASS + +FixStyle(rigid/nvt/small,FixRigidNVTSmall) + +#else + +#ifndef LMP_FIX_RIGID_NVT_SMALL_H +#define LMP_FIX_RIGID_NVT_SMALL_H + +#include "fix_rigid_nh_small.h" + +namespace LAMMPS_NS { + +class FixRigidNVTSmall : public FixRigidNHSmall { + public: + FixRigidNVTSmall(class LAMMPS *, int, char **); + ~FixRigidNVTSmall() {} +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Did not set temp for fix rigid/nvt/small + +UNDOCUMENTED + +E: Target temperature for fix rigid/nvt/small cannot be 0.0 + +UNDOCUMENTED + +E: Fix rigid/nvt/small period must be > 0.0 + +UNDOCUMENTED + +E: Illegal ... command + +UNDOCUMENTED + +E: Fix_modify order must be 3 or 5 + +UNDOCUMENTED + +*/ diff --git a/src/RIGID/fix_rigid_small.cpp b/src/RIGID/fix_rigid_small.cpp index d82f17d0c0..285cd6e710 100644 --- a/src/RIGID/fix_rigid_small.cpp +++ b/src/RIGID/fix_rigid_small.cpp @@ -59,6 +59,9 @@ FixRigidSmall *FixRigidSmall::frsptr; #define DELTA_BODY 10000 +enum{NONE,XYZ,XY,YZ,XZ}; // same as in FixRigid +enum{ISO,ANISO,TRICLINIC}; // same as in FixRigid + enum{FULL_BODY,INITIAL,FINAL,FORCE_TORQUE,VCM_ANGMOM,XCM_MASS,ITENSOR,DOF}; /* ---------------------------------------------------------------------- */ @@ -125,6 +128,23 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) : infile = NULL; onemol = NULL; + tstat_flag = 0; + pstat_flag = 0; + allremap = 1; + id_dilate = NULL; + t_chain = 10; + t_iter = 1; + t_order = 3; + p_chain = 10; + + pcouple = NONE; + pstyle = ANISO; + + for (int i = 0; i < 3; i++) { + p_start[i] = p_stop[i] = p_period[i] = 0.0; + p_flag[i] = 0; + } + int iarg = 4; while (iarg < narg) { if (strcmp(arg[iarg],"langevin") == 0) { @@ -159,6 +179,121 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) : "fix rigid/small has multiple molecules"); onemol = atom->molecules[imol]; iarg += 2; + + } else if (strcmp(arg[iarg],"temp") == 0) { + if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid/small command"); + if (strcmp(style,"rigid/nvt/small") != 0 && + strcmp(style,"rigid/npt/small") != 0) + error->all(FLERR,"Illegal fix rigid command"); + tstat_flag = 1; + t_start = force->numeric(FLERR,arg[iarg+1]); + t_stop = force->numeric(FLERR,arg[iarg+2]); + t_period = force->numeric(FLERR,arg[iarg+3]); + iarg += 4; + + } else if (strcmp(arg[iarg],"iso") == 0) { + if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid/small command"); + if (strcmp(style,"rigid/npt/small") != 0 && + strcmp(style,"rigid/nph/small") != 0) + error->all(FLERR,"Illegal fix rigid/small command"); + pcouple = XYZ; + p_start[0] = p_start[1] = p_start[2] = force->numeric(FLERR,arg[iarg+1]); + p_stop[0] = p_stop[1] = p_stop[2] = force->numeric(FLERR,arg[iarg+2]); + p_period[0] = p_period[1] = p_period[2] = + force->numeric(FLERR,arg[iarg+3]); + p_flag[0] = p_flag[1] = p_flag[2] = 1; + if (domain->dimension == 2) { + p_start[2] = p_stop[2] = p_period[2] = 0.0; + p_flag[2] = 0; + } + iarg += 4; + + } else if (strcmp(arg[iarg],"aniso") == 0) { + if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid/small command"); + if (strcmp(style,"rigid/npt/small") != 0 && + strcmp(style,"rigid/nph/small") != 0) + error->all(FLERR,"Illegal fix rigid/small command"); + p_start[0] = p_start[1] = p_start[2] = force->numeric(FLERR,arg[iarg+1]); + p_stop[0] = p_stop[1] = p_stop[2] = force->numeric(FLERR,arg[iarg+2]); + p_period[0] = p_period[1] = p_period[2] = + force->numeric(FLERR,arg[iarg+3]); + p_flag[0] = p_flag[1] = p_flag[2] = 1; + if (domain->dimension == 2) { + p_start[2] = p_stop[2] = p_period[2] = 0.0; + p_flag[2] = 0; + } + iarg += 4; + + } else if (strcmp(arg[iarg],"x") == 0) { + if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid/small command"); + p_start[0] = force->numeric(FLERR,arg[iarg+1]); + p_stop[0] = force->numeric(FLERR,arg[iarg+2]); + p_period[0] = force->numeric(FLERR,arg[iarg+3]); + p_flag[0] = 1; + iarg += 4; + + } else if (strcmp(arg[iarg],"y") == 0) { + if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid/small command"); + p_start[1] = force->numeric(FLERR,arg[iarg+1]); + p_stop[1] = force->numeric(FLERR,arg[iarg+2]); + p_period[1] = force->numeric(FLERR,arg[iarg+3]); + p_flag[1] = 1; + iarg += 4; + + } else if (strcmp(arg[iarg],"z") == 0) { + if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid/small command"); + p_start[2] = force->numeric(FLERR,arg[iarg+1]); + p_stop[2] = force->numeric(FLERR,arg[iarg+2]); + p_period[2] = force->numeric(FLERR,arg[iarg+3]); + p_flag[2] = 1; + iarg += 4; + + } else if (strcmp(arg[iarg],"couple") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix rigid/small command"); + if (strcmp(arg[iarg+1],"xyz") == 0) pcouple = XYZ; + else if (strcmp(arg[iarg+1],"xy") == 0) pcouple = XY; + else if (strcmp(arg[iarg+1],"yz") == 0) pcouple = YZ; + else if (strcmp(arg[iarg+1],"xz") == 0) pcouple = XZ; + else if (strcmp(arg[iarg+1],"none") == 0) pcouple = NONE; + else error->all(FLERR,"Illegal fix rigid/small command"); + iarg += 2; + + } else if (strcmp(arg[iarg],"dilate") == 0) { + if (iarg+2 > narg) + error->all(FLERR,"Illegal fix rigid/small nvt/npt/nph command"); + if (strcmp(arg[iarg+1],"all") == 0) allremap = 1; + else { + allremap = 0; + delete [] id_dilate; + int n = strlen(arg[iarg+1]) + 1; + id_dilate = new char[n]; + strcpy(id_dilate,arg[iarg+1]); + int idilate = group->find(id_dilate); + if (idilate == -1) + error->all(FLERR,"Fix rigid/small nvt/npt/nph dilate group ID " + "does not exist"); + } + iarg += 2; + + } else if (strcmp(arg[iarg],"tparam") == 0) { + if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid/small command"); + if (strcmp(style,"rigid/nvt/small") != 0 && + strcmp(style,"rigid/npt/small") != 0) + error->all(FLERR,"Illegal fix rigid/small command"); + t_chain = force->numeric(FLERR,arg[iarg+1]); + t_iter = force->numeric(FLERR,arg[iarg+2]); + t_order = force->numeric(FLERR,arg[iarg+3]); + iarg += 4; + + } else if (strcmp(arg[iarg],"pchain") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix rigid/small command"); + if (strcmp(style,"rigid/npt/small") != 0 && + strcmp(style,"rigid/nph/small") != 0) + error->all(FLERR,"Illegal fix rigid/small command"); + p_chain = force->numeric(FLERR,arg[iarg+1]); + iarg += 2; + + } else error->all(FLERR,"Illegal fix rigid/small command"); } diff --git a/src/RIGID/fix_rigid_small.h b/src/RIGID/fix_rigid_small.h index e3df283a59..81cdf463ba 100644 --- a/src/RIGID/fix_rigid_small.h +++ b/src/RIGID/fix_rigid_small.h @@ -97,7 +97,8 @@ class FixRigidSmall : public Fix { double ez_space[3]; double angmom[3]; // space-frame angular momentum of body double omega[3]; // space-frame omega of body - imageint image; // image flags of xcm + double conjqm[4]; // conjugate quaternion momentum + imageint image; // image flags of xcm int remapflag[4]; // PBC remap flags int ilocal; // index of owning atom }; @@ -151,6 +152,22 @@ class FixRigidSmall : public Fix { int maxlang; // max size of langextra class RanMars *random; // RNG + int tstat_flag,pstat_flag; // 0/1 = no/yes thermostat/barostat + + int t_chain,t_iter,t_order; + + double p_start[3],p_stop[3]; + double p_period[3],p_freq[3]; + int p_flag[3]; + int pcouple,pstyle; + int p_chain; + + int allremap; // remap all atoms + int dilate_group_bit; // mask for dilation group + char *id_dilate; // group name to dilate + + double p_current[3],p_target[3]; + // molecules added on-the-fly as rigid bodies class Molecule *onemol; From 3f09726907c952f98a814dcc35f74123d1f58d79 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Tue, 29 Apr 2014 22:23:31 +0000 Subject: [PATCH 04/10] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@11838 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- doc/Section_commands.html | 8 +++--- doc/Section_commands.txt | 4 +++ doc/compute_voronoi_atom.html | 18 ++++++------- doc/fix_rigid.html | 49 ++++++++++++++++++++--------------- doc/fix_rigid.txt | 45 +++++++++++++++++--------------- 5 files changed, 69 insertions(+), 55 deletions(-) diff --git a/doc/Section_commands.html b/doc/Section_commands.html index 38c83ae038..2bec1769e9 100644 --- a/doc/Section_commands.html +++ b/doc/Section_commands.html @@ -402,10 +402,10 @@ of each style or click on the style itself for a full description: nve/spherenve/trinvtnvt/aspherenvt/sllodnvt/sphereonewayorient/fcc planeforcepoemspourpress/berendsenprintproperty/atomqeq/combreax/bonds recenterrestrainrigidrigid/nphrigid/nptrigid/nverigid/nvtrigid/small -setforceshakespringspring/rgspring/selfsrdstore/forcestore/state -temp/berendsentemp/csvrtemp/rescalethermal/conductivitytmdttmtune/kspacevector -viscosityviscouswall/colloidwall/granwall/harmonicwall/lj1043wall/lj126wall/lj93 -wall/pistonwall/reflectwall/regionwall/srd +rigid/small/nphrigid/small/nptrigid/small/nverigid/small/nvtsetforceshakespringspring/rg +spring/selfsrdstore/forcestore/statetemp/berendsentemp/csvrtemp/rescalethermal/conductivity +tmdttmtune/kspacevectorviscosityviscouswall/colloidwall/gran +wall/harmonicwall/lj1043wall/lj126wall/lj93wall/pistonwall/reflectwall/regionwall/srd

These are fix styles contributed by users, which can be used if diff --git a/doc/Section_commands.txt b/doc/Section_commands.txt index f236c19a8d..e6acf39b18 100644 --- a/doc/Section_commands.txt +++ b/doc/Section_commands.txt @@ -531,6 +531,10 @@ of each style or click on the style itself for a full description: "rigid/nve"_fix_rigid.html, "rigid/nvt"_fix_rigid.html, "rigid/small"_fix_rigid.html, +"rigid/small/nph"_fix_rigid.html, +"rigid/small/npt"_fix_rigid.html, +"rigid/small/nve"_fix_rigid.html, +"rigid/small/nvt"_fix_rigid.html, "setforce"_fix_setforce.html, "shake"_fix_shake.html, "spring"_fix_spring.html, diff --git a/doc/compute_voronoi_atom.html b/doc/compute_voronoi_atom.html index 2b086741ea..97c360769e 100644 --- a/doc/compute_voronoi_atom.html +++ b/doc/compute_voronoi_atom.html @@ -113,16 +113,16 @@ performed for the first invocation of the compute and then stored. For all following invocations of the compute the number of atoms in each Voronoi cell in the stored tessellation is counted. In this mode the compute returns a per-atom array with 2 columns. The first column -is the number of atoms currently in the Voronoi volume defined by this -atom at the time of the first invocation of the compute (note that the +is the number of atoms currently in the Voronoi volume defined by this +atom at the time of the first invocation of the compute (note that the atom may have moved significantly). The second column contains the -total number of atoms sharing the Voronoi cell of the stored -tessellation at the location of the current atom. Numbers in column one -can be any positive integer including zero, while column two values will -always be greater than zero. Column one data can be used to locate -vacancies (the coordinates are given by the atom coordinates at the -time step when the compute was first invoked), while column two data -can be used to identify interstitial atoms. +total number of atoms sharing the Voronoi cell of the stored +tessellation at the location of the current atom. Numbers in column +one can be any positive integer including zero, while column two +values will always be greater than zero. Column one data can be used +to locate vacancies (the coordinates are given by the atom coordinates +at the time step when the compute was first invoked), while column two +data can be used to identify interstitial atoms.

diff --git a/doc/fix_rigid.html b/doc/fix_rigid.html index 1d8549a090..a71ffde390 100644 --- a/doc/fix_rigid.html +++ b/doc/fix_rigid.html @@ -21,13 +21,21 @@

fix rigid/small command

+

fix rigid/nve/small command +

+

fix rigid/nvt/small command +

+

fix rigid/npt/small command +

+

fix rigid/nph/small command +

Syntax: 

fix ID group-ID style bodystyle args keyword values ...
  • ID, group-ID are documented in fix command -
  • style = rigid or rigid/nve or rigid/nvt or rigid/npt or rigid/nph or rigid/small +
  • style = rigid or rigid/nve or rigid/nvt or rigid/npt or rigid/nph or rigid/small or rigid/nve/small or rigid/nvt/small or rigid/npt/small or rigid/nph/small
  • bodystyle = single or molecule or group @@ -39,7 +47,7 @@
  • zero or more keyword/value pairs may be appended -
  • keyword = langevin or temp or iso or aniso or x or y or z or couple or tparam or pchain or dilate or force or torque or infile or mol +
  • keyword = langevin or temp or iso or aniso or x or y or z or couple or tparam or pchain or dilate or force or torque or infile 
      langevin values = Tstart Tstop Tperiod seed
     Tstart,Tstop = desired temperature at start/stop of run (temperature units)
@@ -70,7 +78,7 @@
     M = which rigid body from 1-Nbody (see asterisk form below)
     xflag,yflag,zflag = off/on if component of center-of-mass torque is active
   infile filename
-    filename = file with per-body values of mass, center-of-mass, moments of inertia
+    filename = file with per-body values of mass, center-of-mass, moments of inertia 
   mol value = template-ID
     template-ID = ID of molecule template specified in a separate molecule command
    @@ -87,7 +95,8 @@ fix 1 polychains rigid/small molecule langevin 1.0 1.0 1.0 428984 fix 2 fluid rigid group 3 clump1 clump2 clump3 torque * off off off fix 1 rods rigid/npt molecule temp 300.0 300.0 100.0 iso 0.5 0.5 10.0 fix 1 particles rigid/npt molecule temp 1.0 1.0 5.0 x 0.5 0.5 1.0 z 0.5 0.5 1.0 couple xz -fix 1 water rigid/nph molecule iso 0.5 0.5 1.0 +fix 1 water rigid/nph molecule iso 0.5 0.5 1.0 +fix 1 particles rigid/npt/small molecule temp 1.0 1.0 1.0 iso 0.5 0.5 1.0

    Description:

    @@ -292,17 +301,17 @@ iterations. The rigid/nve style uses the methods described in the paper by Miller, which are thought to provide better energy conservation than an iterative approach.

    -

    The rigid/nvt style performs constant NVT integration using a -Nose/Hoover thermostat with chains as described originally in -(Hoover) and (Martyna), which thermostats both -the translational and rotational degrees of freedom of the rigid -bodies. The rigid-body algorithm used by rigid/nvt is described in -the paper by Kamberaj. +

    The rigid/nvt and rigid/nvt/small styles performs constant NVT +integration using a Nose/Hoover thermostat with chains as described +originally in (Hoover) and (Martyna), which +thermostats both the translational and rotational degrees of freedom +of the rigid bodies. The rigid-body algorithm used by rigid/nvt +is described in the paper by Kamberaj.

    -

    The rigid/npt and rigid/nph styles perform constant NPT or NPH -integration using a Nose/Hoover barostat with chains. For the NPT -case, the same Nose/Hoover thermostat is also used as with -rigid/nvt. +

    The rigid/npt and rigid/nph (and their /small counterparts) styles +perform constant NPT or NPH integration using a Nose/Hoover barostat +with chains. For the NPT case, the same Nose/Hoover thermostat is also +used as with rigid/nvt.

    The barostat parameters are specified using one or more of the iso, aniso, x, y, z and couple keywords. These keywords give you @@ -312,8 +321,8 @@ they represent are varied together during a constant-pressure simulation. The effects of these keywords are similar to those defined in fix npt/nph

    -

    NOTE: Currently the rigid/npt and rigid/nph styles do not support -triclinic (non-orthongonal) boxes. +

    NOTE: Currently the rigid/npt and rigid/nph (and their /small +counterparts) styles do not support triclinic (non-orthongonal) boxes.

    The target pressures for each of the 6 components of the stress tensor can be specified independently via the x, y, z keywords, which @@ -483,12 +492,10 @@ have to list attributes for every rigid body integrated by fix rigid. Only bodies which the file specifies will have their computed attributes overridden. The file can contain initial blank lines or comment lines starting with "#" which are ignored. The first -non-blank, non-comment line should list N, which is the number of -lines to follow. The N successive lines contain the following -information: +non-blank, non-comment line should list N = the number of lines to +follow. The N successive lines contain the following information:

    -
    N
-ID1 masstotal xcm ycm zcm ixx iyy izz ixy ixz iyz
+
    ID1 masstotal xcm ycm zcm ixx iyy izz ixy ixz iyz
 ID2 masstotal xcm ycm zcm ixx iyy izz ixy ixz iyz
 ...
 IDN masstotal xcm ycm zcm ixx iyy izz ixy ixz iyz 
diff --git a/doc/fix_rigid.txt b/doc/fix_rigid.txt
index ca23c5a74a..62bd5e2ee3 100644
--- a/doc/fix_rigid.txt
+++ b/doc/fix_rigid.txt
@@ -12,13 +12,17 @@ fix rigid/nvt command :h3
 fix rigid/npt command :h3
 fix rigid/nph command :h3
 fix rigid/small command :h3
+fix rigid/nve/small command :h3
+fix rigid/nvt/small command :h3
+fix rigid/npt/small command :h3
+fix rigid/nph/small command :h3
 
 [Syntax:]
 
 fix ID group-ID style bodystyle args keyword values ... :pre
 
 ID, group-ID are documented in "fix"_fix.html command :ulb,l
-style = {rigid} or {rigid/nve} or {rigid/nvt} or {rigid/npt} or {rigid/nph} or {rigid/small} :l
+style = {rigid} or {rigid/nve} or {rigid/nvt} or {rigid/npt} or {rigid/nph} or {rigid/small} or {rigid/nve/small} or {rigid/nvt/small} or {rigid/npt/small} or {rigid/nph/small} :l
 bodystyle = {single} or {molecule} or {group} :l
   {single} args = none
   {molecule} args = none
@@ -27,7 +31,7 @@ bodystyle = {single} or {molecule} or {group} :l
     groupID1, groupID2, ... = list of N group IDs :pre
 
 zero or more keyword/value pairs may be appended :l
-keyword = {langevin} or {temp} or {iso} or {aniso} or {x} or {y} or {z} or {couple} or {tparam} or {pchain} or {dilate} or {force} or {torque} or {infile} or {mol} :l
+keyword = {langevin} or {temp} or {iso} or {aniso} or {x} or {y} or {z} or {couple} or {tparam} or {pchain} or {dilate} or {force} or {torque} or {infile} :l
   {langevin} values = Tstart Tstop Tperiod seed
     Tstart,Tstop = desired temperature at start/stop of run (temperature units)
     Tdamp = temperature damping parameter (time units)
@@ -57,7 +61,7 @@ keyword = {langevin} or {temp} or {iso} or {aniso} or {x} or {y} or {z} or {coup
     M = which rigid body from 1-Nbody (see asterisk form below)
     xflag,yflag,zflag = off/on if component of center-of-mass torque is active
   {infile} filename
-    filename = file with per-body values of mass, center-of-mass, moments of inertia
+    filename = file with per-body values of mass, center-of-mass, moments of inertia 
   {mol} value = template-ID
     template-ID = ID of molecule template specified in a separate "molecule"_molecule.html command :pre
 :ule
@@ -73,8 +77,9 @@ fix 1 polychains rigid/small molecule langevin 1.0 1.0 1.0 428984
 fix 2 fluid rigid group 3 clump1 clump2 clump3 torque * off off off 
 fix 1 rods rigid/npt molecule temp 300.0 300.0 100.0 iso 0.5 0.5 10.0
 fix 1 particles rigid/npt molecule temp 1.0 1.0 5.0 x 0.5 0.5 1.0 z 0.5 0.5 1.0 couple xz
-fix 1 water rigid/nph molecule iso 0.5 0.5 1.0 :pre
-
+fix 1 water rigid/nph molecule iso 0.5 0.5 1.0
+fix 1 particles rigid/npt/small molecule temp 1.0 1.0 1.0 iso 0.5 0.5 1.0 :pre
+	
 [Description:]
 
 Treat one or more sets of atoms as independent rigid bodies.  This
@@ -278,17 +283,17 @@ iterations.  The {rigid/nve} style uses the methods described in the
 paper by "Miller"_#Miller, which are thought to provide better energy
 conservation than an iterative approach.
 
-The {rigid/nvt} style performs constant NVT integration using a
-Nose/Hoover thermostat with chains as described originally in
-"(Hoover)"_#Hoover and "(Martyna)"_#Martyna, which thermostats both
-the translational and rotational degrees of freedom of the rigid
-bodies.  The rigid-body algorithm used by {rigid/nvt} is described in
-the paper by "Kamberaj"_#Kamberaj.
+The {rigid/nvt} and {rigid/nvt/small} styles performs constant NVT 
+integration using a Nose/Hoover thermostat with chains as described 
+originally in "(Hoover)"_#Hoover and "(Martyna)"_#Martyna, which 
+thermostats both the translational and rotational degrees of freedom 
+of the rigid bodies.  The rigid-body algorithm used by {rigid/nvt} 
+is described in the paper by "Kamberaj"_#Kamberaj.
 
-The {rigid/npt} and {rigid/nph} styles perform constant NPT or NPH
-integration using a Nose/Hoover barostat with chains.  For the NPT
-case, the same Nose/Hoover thermostat is also used as with
-{rigid/nvt}.
+The {rigid/npt} and {rigid/nph} (and their /small counterparts) styles 
+perform constant NPT or NPH integration using a Nose/Hoover barostat 
+with chains.  For the NPT case, the same Nose/Hoover thermostat is also 
+used as with {rigid/nvt}.
 
 The barostat parameters are specified using one or more of the {iso},
 {aniso}, {x}, {y}, {z} and {couple} keywords.  These keywords give you
@@ -298,8 +303,8 @@ they represent are varied together during a constant-pressure
 simulation.  The effects of these keywords are similar to those
 defined in "fix npt/nph"_fix_nh.html
 
-NOTE: Currently the {rigid/npt} and {rigid/nph} styles do not support
-triclinic (non-orthongonal) boxes.
+NOTE: Currently the {rigid/npt} and {rigid/nph} (and their /small 
+counterparts) styles do not support triclinic (non-orthongonal) boxes.
 
 The target pressures for each of the 6 components of the stress tensor
 can be specified independently via the {x}, {y}, {z} keywords, which
@@ -469,11 +474,9 @@ have to list attributes for every rigid body integrated by fix rigid.
 Only bodies which the file specifies will have their computed
 attributes overridden.  The file can contain initial blank lines or
 comment lines starting with "#" which are ignored.  The first
-non-blank, non-comment line should list N, which is the number of
-lines to follow.  The N successive lines contain the following
-information:
+non-blank, non-comment line should list N = the number of lines to
+follow.  The N successive lines contain the following information:
 
-N
 ID1 masstotal xcm ycm zcm ixx iyy izz ixy ixz iyz
 ID2 masstotal xcm ycm zcm ixx iyy izz ixy ixz iyz
 ...

From ce8a12bbdc0c487be1ba7ff994a617715bf34ba7 Mon Sep 17 00:00:00 2001
From: sjplimp 
Date: Tue, 29 Apr 2014 22:26:56 +0000
Subject: [PATCH 05/10] git-svn-id:
 svn://svn.icms.temple.edu/lammps-ro/trunk@11840
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 doc/fix_vector.html | 2 +-
 doc/fix_vector.txt  | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/doc/fix_vector.html b/doc/fix_vector.html
index 9acb53eab8..af0bd770cc 100644
--- a/doc/fix_vector.html
+++ b/doc/fix_vector.html
@@ -53,7 +53,7 @@ For example the velocity auto-correlation function (VACF) can be
 time-integrated, to yield a diffusion coefficient, as follows:
 
    
 
    compute         2 all vacf
-fix             5 all vector 1 c_24
+fix             5 all vector 1 c_2[4]
 variable        diff equal dt*trap(f_5)
 thermo_style    custom step v_diff 
 
    
diff --git a/doc/fix_vector.txt b/doc/fix_vector.txt
index 5dd73c4fba..f26d0b4359 100644
--- a/doc/fix_vector.txt
+++ b/doc/fix_vector.txt
@@ -44,7 +44,7 @@ For example the velocity auto-correlation function (VACF) can be
 time-integrated, to yield a diffusion coefficient, as follows:
 
 compute         2 all vacf
-fix             5 all vector 1 c_2[4]
+fix             5 all vector 1 c_2\[4\]
 variable        diff equal dt*trap(f_5)
 thermo_style    custom step v_diff :pre
 

From f3e58c8d0aa7f5b0c952776dcc70081121285341 Mon Sep 17 00:00:00 2001
From: sjplimp 
Date: Tue, 29 Apr 2014 22:43:28 +0000
Subject: [PATCH 06/10] git-svn-id:
 svn://svn.icms.temple.edu/lammps-ro/trunk@11841
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 src/force.cpp | 48 +++++++++++++++++++++++++++++++++++++++++++++++-
 src/force.h   |  1 +
 src/group.cpp | 20 ++++++++------------
 3 files changed, 56 insertions(+), 13 deletions(-)

diff --git a/src/force.cpp b/src/force.cpp
index cef61fc557..44c369acca 100644
--- a/src/force.cpp
+++ b/src/force.cpp
@@ -747,7 +747,14 @@ void Force::boundsbig(char *str, bigint nmax, bigint &nlo, bigint &nhi,
 
 double Force::numeric(const char *file, int line, char *str)
 {
+  if (!str)
+    error->all(file,line,"Expected floating point parameter "
+               "in input script or data file");
   int n = strlen(str);
+  if (n == 0)
+    error->all(file,line,"Expected floating point parameter "
+               "in input script or data file");
+
   for (int i = 0; i < n; i++) {
     if (isdigit(str[i])) continue;
     if (str[i] == '-' || str[i] == '+' || str[i] == '.') continue;
@@ -767,7 +774,14 @@ double Force::numeric(const char *file, int line, char *str)
 
 int Force::inumeric(const char *file, int line, char *str)
 {
+  if (!str) 
+    error->all(file,line,
+               "Expected integer parameter in input script or data file");
   int n = strlen(str);
+  if (n == 0) 
+    error->all(file,line,
+               "Expected integer parameter in input script or data file");
+
   for (int i = 0; i < n; i++) {
     if (isdigit(str[i]) || str[i] == '-' || str[i] == '+') continue;
     error->all(file,line,
@@ -785,14 +799,46 @@ int Force::inumeric(const char *file, int line, char *str)
 
 bigint Force::bnumeric(const char *file, int line, char *str)
 {
+  if (!str) 
+    error->all(file,line,
+               "Expected integer parameter in input script or data file");
   int n = strlen(str);
+  if (n == 0) 
+    error->all(file,line,
+               "Expected integer parameter in input script or data file");
+
   for (int i = 0; i < n; i++) {
     if (isdigit(str[i]) || str[i] == '-' || str[i] == '+') continue;
     error->all(file,line,
                "Expected integer parameter in input script or data file");
   }
 
-  return ATOLL(str);
+  return ATOBIGINT(str);
+}
+
+/* ----------------------------------------------------------------------
+   read a tag integer value from a string
+   generate an error if not a legitimate integer value
+   called by various commands to check validity of their arguments
+------------------------------------------------------------------------- */
+
+tagint Force::tnumeric(const char *file, int line, char *str)
+{
+  if (!str) 
+    error->all(file,line,
+               "Expected integer parameter in input script or data file");
+  int n = strlen(str);
+  if (n == 0) 
+    error->all(file,line,
+               "Expected integer parameter in input script or data file");
+
+  for (int i = 0; i < n; i++) {
+    if (isdigit(str[i]) || str[i] == '-' || str[i] == '+') continue;
+    error->all(file,line,
+               "Expected integer parameter in input script or data file");
+  }
+
+  return ATOTAGINT(str);
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/force.h b/src/force.h
index a9efb08fbf..bf364f253b 100644
--- a/src/force.h
+++ b/src/force.h
@@ -105,6 +105,7 @@ class Force : protected Pointers {
   double numeric(const char *, int, char *);
   int inumeric(const char *, int, char *);
   bigint bnumeric(const char *, int, char *);
+  tagint tnumeric(const char *, int, char *);
 
   FILE *open_potential(const char *);
   const char *potname(const char *);
diff --git a/src/group.cpp b/src/group.cpp
index 4144108f24..3b494a0b42 100644
--- a/src/group.cpp
+++ b/src/group.cpp
@@ -198,18 +198,12 @@ void Group::assign(int narg, char **arg)
       else error->all(FLERR,"Illegal group command");
       
       tagint bound1,bound2;
-      if (sizeof(tagint) == sizeof(smallint))
-        bound1 = force->inumeric(FLERR,arg[3]);
-      else
-        bound1 = force->bnumeric(FLERR,arg[3]);
+      bound1 = force->tnumeric(FLERR,arg[3]);
       bound2 = -1;
 
       if (condition == BETWEEN) {
         if (narg != 5) error->all(FLERR,"Illegal group command");
-        if (sizeof(tagint) == sizeof(smallint))
-          bound2 = force->inumeric(FLERR,arg[4]);
-        else
-          bound2 = force->bnumeric(FLERR,arg[4]);
+        bound2 = force->tnumeric(FLERR,arg[4]);
       } else if (narg != 4) error->all(FLERR,"Illegal group command");
 
       int *attribute = NULL;
@@ -284,13 +278,15 @@ void Group::assign(int narg, char **arg)
 
       for (int iarg = 2; iarg < narg; iarg++) {
         if (strchr(arg[iarg],':')) {
-          start = ATOTAGINT(strtok(arg[iarg],":")); 
-          stop = ATOTAGINT(strtok(NULL,":"));
+          ptr = strtok(arg[iarg],":"); 
+          start = force->tnumeric(FLERR,ptr); 
+          ptr = strtok(NULL,":"); 
+          stop = force->tnumeric(FLERR,ptr); 
           ptr = strtok(NULL,":");
-          if (ptr) delta = ATOTAGINT(ptr);
+          if (ptr) delta = force->tnumeric(FLERR,ptr);
           else delta = 1;
         } else {
-          start = stop = ATOTAGINT(arg[iarg]);
+          start = stop = force->tnumeric(FLERR,arg[iarg]);
           delta = 1;
         }
 

From 99f44f64cc10fc394b041623b3daab2c761ef062 Mon Sep 17 00:00:00 2001
From: sjplimp 
Date: Tue, 29 Apr 2014 23:10:42 +0000
Subject: [PATCH 07/10] git-svn-id:
 svn://svn.icms.temple.edu/lammps-ro/trunk@11842
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 lib/gpu/lal_sw.cpp |  8 +++--
 lib/gpu/lal_sw.cu  | 79 +++++++++++++++++++++++-----------------------
 2 files changed, 46 insertions(+), 41 deletions(-)

diff --git a/lib/gpu/lal_sw.cpp b/lib/gpu/lal_sw.cpp
index 27974874b0..f14b0a3438 100644
--- a/lib/gpu/lal_sw.cpp
+++ b/lib/gpu/lal_sw.cpp
@@ -112,8 +112,12 @@ int SWT::init(const int ntypes, const int nlocal, const int nall, const int max_
   sw3.alloc(nparams,*(this->ucl_device),UCL_READ_ONLY);
   
   for (int i=0; i(cut[i]);
-    dview[i].y=static_cast(cutsq[i]);
+    double sw_cut = cut[i];
+    double sw_cutsq = cutsq[i];
+    if (sw_cutsq>=sw_cut*sw_cut)
+      sw_cutsq=sw_cut*sw_cut-1e-4;
+    dview[i].x=static_cast(sw_cut);
+    dview[i].y=static_cast(sw_cutsq);
     dview[i].z=static_cast(costheta[i]);
     dview[i].w=(numtyp)0;
   }
diff --git a/lib/gpu/lal_sw.cu b/lib/gpu/lal_sw.cu
index 71433cd5c2..c99aa3f675 100644
--- a/lib/gpu/lal_sw.cu
+++ b/lib/gpu/lal_sw.cu
@@ -195,17 +195,17 @@ __kernel void k_sw(const __global numtyp4 *restrict x_,
         numtyp sw_cut=sw3_ijparam.x;
         numtyp sw_cutsq=sw3_ijparam.y;
         numtyp pre_sw_c1=sw_biga*sw_epsilon*sw_powerp*sw_bigb*
-            pow(sw_sigma,sw_powerp);
+            ucl_powr(sw_sigma,sw_powerp);
         numtyp pre_sw_c2=sw_biga*sw_epsilon*sw_powerq*
-            pow(sw_sigma,sw_powerq);
+            ucl_powr(sw_sigma,sw_powerq);
         numtyp pre_sw_c3=sw_biga*sw_epsilon*sw_bigb*
-            pow(sw_sigma,sw_powerp+(numtyp)1.0);
+            ucl_powr(sw_sigma,sw_powerp+(numtyp)1.0);
         numtyp pre_sw_c4=sw_biga*sw_epsilon*
-            pow(sw_sigma,sw_powerq+(numtyp)1.0);
+            ucl_powr(sw_sigma,sw_powerq+(numtyp)1.0);
         numtyp pre_sw_c5=sw_biga*sw_epsilon*sw_bigb*
-            pow(sw_sigma,sw_powerp);
+            ucl_powr(sw_sigma,sw_powerp);
         numtyp pre_sw_c6=sw_biga*sw_epsilon*
-            pow(sw_sigma,sw_powerq);
+            ucl_powr(sw_sigma,sw_powerq);
 
         numtyp r=ucl_sqrt(rsq);
         numtyp rp=ucl_powr(r,-sw_powerp);
@@ -343,7 +343,6 @@ __kernel void k_sw_three_center(const __global numtyp4 *restrict x_,
                                 const int t_per_atom, const int evatom) {
   __local int tpa_sq, n_stride;
   tpa_sq=fast_mul(t_per_atom,t_per_atom);
-  numtyp sw_epsilon, sw_sigma, sw_lambda, sw_gamma;
   numtyp sw_sigma_gamma_ij, sw_cut_ij, sw_sigma_gamma_ik, sw_cut_ik;
   numtyp sw_costheta_ijk, sw_lambda_epsilon_ijk, sw_lambda_epsilon2_ijk;
 
@@ -467,7 +466,6 @@ __kernel void k_sw_three_end(const __global numtyp4 *restrict x_,
                              const int t_per_atom) {
   __local int tpa_sq, n_stride;
   tpa_sq=fast_mul(t_per_atom,t_per_atom);
-  numtyp sw_epsilon, sw_sigma, sw_lambda, sw_gamma;
   numtyp sw_sigma_gamma_ij, sw_cut_ij, sw_sigma_gamma_ik, sw_cut_ik;
   numtyp sw_costheta_ijk, sw_lambda_epsilon_ijk, sw_lambda_epsilon2_ijk;
 
@@ -515,12 +513,14 @@ __kernel void k_sw_three_end(const __global numtyp4 *restrict x_,
 
       if (rsq1 > sw3_ijparam.y) continue;
 
-      numtyp4 sw1_ijparam; fetch4(sw1_ijparam,ijparam,sw1_tex);
-      sw_sigma_gamma_ij=sw1_ijparam.y*sw1_ijparam.w; //sw_sigma*sw_gamma;
-      sw_cut_ij=sw3_ijparam.x;
+      int jiparam=elem2param[jtype*nelements*nelements+itype*nelements+itype]; 
+      numtyp4 sw1_jiparam; fetch4(sw1_jiparam,jiparam,sw1_tex);
+      numtyp4 sw3_jiparam; fetch4(sw3_jiparam,jiparam,sw3_tex);
+      sw_sigma_gamma_ij=sw1_jiparam.y*sw1_jiparam.w; //sw_sigma*sw_gamma;
+      sw_cut_ij=sw3_jiparam.x;
 
       const __global int *nbor_k=dev_nbor+j+nbor_pitch;
-      int numk=*nbor_k; 
+      int numk=*nbor_k;
       if (dev_nbor==dev_packed) {
         nbor_k+=nbor_pitch+fast_mul(j,t_per_atom-1);
         k_end=nbor_k+fast_mul(numk/t_per_atom,n_stride)+(numk & (t_per_atom-1));
@@ -541,25 +541,25 @@ __kernel void k_sw_three_end(const __global numtyp4 *restrict x_,
         numtyp4 kx; fetch4(kx,k,pos_tex);
         int ktype=kx.w;
         ktype=map[ktype];
-        int ikparam=elem2param[itype*nelements*nelements+ktype*nelements+ktype]; 
+        int jkparam=elem2param[jtype*nelements*nelements+ktype*nelements+ktype]; 
 
         numtyp delr2x = kx.x - jx.x;
         numtyp delr2y = kx.y - jx.y;
         numtyp delr2z = kx.z - jx.z;
         numtyp rsq2 = delr2x*delr2x + delr2y*delr2y + delr2z*delr2z;
-        numtyp4 sw3_ikparam; fetch4(sw3_ikparam,ikparam,sw3_tex);
+        numtyp4 sw3_jkparam; fetch4(sw3_jkparam,jkparam,sw3_tex);
 
-        if (rsq2 < sw3_ikparam.y) {
-          numtyp4 sw1_ikparam; fetch4(sw1_ikparam,ikparam,sw1_tex);
-          sw_sigma_gamma_ik=sw1_ikparam.y*sw1_ikparam.w; //sw_sigma*sw_gamma;
-          sw_cut_ik=sw3_ikparam.x;
+        if (rsq2 < sw3_jkparam.y) {
+          numtyp4 sw1_jkparam; fetch4(sw1_jkparam,jkparam,sw1_tex);
+          sw_sigma_gamma_ik=sw1_jkparam.y*sw1_jkparam.w; //sw_sigma*sw_gamma;
+          sw_cut_ik=sw3_jkparam.x;
 
-          int ijkparam=elem2param[itype*nelements*nelements+jtype*nelements+ktype];
-          numtyp4 sw1_ijkparam; fetch4(sw1_ijkparam,ijkparam,sw1_tex);
-          sw_lambda_epsilon_ijk=sw1_ijkparam.x*sw1_ijkparam.z; //sw_lambda*sw_epsilon;
+          int jikparam=elem2param[jtype*nelements*nelements+itype*nelements+ktype];
+          numtyp4 sw1_jikparam; fetch4(sw1_jikparam,jikparam,sw1_tex);
+          sw_lambda_epsilon_ijk=sw1_jikparam.x*sw1_jikparam.z; //sw_lambda*sw_epsilon;
           sw_lambda_epsilon2_ijk=(numtyp)2.0*sw_lambda_epsilon_ijk;
-          numtyp4 sw3_ijkparam; fetch4(sw3_ijkparam,ijkparam,sw3_tex);
-          sw_costheta_ijk=sw3_ijkparam.z;
+          numtyp4 sw3_jikparam; fetch4(sw3_jikparam,jikparam,sw3_tex);
+          sw_costheta_ijk=sw3_jikparam.z;
 
           numtyp fjx, fjy, fjz;
           //if (evatom==0) {
@@ -602,7 +602,6 @@ __kernel void k_sw_three_end_vatom(const __global numtyp4 *restrict x_,
                              const int t_per_atom) {
   __local int tpa_sq, n_stride;
   tpa_sq=fast_mul(t_per_atom,t_per_atom);
-  numtyp sw_epsilon, sw_sigma, sw_lambda, sw_gamma;
   numtyp sw_sigma_gamma_ij, sw_cut_ij, sw_sigma_gamma_ik, sw_cut_ik;
   numtyp sw_costheta_ijk, sw_lambda_epsilon_ijk, sw_lambda_epsilon2_ijk;
 
@@ -650,10 +649,12 @@ __kernel void k_sw_three_end_vatom(const __global numtyp4 *restrict x_,
 
       if (rsq1 > sw3_ijparam.y) continue;
 
-      numtyp4 sw1_ijparam; fetch4(sw1_ijparam,ijparam,sw1_tex);
-      sw_sigma_gamma_ij=sw1_ijparam.y*sw1_ijparam.w; //sw_sigma*sw_gamma;
-      sw_cut_ij=sw3_ijparam.x;
-        
+      int jiparam=elem2param[jtype*nelements*nelements+itype*nelements+itype]; 
+      numtyp4 sw1_jiparam; fetch4(sw1_jiparam,jiparam,sw1_tex);
+      numtyp4 sw3_jiparam; fetch4(sw3_jiparam,jiparam,sw3_tex);
+      sw_sigma_gamma_ij=sw1_jiparam.y*sw1_jiparam.w; //sw_sigma*sw_gamma;
+      sw_cut_ij=sw3_jiparam.x;
+       
       const __global int *nbor_k=dev_nbor+j+nbor_pitch;
       int numk=*nbor_k;
       if (dev_nbor==dev_packed) {
@@ -676,25 +677,25 @@ __kernel void k_sw_three_end_vatom(const __global numtyp4 *restrict x_,
         numtyp4 kx; fetch4(kx,k,pos_tex);
         int ktype=kx.w;
         ktype=map[ktype];
-        int ikparam=elem2param[itype*nelements*nelements+ktype*nelements+ktype]; 
-        numtyp4 sw3_ikparam; fetch4(sw3_ikparam,ikparam,sw3_tex);
+        int jkparam=elem2param[jtype*nelements*nelements+ktype*nelements+ktype]; 
+        numtyp4 sw3_jkparam; fetch4(sw3_jkparam,jkparam,sw3_tex);
 
         numtyp delr2x = kx.x - jx.x;
         numtyp delr2y = kx.y - jx.y;
         numtyp delr2z = kx.z - jx.z;
         numtyp rsq2 = delr2x*delr2x + delr2y*delr2y + delr2z*delr2z;
 
-        if (rsq2 < sw3_ikparam.y) {
-          numtyp4 sw1_ikparam; fetch4(sw1_ikparam,ikparam,sw1_tex);
-          sw_sigma_gamma_ik=sw1_ikparam.y*sw1_ikparam.w; //sw_sigma*sw_gamma;
-          sw_cut_ik=sw3_ikparam.x;
+        if (rsq2 < sw3_jkparam.y) {
+          numtyp4 sw1_jkparam; fetch4(sw1_jkparam,jkparam,sw1_tex);
+          sw_sigma_gamma_ik=sw1_jkparam.y*sw1_jkparam.w; //sw_sigma*sw_gamma;
+          sw_cut_ik=sw3_jkparam.x;
 
-          int ijkparam=elem2param[itype*nelements*nelements+jtype*nelements+ktype];
-          numtyp4 sw1_ijkparam; fetch4(sw1_ijkparam,ijkparam,sw1_tex);
-          sw_lambda_epsilon_ijk=sw1_ijkparam.x*sw1_ijkparam.z; //sw_lambda*sw_epsilon;
+          int jikparam=elem2param[jtype*nelements*nelements+itype*nelements+ktype];
+          numtyp4 sw1_jikparam; fetch4(sw1_jikparam,jikparam,sw1_tex);
+          sw_lambda_epsilon_ijk=sw1_jikparam.x*sw1_jikparam.z; //sw_lambda*sw_epsilon;
           sw_lambda_epsilon2_ijk=(numtyp)2.0*sw_lambda_epsilon_ijk;
-          numtyp4 sw3_ijkparam; fetch4(sw3_ijkparam,ijkparam,sw3_tex);
-          sw_costheta_ijk=sw3_ijkparam.z;
+          numtyp4 sw3_jikparam; fetch4(sw3_jikparam,jikparam,sw3_tex);
+          sw_costheta_ijk=sw3_jikparam.z;
 
           numtyp fjx, fjy, fjz, fkx, fky, fkz;
           threebody(delr1x,delr1y,delr1z,eflag,energy);

From a9ccfa94bf20fbf70b27f5861e739284f84c7202 Mon Sep 17 00:00:00 2001
From: sjplimp 
Date: Tue, 29 Apr 2014 23:25:54 +0000
Subject: [PATCH 08/10] git-svn-id:
 svn://svn.icms.temple.edu/lammps-ro/trunk@11843
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 src/RIGID/fix_rigid_small.cpp | 33 +++++++++++++++++++++++++++++++--
 1 file changed, 31 insertions(+), 2 deletions(-)

diff --git a/src/RIGID/fix_rigid_small.cpp b/src/RIGID/fix_rigid_small.cpp
index 285cd6e710..b2974de386 100644
--- a/src/RIGID/fix_rigid_small.cpp
+++ b/src/RIGID/fix_rigid_small.cpp
@@ -313,6 +313,15 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) :
     onemol->compute_inertia();
   }
 
+  // set pstat_flag
+
+  pstat_flag = 0;
+  for (int i = 0; i < 3; i++) 
+    if (p_flag[i]) pstat_flag = 1;
+
+  if (pcouple == XYZ || (domain->dimension == 2 && pcouple == XY)) pstyle = ISO;
+  else pstyle = ANISO;
+
   // create rigid bodies based on molecule ID
   // sets bodytag for owned atoms
   // body attributes are computed later by setup_bodies()
@@ -2742,7 +2751,7 @@ int FixRigidSmall::pack_comm(int n, int *list, double *buf,
                              int pbc_flag, int *pbc)
 {
   int i,j;
-  double *xcm,*vcm,*quat,*omega,*ex_space,*ey_space,*ez_space;
+  double *xcm,*vcm,*quat,*omega,*ex_space,*ey_space,*ez_space,*conjqm;
 
   int m = 0;
 
@@ -2779,6 +2788,11 @@ int FixRigidSmall::pack_comm(int n, int *list, double *buf,
       buf[m++] = ez_space[0];
       buf[m++] = ez_space[1];
       buf[m++] = ez_space[2];
+      conjqm = body[bodyown[j]].conjqm;
+      buf[m++] = conjqm[0];
+      buf[m++] = conjqm[1];
+      buf[m++] = conjqm[2];
+      buf[m++] = conjqm[3];
     }
 
   } else if (commflag == FINAL) {
@@ -2793,6 +2807,11 @@ int FixRigidSmall::pack_comm(int n, int *list, double *buf,
       buf[m++] = omega[0];
       buf[m++] = omega[1];
       buf[m++] = omega[2];
+      conjqm = body[bodyown[j]].conjqm;
+      buf[m++] = conjqm[0];
+      buf[m++] = conjqm[1];
+      buf[m++] = conjqm[2];
+      buf[m++] = conjqm[3];
     }
 
   } else if (commflag == FULL_BODY) {
@@ -2819,7 +2838,7 @@ int FixRigidSmall::pack_comm(int n, int *list, double *buf,
 void FixRigidSmall::unpack_comm(int n, int first, double *buf)
 {
   int i,j,last;
-  double *xcm,*vcm,*quat,*omega,*ex_space,*ey_space,*ez_space;
+  double *xcm,*vcm,*quat,*omega,*ex_space,*ey_space,*ez_space,*conjqm;
 
   int m = 0;
   last = first + n;
@@ -2856,6 +2875,11 @@ void FixRigidSmall::unpack_comm(int n, int first, double *buf)
       ez_space[0] = buf[m++];
       ez_space[1] = buf[m++];
       ez_space[2] = buf[m++];
+      conjqm = body[bodyown[i]].conjqm;
+      conjqm[0] = buf[m++];
+      conjqm[1] = buf[m++];
+      conjqm[2] = buf[m++];
+      conjqm[3] = buf[m++];
     }
 
   } else if (commflag == FINAL) {
@@ -2869,6 +2893,11 @@ void FixRigidSmall::unpack_comm(int n, int first, double *buf)
       omega[0] = buf[m++];
       omega[1] = buf[m++];
       omega[2] = buf[m++];
+      conjqm = body[bodyown[i]].conjqm;
+      conjqm[0] = buf[m++];
+      conjqm[1] = buf[m++];
+      conjqm[2] = buf[m++];
+      conjqm[3] = buf[m++];
     }
 
   } else if (commflag == FULL_BODY) {

From 4f6dcceebfdb8153af4c6706f321c23e308f987c Mon Sep 17 00:00:00 2001
From: sjplimp 
Date: Tue, 29 Apr 2014 23:26:30 +0000
Subject: [PATCH 09/10] git-svn-id:
 svn://svn.icms.temple.edu/lammps-ro/trunk@11844
 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 74bdfbf251..349cbaee43 100644
--- a/src/version.h
+++ b/src/version.h
@@ -1 +1 @@
-#define LAMMPS_VERSION "18 Apr 2014"
+#define LAMMPS_VERSION "29 Apr 2014"

From 2cf0a3bafa5a5ca1a7e685ec56afca6fab48d1ef Mon Sep 17 00:00:00 2001
From: sjplimp 
Date: Tue, 29 Apr 2014 23:26:34 +0000
Subject: [PATCH 10/10] git-svn-id:
 svn://svn.icms.temple.edu/lammps-ro/trunk@11845
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 doc/Manual.html | 4 ++--
 doc/Manual.txt  | 4 ++--
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/doc/Manual.html b/doc/Manual.html
index 34f552dbb6..9ede33e273 100644
--- a/doc/Manual.html
+++ b/doc/Manual.html
@@ -1,7 +1,7 @@
 
 
 LAMMPS Users Manual
-
+
 
 
 
@@ -22,7 +22,7 @@
 
 
    LAMMPS Documentation 
 
    
-
    18 Apr 2014 version 
+
    29 Apr 2014 version 
 
    
 
    Version info: 
 
    
diff --git a/doc/Manual.txt b/doc/Manual.txt
index f82d267a0b..006864a8d2 100644
--- a/doc/Manual.txt
+++ b/doc/Manual.txt
@@ -1,6 +1,6 @@
 
 LAMMPS Users Manual
-
+
 
 
 
@@ -18,7 +18,7 @@
 
    
 
 LAMMPS Documentation :c,h3
-18 Apr 2014 version :c,h4
+29 Apr 2014 version :c,h4
 
 Version info: :h4