diff --git a/src/fix_rigid.cpp b/src/fix_rigid.cpp index 5bb655e834..5514a48482 100644 --- a/src/fix_rigid.cpp +++ b/src/fix_rigid.cpp @@ -39,6 +39,8 @@ using namespace FixConst; using namespace MathConst; enum{SINGLE,MOLECULE,GROUP}; +enum{NONE,XYZ,XY,YZ,XZ}; +enum{ISO,ANISO,TRICLINIC}; #define MAXLINE 256 #define CHUNK 1024 @@ -241,14 +243,25 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) : int seed; langflag = 0; - tempflag = 0; - pressflag = 0; + tstat_flag = 0; + pstat_flag = 0; + allremap = 1; + id_dilate = NULL; t_chain = 10; t_iter = 1; t_order = 3; p_chain = 10; infile = NULL; + pcouple = NONE; + pstyle = ANISO; + dimension = domain->dimension; + + for (int i = 0; i < 3; i++) { + p_start[i] = p_stop[i] = p_period[i] = 0.0; + p_flag[i] = 0; + } + while (iarg < narg) { if (strcmp(arg[iarg],"force") == 0) { if (iarg+5 > narg) error->all(FLERR,"Illegal fix rigid command"); @@ -330,34 +343,104 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) : if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid command"); if (strcmp(style,"rigid/nvt") != 0 && strcmp(style,"rigid/npt") != 0) error->all(FLERR,"Illegal fix rigid command"); - tempflag = 1; + tstat_flag = 1; t_start = atof(arg[iarg+1]); t_stop = atof(arg[iarg+2]); t_period = atof(arg[iarg+3]); iarg += 4; - } else if (strcmp(arg[iarg],"press") == 0) { + } else if (strcmp(arg[iarg],"iso") == 0) { if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid command"); - if (strcmp(style,"rigid/npt") != 0) - error->all(FLERR,"Illegal fix rigid command"); - pressflag = 1; - p_start = atof(arg[iarg+1]); - p_stop = atof(arg[iarg+2]); - p_period = atof(arg[iarg+3]); + if (strcmp(style,"rigid/npt") != 0 && strcmp(style,"rigid/nph") != 0) + error->all(FLERR,"Illegal fix rigid command"); + pcouple = XYZ; + p_start[0] = p_start[1] = p_start[2] = atof(arg[iarg+1]); + p_stop[0] = p_stop[1] = p_stop[2] = atof(arg[iarg+2]); + p_period[0] = p_period[1] = p_period[2] = atof(arg[iarg+3]); + p_flag[0] = p_flag[1] = p_flag[2] = 1; + if (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 command"); + if (strcmp(style,"rigid/npt") != 0 && strcmp(style,"rigid/nph") != 0) + error->all(FLERR,"Illegal fix rigid command"); + p_start[0] = p_start[1] = p_start[2] = atof(arg[iarg+1]); + p_stop[0] = p_stop[1] = p_stop[2] = atof(arg[iarg+2]); + p_period[0] = p_period[1] = p_period[2] = atof(arg[iarg+3]); + p_flag[0] = p_flag[1] = p_flag[2] = 1; + if (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 command"); + p_start[0] = atof(arg[iarg+1]); + p_stop[0] = atof(arg[iarg+2]); + p_period[0] = atof(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 command"); + p_start[1] = atof(arg[iarg+1]); + p_stop[1] = atof(arg[iarg+2]); + p_period[1] = atof(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 command"); + p_start[2] = atof(arg[iarg+1]); + p_stop[2] = atof(arg[iarg+2]); + p_period[2] = atof(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 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 command"); + iarg += 2; + + } else if (strcmp(arg[iarg],"dilate") == 0) { + if (iarg+2 > narg) + error->all(FLERR,"Illegal fix rigid 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 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 command"); - if (strcmp(style,"rigid/nvt") != 0) + if (strcmp(style,"rigid/nvt") != 0 && strcmp(style,"rigid/npt") != 0) error->all(FLERR,"Illegal fix rigid command"); t_chain = atoi(arg[iarg+1]); t_iter = atoi(arg[iarg+2]); t_order = atoi(arg[iarg+3]); iarg += 4; - } else if (strcmp(arg[iarg],"pparam") == 0) { + } else if (strcmp(arg[iarg],"pchain") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix rigid command"); - if (strcmp(style,"rigid/npt") != 0) + if (strcmp(style,"rigid/npt") != 0 && strcmp(style,"rigid/nph") != 0) error->all(FLERR,"Illegal fix rigid command"); p_chain = atoi(arg[iarg+1]); iarg += 2; @@ -372,8 +455,16 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) : } else error->all(FLERR,"Illegal fix rigid command"); } - t_target = t_start; + + // set pstat_flag + pstat_flag = 0; + for (int i = 0; i < 3; i++) + if (p_flag[i]) pstat_flag = 1; + + if (pcouple == XYZ || (dimension == 2 && pcouple == XY)) pstyle = ISO; + else pstyle = ANISO; + // initialize Marsaglia RNG with processor-unique seed if (langflag) random = new RanMars(lmp,seed + me); @@ -1651,7 +1742,7 @@ void FixRigid::setup_bodies() // then remap the xcm of each body back into simulation box if needed for (ibody = 0; ibody < nbody; ibody++) - imagebody[ibody] = ((tagint) IMGMAX << IMG2BITS) | + imagebody[ibody] = ((tagint) IMGMAX << IMG2BITS) | ((tagint) IMGMAX << IMGBITS) | IMGMAX; pre_neighbor(); @@ -2003,7 +2094,7 @@ void FixRigid::readfile(int which, double *vec, double **array, int *inbody) FILE *fp; char *eof,*start,*next,*buf; char line[MAXLINE]; - + if (me == 0) { fp = fopen(infile,"r"); if (fp == NULL) { @@ -2052,7 +2143,7 @@ void FixRigid::readfile(int which, double *vec, double **array, int *inbody) if (nwords != ATTRIBUTE_PERBODY) error->all(FLERR,"Incorrect rigid body format in fix rigid file"); - + // loop over lines of rigid body attributes // tokenize the line into values // id = rigid body ID @@ -2062,19 +2153,19 @@ void FixRigid::readfile(int which, double *vec, double **array, int *inbody) for (int i = 0; i < nchunk; i++) { next = strchr(buf,'\n'); - + values[0] = strtok(buf," \t\n\r\f"); for (j = 1; j < nwords; j++) values[j] = strtok(NULL," \t\n\r\f"); - + id = atoi(values[0]); if (rstyle == MOLECULE) { - if (id <= 0 || id > maxmol) + if (id <= 0 || id > maxmol) error->all(FLERR,"Invalid rigid body ID in fix rigid file"); id = mol2body[id]; } else id--; - if (id < 0 || id >= nbody) + if (id < 0 || id >= nbody) error->all(FLERR,"Invalid rigid body ID in fix rigid file"); inbody[id] = 1; @@ -2094,7 +2185,7 @@ void FixRigid::readfile(int which, double *vec, double **array, int *inbody) buf = next + 1; } - + nread += nchunk; } diff --git a/src/fix_rigid.h b/src/fix_rigid.h index 3888d88bd1..9c0dcf7d2e 100644 --- a/src/fix_rigid.h +++ b/src/fix_rigid.h @@ -41,7 +41,8 @@ class FixRigid : public Fix { void initial_integrate_respa(int, int, int); void final_integrate_respa(int, int); virtual double compute_scalar(); - + virtual int modify_param(int, char **) {return 0;} + double memory_usage(); void grow_arrays(int); void copy_arrays(int, int); @@ -55,7 +56,7 @@ class FixRigid : public Fix { void reset_dt(); virtual void *extract(const char*,int &); double compute_array(int, int); - + protected: int me,nprocs; double dtv,dtf,dtq; @@ -67,6 +68,7 @@ class FixRigid : public Fix { int firstflag; // 1 for first-time setup of rigid bodies char *infile; // file to read rigid body attributes from + int dimension; // # of dimensions int nbody; // # of rigid bodies int *nrigid; // # of atoms in each rigid body int *mol2body; // convert mol-ID to rigid body index @@ -103,16 +105,22 @@ class FixRigid : public Fix { double tfactor; // scale factor on temperature of rigid bodies int langflag; // 0/1 = no/yes Langevin thermostat - int tempflag; // NVT settings + int tstat_flag; // NVT settings double t_start,t_stop,t_target; double t_period,t_freq; int t_chain,t_iter,t_order; - int pressflag; // NPT settings - double p_start,p_stop; - double p_period,p_freq; + int pstat_flag; // NPT settings + 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 + class RanMars *random; class AtomVecEllipsoid *avec_ellipsoid; class AtomVecLine *avec_line; diff --git a/src/fix_rigid_nh.cpp b/src/fix_rigid_nh.cpp new file mode 100644 index 0000000000..dbedfe1530 --- /dev/null +++ b/src/fix_rigid_nh.cpp @@ -0,0 +1,1436 @@ +/* ---------------------------------------------------------------------- + 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: Tony Sheh (U Michigan), Trung Dac Nguyen (U Michigan) + 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.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 + +/* ---------------------------------------------------------------------- */ + +FixRigidNH::FixRigidNH(LAMMPS *lmp, int narg, char **arg) : + FixRigid(lmp, narg, arg) +{ + // error checks: could be moved up to FixRigid + + 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 npt/nph period must be > 0.0"); + if (domain->nonperiodic) + error->all(FLERR, + "Cannot use fix rigid npt/nph with non-periodic dimension"); + + if (dimension == 2 && p_flag[2]) + error->all(FLERR,"Invalid fix rigid npt/nph command for a 2d simulation"); + if (dimension == 2 && (pcouple == YZ || pcouple == XZ)) + error->all(FLERR,"Invalid fix rigid npt/nph command for a 2d simulation"); + + if (pcouple == XYZ && (p_flag[0] == 0 || p_flag[1] == 0)) + error->all(FLERR,"Invalid fix rigid npt/nph command pressure settings"); + if (pcouple == XYZ && dimension == 3 && p_flag[2] == 0) + error->all(FLERR,"Invalid fix rigid npt/nph command pressure settings"); + if (pcouple == XY && (p_flag[0] == 0 || p_flag[1] == 0)) + error->all(FLERR,"Invalid fix rigid npt/nph command pressure settings"); + if (pcouple == YZ && (p_flag[1] == 0 || p_flag[2] == 0)) + error->all(FLERR,"Invalid fix rigid npt/nph command pressure settings"); + if (pcouple == XZ && (p_flag[0] == 0 || p_flag[2] == 0)) + error->all(FLERR,"Invalid fix rigid npt/nph command pressure settings"); + + // require periodicity in tensile dimension + + if (p_flag[0] && domain->xperiodic == 0) + error->all(FLERR, + "Cannot use fix rigid npt/nph on a non-periodic dimension"); + if (p_flag[1] && domain->yperiodic == 0) + error->all(FLERR, + "Cannot use fix rigid npt/nph on a non-periodic dimension"); + if (p_flag[2] && domain->zperiodic == 0) + error->all(FLERR, + "Cannot use fix rigid npt/nph on a non-periodic dimension"); + + if (pcouple == XYZ && 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 npt/nph pressure settings"); + if (pcouple == XYZ && 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 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 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 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 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 nvt/npt/nph damping parameters must be > 0.0"); + + // memory allocation and initialization + + memory->create(conjqm,nbody,4,"rigid_nh:conjqm"); + 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; +} + +/* ---------------------------------------------------------------------- */ + +FixRigidNH::~FixRigidNH() +{ + memory->destroy(conjqm); + 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 FixRigidNH::setmask() +{ + int mask = 0; + mask = FixRigid::setmask(); + if (tstat_flag || pstat_flag) mask |= THERMO_ENERGY; + + return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixRigidNH::init() +{ + FixRigid::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; + + if (force->kspace) kspace_flag = 1; + else kspace_flag = 0; + + nf_t = nf_r = dimension * nbody; + for (int ibody = 0; ibody < nbody; ibody++) + for (int k = 0; k < domain->dimension; k++) + if (fabs(inertia[ibody][k]) < EPSILON) nf_r--; + + // 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]; + } + } + + g_f = nf_t + nf_r; + onednft = 1.0 + (double)(dimension) / (double)g_f; + onednfr = (double) (dimension) / (double)g_f; + + 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 FixRigidNH::setup(int vflag) +{ + FixRigid::setup(vflag); + + double mbody[3]; + akin_t = akin_r = 0.0; + for (int ibody = 0; ibody < nbody; ibody++) { + MathExtra::transpose_matvec(ex_space[ibody],ey_space[ibody],ez_space[ibody], + angmom[ibody],mbody); + MathExtra::quatvec(quat[ibody],mbody,conjqm[ibody]); + conjqm[ibody][0] *= 2.0; + conjqm[ibody][1] *= 2.0; + conjqm[ibody][2] *= 2.0; + conjqm[ibody][3] *= 2.0; + + if (tstat_flag || pstat_flag) { + akin_t += masstotal[ibody]*(vcm[ibody][0]*vcm[ibody][0] + + vcm[ibody][1]*vcm[ibody][1] + vcm[ibody][2]*vcm[ibody][2]); + akin_r += angmom[ibody][0]*omega[ibody][0] + + angmom[ibody][1]*omega[ibody][1] + angmom[ibody][2]*omega[ibody][2]; + } + } + + // 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 FixRigidNH::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) { + akin_t = akin_r = 0.0; + tmp = exp(-1.0 * dtq * eta_dot_t[0]); + scale_t[0] = scale_t[1] = scale_t[2] = tmp; + tmp = exp(-1.0 * dtq * eta_dot_r[0]); + scale_r = tmp; + } + + if (pstat_flag) { + akin_t = akin_r = 0.0; + 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 < nbody; ibody++) { + + // step 1.1 - update vcm by 1/2 step + + dtfm = dtf / masstotal[ibody]; + vcm[ibody][0] += dtfm * fcm[ibody][0] * fflag[ibody][0]; + vcm[ibody][1] += dtfm * fcm[ibody][1] * fflag[ibody][1]; + vcm[ibody][2] += dtfm * fcm[ibody][2] * fflag[ibody][2]; + + if (tstat_flag || pstat_flag) { + vcm[ibody][0] *= scale_t[0]; + vcm[ibody][1] *= scale_t[1]; + vcm[ibody][2] *= scale_t[2]; + + tmp = vcm[ibody][0]*vcm[ibody][0] + vcm[ibody][1]*vcm[ibody][1] + + vcm[ibody][2]*vcm[ibody][2]; + akin_t += masstotal[ibody]*tmp; + } + + // step 1.2 - update xcm by full step + + if (!pstat_flag) { + xcm[ibody][0] += dtv * vcm[ibody][0]; + xcm[ibody][1] += dtv * vcm[ibody][1]; + xcm[ibody][2] += dtv * vcm[ibody][2]; + } else { + xcm[ibody][0] += scale_v[0] * vcm[ibody][0]; + xcm[ibody][1] += scale_v[1] * vcm[ibody][1]; + xcm[ibody][2] += scale_v[2] * vcm[ibody][2]; + } + + // step 1.3 - apply torque (body coords) to quaternion momentum + + torque[ibody][0] *= tflag[ibody][0]; + torque[ibody][1] *= tflag[ibody][1]; + torque[ibody][2] *= tflag[ibody][2]; + + MathExtra::transpose_matvec(ex_space[ibody],ey_space[ibody],ez_space[ibody], + torque[ibody],tbody); + MathExtra::quatvec(quat[ibody],tbody,fquat); + + conjqm[ibody][0] += dtf2 * fquat[0]; + conjqm[ibody][1] += dtf2 * fquat[1]; + conjqm[ibody][2] += dtf2 * fquat[2]; + conjqm[ibody][3] += dtf2 * fquat[3]; + + if (tstat_flag || pstat_flag) { + conjqm[ibody][0] *= scale_r; + conjqm[ibody][1] *= scale_r; + conjqm[ibody][2] *= scale_r; + conjqm[ibody][3] *= scale_r; + } + + // step 1.4 to 1.13 - use no_squish rotate to update p and q + + no_squish_rotate(3,conjqm[ibody],quat[ibody],inertia[ibody],dtq); + no_squish_rotate(2,conjqm[ibody],quat[ibody],inertia[ibody],dtq); + no_squish_rotate(1,conjqm[ibody],quat[ibody],inertia[ibody],dtv); + no_squish_rotate(2,conjqm[ibody],quat[ibody],inertia[ibody],dtq); + no_squish_rotate(3,conjqm[ibody],quat[ibody],inertia[ibody],dtq); + + // update exyz_space + // transform p back to angmom + // update angular velocity + + MathExtra::q_to_exyz(quat[ibody],ex_space[ibody],ey_space[ibody], + ez_space[ibody]); + MathExtra::invquatvec(quat[ibody],conjqm[ibody],mbody); + MathExtra::matvec(ex_space[ibody],ey_space[ibody],ez_space[ibody], + mbody,angmom[ibody]); + + angmom[ibody][0] *= 0.5; + angmom[ibody][1] *= 0.5; + angmom[ibody][2] *= 0.5; + + MathExtra::angmom_to_omega(angmom[ibody],ex_space[ibody],ey_space[ibody], + ez_space[ibody],inertia[ibody],omega[ibody]); + + if (tstat_flag || pstat_flag) { + akin_r += angmom[ibody][0]*omega[ibody][0] + + angmom[ibody][1]*omega[ibody][1] + angmom[ibody][2]*omega[ibody][2]; + } + } + + // virial setup before call to set_xv + + if (vflag) v_setup(vflag); + else evflag = 0; + + // 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 FixRigidNH::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)); + + akin_t = akin_r = 0.0; + } + + // sum over atoms to get force and torque on rigid body + + int *image = atom->image; + double **x = atom->x; + double **f = atom->f; + int nlocal = atom->nlocal; + + double xprd = domain->xprd; + double yprd = domain->yprd; + double zprd = domain->zprd; + if (triclinic) { + xy = domain->xy; + xz = domain->xz; + yz = domain->yz; + } + + int xbox,ybox,zbox; + double xunwrap,yunwrap,zunwrap,dx,dy,dz; + for (ibody = 0; ibody < nbody; ibody++) + for (i = 0; i < 6; i++) sum[ibody][i] = 0.0; + + for (i = 0; i < nlocal; i++) { + if (body[i] < 0) continue; + ibody = body[i]; + + sum[ibody][0] += f[i][0]; + sum[ibody][1] += f[i][1]; + sum[ibody][2] += f[i][2]; + + xbox = (image[i] & IMGMASK) - IMGMAX; + ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; + zbox = (image[i] >> IMG2BITS) - IMGMAX; + + if (triclinic == 0) { + xunwrap = x[i][0] + xbox*xprd; + yunwrap = x[i][1] + ybox*yprd; + zunwrap = x[i][2] + zbox*zprd; + } else { + xunwrap = x[i][0] + xbox*xprd + ybox*xy + zbox*xz; + yunwrap = x[i][1] + ybox*yprd + zbox*yz; + zunwrap = x[i][2] + zbox*zprd; + } + + dx = xunwrap - xcm[ibody][0]; + dy = yunwrap - xcm[ibody][1]; + dz = zunwrap - xcm[ibody][2]; + + sum[ibody][3] += dy*f[i][2] - dz*f[i][1]; + sum[ibody][4] += dz*f[i][0] - dx*f[i][2]; + sum[ibody][5] += dx*f[i][1] - dy*f[i][0]; + } + + // extended particles add their torque to torque of body + + if (extended) { + double **torque_one = atom->torque; + + for (i = 0; i < nlocal; i++) { + if (body[i] < 0) continue; + ibody = body[i]; + + if (eflags[i] & TORQUE) { + sum[ibody][3] += torque_one[i][0]; + sum[ibody][4] += torque_one[i][1]; + sum[ibody][5] += torque_one[i][2]; + } + } + } + + MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world); + + // update vcm and angmom + // include Langevin thermostat forces + // fflag,tflag = 0 for some dimensions in 2d + + for (ibody = 0; ibody < nbody; ibody++) { + fcm[ibody][0] = all[ibody][0] + langextra[ibody][0]; + fcm[ibody][1] = all[ibody][1] + langextra[ibody][1]; + fcm[ibody][2] = all[ibody][2] + langextra[ibody][2]; + torque[ibody][0] = all[ibody][3] + langextra[ibody][3]; + torque[ibody][1] = all[ibody][4] + langextra[ibody][4]; + torque[ibody][2] = all[ibody][5] + langextra[ibody][5]; + + // update vcm by 1/2 step + + dtfm = dtf / masstotal[ibody]; + if (tstat_flag || pstat_flag) { + vcm[ibody][0] *= scale_t[0]; + vcm[ibody][1] *= scale_t[1]; + vcm[ibody][2] *= scale_t[2]; + } + + vcm[ibody][0] += dtfm * fcm[ibody][0] * fflag[ibody][0]; + vcm[ibody][1] += dtfm * fcm[ibody][1] * fflag[ibody][1]; + vcm[ibody][2] += dtfm * fcm[ibody][2] * fflag[ibody][2]; + + if (pstat_flag) { + tmp = vcm[ibody][0]*vcm[ibody][0] + vcm[ibody][1]*vcm[ibody][1] + + vcm[ibody][2]*vcm[ibody][2]; + akin_t += masstotal[ibody]*tmp; + } + + // update conjqm, then transform to angmom, set velocity again + // virial is already setup from initial_integrate + + torque[ibody][0] *= tflag[ibody][0]; + torque[ibody][1] *= tflag[ibody][1]; + torque[ibody][2] *= tflag[ibody][2]; + + MathExtra::transpose_matvec(ex_space[ibody],ey_space[ibody], + ez_space[ibody],torque[ibody],tbody); + MathExtra::quatvec(quat[ibody],tbody,fquat); + + if (tstat_flag || pstat_flag) { + conjqm[ibody][0] = scale_r * conjqm[ibody][0] + dtf2 * fquat[0]; + conjqm[ibody][1] = scale_r * conjqm[ibody][1] + dtf2 * fquat[1]; + conjqm[ibody][2] = scale_r * conjqm[ibody][2] + dtf2 * fquat[2]; + conjqm[ibody][3] = scale_r * conjqm[ibody][3] + dtf2 * fquat[3]; + } else { + conjqm[ibody][0] += dtf2 * fquat[0]; + conjqm[ibody][1] += dtf2 * fquat[1]; + conjqm[ibody][2] += dtf2 * fquat[2]; + conjqm[ibody][3] += dtf2 * fquat[3]; + } + + MathExtra::invquatvec(quat[ibody],conjqm[ibody],mbody); + MathExtra::matvec(ex_space[ibody],ey_space[ibody],ez_space[ibody], + mbody,angmom[ibody]); + + angmom[ibody][0] *= 0.5; + angmom[ibody][1] *= 0.5; + angmom[ibody][2] *= 0.5; + + MathExtra::angmom_to_omega(angmom[ibody],ex_space[ibody],ey_space[ibody], + ez_space[ibody],inertia[ibody],omega[ibody]); + + if (pstat_flag) { + akin_r += angmom[ibody][0]*omega[ibody][0] + + angmom[ibody][1]*omega[ibody][1] + + angmom[ibody][2]*omega[ibody][2]; + } + } + + // 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 FixRigidNH::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 FixRigidNH::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 FixRigidNH::compute_scalar() +{ + int i,k,ibody; + double kt = boltz * t_target; + double energy,ke_t,ke_q,tmp,Pkq[4]; + + // compute the kinetic parts of H_NVE in Kameraj et al (JCP 2005, pp 224114) + + // translational kinetic energy + + ke_t = 0.0; + for (ibody = 0; ibody < nbody; ibody++) + ke_t += 0.5 * masstotal[ibody] * (vcm[ibody][0]*vcm[ibody][0] + + vcm[ibody][1]*vcm[ibody][1] + + vcm[ibody][2]*vcm[ibody][2]); + + // rotational kinetic energy + + ke_q = 0.0; + for (ibody = 0; ibody < nbody; ibody++) { + for (k = 1; k < 4; k++) { + if (k == 1) { + Pkq[0] = -quat[ibody][1]; + Pkq[1] = quat[ibody][0]; + Pkq[2] = quat[ibody][3]; + Pkq[3] = -quat[ibody][2]; + } else if (k == 2) { + Pkq[0] = -quat[ibody][2]; + Pkq[1] = -quat[ibody][3]; + Pkq[2] = quat[ibody][0]; + Pkq[3] = quat[ibody][1]; + } else if (k == 3) { + Pkq[0] = -quat[ibody][3]; + Pkq[1] = quat[ibody][2]; + Pkq[2] = -quat[ibody][1]; + Pkq[3] = quat[ibody][0]; + } + + tmp = conjqm[ibody][0]*Pkq[0] + conjqm[ibody][1]*Pkq[1] + + conjqm[ibody][2]*Pkq[2] + conjqm[ibody][3]*Pkq[3]; + tmp *= tmp; + + if (fabs(inertia[ibody][k-1]) < 1e-6) tmp = 0.0; + else tmp /= (8.0 * inertia[ibody][k-1]); + ke_q += tmp; + } + } + + 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 FixRigidNH::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 FixRigidNH::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 (nrigid) + 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 (nrigid) + for (i = 0; i< nrigidfix; i++) + modify->fix[rfix[i]]->deform(1); +} + +/* ---------------------------------------------------------------------- + compute target temperature and kinetic energy +-----------------------------------------------------------------------*/ + +void FixRigidNH::compute_temp_target() +{ + double delta = update->ntimestep - update->beginstep; + if (update->endstep > update->beginstep) + delta /= update->endstep - update->beginstep; + else delta = 0.0; + + t_target = t_start + delta * (t_stop-t_start); +} + +/* ---------------------------------------------------------------------- + compute hydrostatic target pressure +-----------------------------------------------------------------------*/ + +void FixRigidNH::compute_press_target() +{ + double delta = update->ntimestep - update->beginstep; + if (update->endstep > update->beginstep) + delta /= update->endstep - update->beginstep; + else delta = 0.0; + + 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; +} + +/* ---------------------------------------------------------------------- + update epsilon_dot +-----------------------------------------------------------------------*/ + +void FixRigidNH::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 FixRigidNH::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 FixRigidNH::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 FixRigidNH::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 FixRigidNH::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 FixRigidNH::reset_target(double t_new) +{ + t_start = t_stop = t_new; +} + +/* ---------------------------------------------------------------------- */ + +void FixRigidNH::allocate_order() +{ + w = new double[t_order]; + wdti1 = new double[t_order]; + wdti2 = new double[t_order]; + wdti4 = new double[t_order]; +} + +/* ---------------------------------------------------------------------- */ + +void FixRigidNH::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 FixRigidNH::deallocate_order() +{ + delete [] w; + delete [] wdti1; + delete [] wdti2; + delete [] wdti4; +} diff --git a/src/fix_rigid_nh.h b/src/fix_rigid_nh.h new file mode 100644 index 0000000000..a348c92aae --- /dev/null +++ b/src/fix_rigid_nh.h @@ -0,0 +1,178 @@ +/* ---------------------------------------------------------------------- + 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_H +#define LMP_FIX_RIGID_NH_H + +#include "fix_rigid.h" + +namespace LAMMPS_NS { + +class FixRigidNH : public FixRigid { + public: + FixRigidNH(class LAMMPS *, int, char **); + virtual ~FixRigidNH(); + 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 **conjqm; // conjugate quaternion momentum + double boltz,nktv2p,mvv2e; // boltzman constant, conversion factors + + 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_current,t_target; + double p_current[3],p_target[3]; + + 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(); + + inline double maclaurin_series(double); +}; + +inline double FixRigidNH::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: 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: Target temperature for fix rigid nvt/npt cannot 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: Fix rigid npt/nph dilate group ID does not exist + +Self-explanatory. + +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: 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: Temperature ID for fix rigid npt/nph does not exist + +Self-explanatory. + +E: Pressure ID for fix rigid npt/nph does not exist + +Self-explanatory. + +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. + +*/ diff --git a/src/fix_rigid_nph.cpp b/src/fix_rigid_nph.cpp new file mode 100644 index 0000000000..49be257e05 --- /dev/null +++ b/src/fix_rigid_nph.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: Tony Sheh (U Michigan), Trung Dac Nguyen (U Michigan) + 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.h" +#include "domain.h" +#include "modify.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +FixRigidNPH::FixRigidNPH(LAMMPS *lmp, int narg, char **arg) : + FixRigidNH(lmp, narg, arg) +{ + // other setting are made by parent + + scalar_flag = 1; + restart_global = 1; + box_change = 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/fix_rigid_nph.h b/src/fix_rigid_nph.h new file mode 100644 index 0000000000..5184952a29 --- /dev/null +++ b/src/fix_rigid_nph.h @@ -0,0 +1,37 @@ +/* ---------------------------------------------------------------------- + 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,FixRigidNPH) + +#else + +#ifndef LMP_FIX_RIGID_NPH_H +#define LMP_FIX_RIGID_NPH_H + +#include "fix_rigid_nh.h" + +namespace LAMMPS_NS { + +class FixRigidNPH : public FixRigidNH { + public: + FixRigidNPH(class LAMMPS *, int, char **); + ~FixRigidNPH() {} +}; + + +} + +#endif +#endif diff --git a/src/fix_rigid_npt.cpp b/src/fix_rigid_npt.cpp new file mode 100644 index 0000000000..235faaba51 --- /dev/null +++ b/src/fix_rigid_npt.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: Tony Sheh (U Michigan), Trung Dac Nguyen (U Michigan) + 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.h" +#include "domain.h" +#include "modify.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +FixRigidNPT::FixRigidNPT(LAMMPS *lmp, int narg, char **arg) : + FixRigidNH(lmp, narg, arg) +{ + // other setting are made by parent + + scalar_flag = 1; + restart_global = 1; + box_change = 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"); + if (t_start <= 0.0 || t_stop <= 0.0) + error->all(FLERR,"Target temperature for fix rigid/npt 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 cannot be 0.0"); + + if (t_period <= 0.0) error->all(FLERR,"Fix rigid/npt 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/fix_rigid_npt.h b/src/fix_rigid_npt.h new file mode 100644 index 0000000000..044a685437 --- /dev/null +++ b/src/fix_rigid_npt.h @@ -0,0 +1,37 @@ +/* ---------------------------------------------------------------------- + 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,FixRigidNPT) + +#else + +#ifndef LMP_FIX_RIGID_NPT_H +#define LMP_FIX_RIGID_NPT_H + +#include "fix_rigid_nh.h" + +namespace LAMMPS_NS { + +class FixRigidNPT : public FixRigidNH { + public: + FixRigidNPT(class LAMMPS *, int, char **); + ~FixRigidNPT() {} +}; + + +} + +#endif +#endif diff --git a/src/fix_rigid_nve.cpp b/src/fix_rigid_nve.cpp index 019fed2876..6cde617e59 100644 --- a/src/fix_rigid_nve.cpp +++ b/src/fix_rigid_nve.cpp @@ -5,7 +5,7 @@ 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 + 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. @@ -14,268 +14,15 @@ /* ---------------------------------------------------------------------- Contributing author: Tony Sheh (U Michigan), Trung Dac Nguyen (U Michigan) references: Kamberaj et al., J. Chem. Phys. 122, 224114 (2005) - Miller et al., J Chem Phys. 116, 8649-8659 + Miller et al., J Chem Phys. 116, 8649-8659 (2002) ------------------------------------------------------------------------- */ -#include "math.h" -#include "stdio.h" -#include "string.h" #include "fix_rigid_nve.h" -#include "math_extra.h" -#include "atom.h" -#include "domain.h" -#include "update.h" -#include "modify.h" -#include "group.h" -#include "comm.h" -#include "force.h" -#include "output.h" -#include "memory.h" -#include "error.h" using namespace LAMMPS_NS; -using namespace FixConst; /* ---------------------------------------------------------------------- */ FixRigidNVE::FixRigidNVE(LAMMPS *lmp, int narg, char **arg) : - FixRigid(lmp, narg, arg) -{ - memory->create(conjqm,nbody,4,"rigid/nve:conjqm"); -} + FixRigidNH(lmp, narg, arg) {} -/* ---------------------------------------------------------------------- */ - -FixRigidNVE::~FixRigidNVE() -{ - memory->destroy(conjqm); -} - -/* ---------------------------------------------------------------------- */ - -void FixRigidNVE::setup(int vflag) -{ - FixRigid::setup(vflag); - - double mbody[3]; - for (int ibody = 0; ibody < nbody; ibody++) { - MathExtra::transpose_matvec(ex_space[ibody],ey_space[ibody],ez_space[ibody], - angmom[ibody],mbody); - MathExtra::quatvec(quat[ibody],mbody,conjqm[ibody]); - conjqm[ibody][0] *= 2.0; - conjqm[ibody][1] *= 2.0; - conjqm[ibody][2] *= 2.0; - conjqm[ibody][3] *= 2.0; - } -} - -/* ---------------------------------------------------------------------- - perform preforce velocity Verlet integration - see Kamberaj paper for step references -------------------------------------------------------------------------- */ - -void FixRigidNVE::initial_integrate(int vflag) -{ - double dtfm,mbody[3],tbody[3],fquat[4]; - double dtf2 = dtf * 2.0; - - for (int ibody = 0; ibody < nbody; ibody++) { - - // step 1.1 - update vcm by 1/2 step - - dtfm = dtf / masstotal[ibody]; - vcm[ibody][0] += dtfm * fcm[ibody][0] * fflag[ibody][0]; - vcm[ibody][1] += dtfm * fcm[ibody][1] * fflag[ibody][1]; - vcm[ibody][2] += dtfm * fcm[ibody][2] * fflag[ibody][2]; - - // step 1.2 - update xcm by full step - - xcm[ibody][0] += dtv * vcm[ibody][0]; - xcm[ibody][1] += dtv * vcm[ibody][1]; - xcm[ibody][2] += dtv * vcm[ibody][2]; - - // step 1.3 - apply torque (body coords) to quaternion momentum - - torque[ibody][0] *= tflag[ibody][0]; - torque[ibody][1] *= tflag[ibody][1]; - torque[ibody][2] *= tflag[ibody][2]; - - MathExtra::transpose_matvec(ex_space[ibody],ey_space[ibody],ez_space[ibody], - torque[ibody],tbody); - MathExtra::quatvec(quat[ibody],tbody,fquat); - - conjqm[ibody][0] += dtf2 * fquat[0]; - conjqm[ibody][1] += dtf2 * fquat[1]; - conjqm[ibody][2] += dtf2 * fquat[2]; - conjqm[ibody][3] += dtf2 * fquat[3]; - - // step 1.4 to 1.13 - use no_squish rotate to update p and q - - no_squish_rotate(3,conjqm[ibody],quat[ibody],inertia[ibody],dtq); - no_squish_rotate(2,conjqm[ibody],quat[ibody],inertia[ibody],dtq); - no_squish_rotate(1,conjqm[ibody],quat[ibody],inertia[ibody],dtv); - no_squish_rotate(2,conjqm[ibody],quat[ibody],inertia[ibody],dtq); - no_squish_rotate(3,conjqm[ibody],quat[ibody],inertia[ibody],dtq); - - // update exyz_space - // transform p back to angmom - // update angular velocity - - MathExtra::q_to_exyz(quat[ibody],ex_space[ibody],ey_space[ibody], - ez_space[ibody]); - MathExtra::invquatvec(quat[ibody],conjqm[ibody],mbody); - MathExtra::matvec(ex_space[ibody],ey_space[ibody],ez_space[ibody], - mbody,angmom[ibody]); - - angmom[ibody][0] *= 0.5; - angmom[ibody][1] *= 0.5; - angmom[ibody][2] *= 0.5; - - MathExtra::angmom_to_omega(angmom[ibody],ex_space[ibody],ey_space[ibody], - ez_space[ibody],inertia[ibody],omega[ibody]); - } - - // virial setup before call to set_xv - - if (vflag) v_setup(vflag); - else evflag = 0; - - // set coords/orient and velocity/rotation of atoms in rigid bodies - // from quarternion and omega - - set_xv(); -} - -/* ---------------------------------------------------------------------- */ - -void FixRigidNVE::final_integrate() -{ - int i,ibody; - double dtfm,xy,xz,yz; - - // sum over atoms to get force and torque on rigid body - - tagint *image = atom->image; - double **x = atom->x; - double **f = atom->f; - int nlocal = atom->nlocal; - - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; - if (triclinic) { - xy = domain->xy; - xz = domain->xz; - yz = domain->yz; - } - - int xbox,ybox,zbox; - double xunwrap,yunwrap,zunwrap,dx,dy,dz; - for (ibody = 0; ibody < nbody; ibody++) - for (i = 0; i < 6; i++) sum[ibody][i] = 0.0; - - for (i = 0; i < nlocal; i++) { - if (body[i] < 0) continue; - ibody = body[i]; - - sum[ibody][0] += f[i][0]; - sum[ibody][1] += f[i][1]; - sum[ibody][2] += f[i][2]; - - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; - - if (triclinic == 0) { - xunwrap = x[i][0] + xbox*xprd; - yunwrap = x[i][1] + ybox*yprd; - zunwrap = x[i][2] + zbox*zprd; - } else { - xunwrap = x[i][0] + xbox*xprd + ybox*xy + zbox*xz; - yunwrap = x[i][1] + ybox*yprd + zbox*yz; - zunwrap = x[i][2] + zbox*zprd; - } - - dx = xunwrap - xcm[ibody][0]; - dy = yunwrap - xcm[ibody][1]; - dz = zunwrap - xcm[ibody][2]; - - sum[ibody][3] += dy*f[i][2] - dz*f[i][1]; - sum[ibody][4] += dz*f[i][0] - dx*f[i][2]; - sum[ibody][5] += dx*f[i][1] - dy*f[i][0]; - } - - // extended particles add their torque to torque of body - - if (extended) { - double **torque_one = atom->torque; - - for (i = 0; i < nlocal; i++) { - if (body[i] < 0) continue; - ibody = body[i]; - - if (eflags[i] & TORQUE) { - sum[ibody][3] += torque_one[i][0]; - sum[ibody][4] += torque_one[i][1]; - sum[ibody][5] += torque_one[i][2]; - } - } - } - - MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world); - - // update vcm and angmom - // include Langevin thermostat forces - // fflag,tflag = 0 for some dimensions in 2d - - double mbody[3],tbody[3],fquat[4]; - double dtf2 = dtf * 2.0; - - for (ibody = 0; ibody < nbody; ibody++) { - fcm[ibody][0] = all[ibody][0] + langextra[ibody][0]; - fcm[ibody][1] = all[ibody][1] + langextra[ibody][1]; - fcm[ibody][2] = all[ibody][2] + langextra[ibody][2]; - torque[ibody][0] = all[ibody][3] + langextra[ibody][3]; - torque[ibody][1] = all[ibody][4] + langextra[ibody][4]; - torque[ibody][2] = all[ibody][5] + langextra[ibody][5]; - - // update vcm by 1/2 step - - dtfm = dtf / masstotal[ibody]; - vcm[ibody][0] += dtfm * fcm[ibody][0] * fflag[ibody][0]; - vcm[ibody][1] += dtfm * fcm[ibody][1] * fflag[ibody][1]; - vcm[ibody][2] += dtfm * fcm[ibody][2] * fflag[ibody][2]; - - // update conjqm, then transform to angmom, set velocity again - // virial is already setup from initial_integrate - - torque[ibody][0] *= tflag[ibody][0]; - torque[ibody][1] *= tflag[ibody][1]; - torque[ibody][2] *= tflag[ibody][2]; - - MathExtra::transpose_matvec(ex_space[ibody],ey_space[ibody],ez_space[ibody], - torque[ibody],tbody); - MathExtra::quatvec(quat[ibody],tbody,fquat); - - conjqm[ibody][0] += dtf2 * fquat[0]; - conjqm[ibody][1] += dtf2 * fquat[1]; - conjqm[ibody][2] += dtf2 * fquat[2]; - conjqm[ibody][3] += dtf2 * fquat[3]; - - MathExtra::invquatvec(quat[ibody],conjqm[ibody],mbody); - MathExtra::matvec(ex_space[ibody],ey_space[ibody],ez_space[ibody], - mbody,angmom[ibody]); - - angmom[ibody][0] *= 0.5; - angmom[ibody][1] *= 0.5; - angmom[ibody][2] *= 0.5; - - MathExtra::angmom_to_omega(angmom[ibody],ex_space[ibody],ey_space[ibody], - ez_space[ibody],inertia[ibody],omega[ibody]); - } - - // set velocity/rotation of atoms in rigid bodies - // virial is already setup from initial_integrate - - set_v(); -} diff --git a/src/fix_rigid_nve.h b/src/fix_rigid_nve.h index ef7105bc83..2f5d854c1b 100644 --- a/src/fix_rigid_nve.h +++ b/src/fix_rigid_nve.h @@ -5,7 +5,7 @@ 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 + 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. @@ -20,26 +20,17 @@ FixStyle(rigid/nve,FixRigidNVE) #ifndef LMP_FIX_RIGID_NVE_H #define LMP_FIX_RIGID_NVE_H -#include "fix_rigid.h" +#include "fix_rigid_nh.h" namespace LAMMPS_NS { -class FixRigidNVE : public FixRigid { +class FixRigidNVE : public FixRigidNH { public: FixRigidNVE(class LAMMPS *, int, char **); - ~FixRigidNVE(); - void setup(int); - void initial_integrate(int); - void final_integrate(); - - private: - double **conjqm; // conjugate quaternion momentum + ~FixRigidNVE() {} }; } #endif #endif -/* ERROR/WARNING messages: - -*/ diff --git a/src/fix_rigid_nvt.cpp b/src/fix_rigid_nvt.cpp index 4598983230..abf681d905 100644 --- a/src/fix_rigid_nvt.cpp +++ b/src/fix_rigid_nvt.cpp @@ -5,7 +5,7 @@ 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 + 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. @@ -14,44 +14,29 @@ /* ---------------------------------------------------------------------- Contributing author: Tony Sheh (U Michigan), Trung Dac Nguyen (U Michigan) references: Kamberaj et al., J. Chem. Phys. 122, 224114 (2005) - Miller et al., J Chem Phys. 116, 8649-8659 + Miller et al., J Chem Phys. 116, 8649-8659 (2002) ------------------------------------------------------------------------- */ -#include "math.h" -#include "stdio.h" -#include "stdlib.h" -#include "string.h" #include "fix_rigid_nvt.h" -#include "math_extra.h" -#include "atom.h" -#include "domain.h" -#include "update.h" -#include "modify.h" -#include "group.h" -#include "comm.h" -#include "force.h" -#include "output.h" -#include "memory.h" #include "error.h" using namespace LAMMPS_NS; -using namespace FixConst; /* ---------------------------------------------------------------------- */ FixRigidNVT::FixRigidNVT(LAMMPS *lmp, int narg, char **arg) : - FixRigid(lmp, narg, arg) -{ - // other settings are made by FixRigid parent + FixRigidNH(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 (tempflag == 0) + if (tstat_flag == 0) error->all(FLERR,"Did not set temp for fix rigid/nvt"); if (t_start < 0.0 || t_stop <= 0.0) error->all(FLERR,"Target temperature for fix rigid/nvt cannot be 0.0"); @@ -60,657 +45,6 @@ FixRigidNVT::FixRigidNVT(LAMMPS *lmp, int narg, char **arg) : 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"); - - allocate_chain(); - allocate_order(); - memory->create(conjqm,nbody,4,"nve_rigid:conjqm"); - - // one-time initialize of thermostat variables - - 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; - } -} - -/* ---------------------------------------------------------------------- */ - -FixRigidNVT::~FixRigidNVT() -{ - deallocate_chain(); - deallocate_order(); - memory->destroy(conjqm); -} - -/* ---------------------------------------------------------------------- */ - -int FixRigidNVT::setmask() -{ - int mask = 0; - mask |= INITIAL_INTEGRATE; - mask |= FINAL_INTEGRATE; - mask |= PRE_NEIGHBOR; - mask |= THERMO_ENERGY; - mask |= INITIAL_INTEGRATE_RESPA; - mask |= FINAL_INTEGRATE_RESPA; - return mask; -} - -/* ---------------------------------------------------------------------- */ - -void FixRigidNVT::init() -{ - FixRigid::init(); - - // 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; - nf_t = nf_r = domain->dimension * nbody; - - for (int ibody = 0; ibody < nbody; ibody++) - for (int k = 0; k < domain->dimension; k++) - if (fabs(inertia[ibody][k]) < 1e-6) nf_r--; - - // see Table 1 in Kamberaj et al - - 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]; - } - - // initialize thermostat settings - - t_target = t_start; - double kt = boltz * t_target; - double 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; - - // initialize thermostat chain positions, velocites, forces - - 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; - f_eta_r[i] = q_r[i-1] * eta_dot_r[i-1] * eta_dot_r[i-1] - kt; - } -} - -/* ---------------------------------------------------------------------- */ - -void FixRigidNVT::setup(int vflag) -{ - FixRigid::setup(vflag); - - t_target = t_start; - - double mbody[3]; - for (int ibody = 0; ibody < nbody; ibody++) { - MathExtra::transpose_matvec(ex_space[ibody],ey_space[ibody],ez_space[ibody], - angmom[ibody],mbody); - MathExtra::quatvec(quat[ibody],mbody,conjqm[ibody]); - conjqm[ibody][0] *= 2.0; - conjqm[ibody][1] *= 2.0; - conjqm[ibody][2] *= 2.0; - conjqm[ibody][3] *= 2.0; - } -} - -/* ---------------------------------------------------------------------- - perform preforce velocity Verlet integration - see Kamberaj paper for step references -------------------------------------------------------------------------- */ - -void FixRigidNVT::initial_integrate(int vflag) -{ - double tmp,akin_t,akin_r,scale_t,scale_r; - double dtfm,mbody[3],tbody[3],fquat[4]; - double dtf2 = dtf * 2.0; - - double delta = update->ntimestep - update->beginstep; - delta /= update->endstep - update->beginstep; - t_target = t_start + delta * (t_stop - t_start); - - // intialize velocity scale for translation and rotation - - akin_t = akin_r = 0.0; - tmp = -1.0 * dtq * eta_dot_t[0]; - scale_t = exp(tmp); - tmp = -1.0 * dtq * eta_dot_r[0]; - scale_r = exp(tmp); - - for (int ibody = 0; ibody < nbody; ibody++) { - - // step 1.1 - update vcm by 1/2 step - - dtfm = dtf / masstotal[ibody]; - vcm[ibody][0] += dtfm * fcm[ibody][0] * fflag[ibody][0]; - vcm[ibody][1] += dtfm * fcm[ibody][1] * fflag[ibody][1]; - vcm[ibody][2] += dtfm * fcm[ibody][2] * fflag[ibody][2]; - vcm[ibody][0] *= scale_t; - vcm[ibody][1] *= scale_t; - vcm[ibody][2] *= scale_t; - - tmp = vcm[ibody][0]*vcm[ibody][0] + vcm[ibody][1]*vcm[ibody][1] + - vcm[ibody][2]*vcm[ibody][2]; - akin_t += masstotal[ibody]*tmp; - - // step 1.2 - update xcm by full step - - xcm[ibody][0] += dtv * vcm[ibody][0]; - xcm[ibody][1] += dtv * vcm[ibody][1]; - xcm[ibody][2] += dtv * vcm[ibody][2]; - - torque[ibody][0] *= tflag[ibody][0]; - torque[ibody][1] *= tflag[ibody][1]; - torque[ibody][2] *= tflag[ibody][2]; - - // step 1.3 - apply torque (body coords) to quaternion momentum - - MathExtra::transpose_matvec(ex_space[ibody],ey_space[ibody],ez_space[ibody], - torque[ibody],tbody); - MathExtra::quatvec(quat[ibody],tbody,fquat); - - conjqm[ibody][0] += dtf2 * fquat[0]; - conjqm[ibody][1] += dtf2 * fquat[1]; - conjqm[ibody][2] += dtf2 * fquat[2]; - conjqm[ibody][3] += dtf2 * fquat[3]; - conjqm[ibody][0] *= scale_r; - conjqm[ibody][1] *= scale_r; - conjqm[ibody][2] *= scale_r; - conjqm[ibody][3] *= scale_r; - - // step 1.4 to 1.8 - use no_squish rotate to update p (i.e. conjqm) and q - - no_squish_rotate(3,conjqm[ibody],quat[ibody],inertia[ibody],dtq); - no_squish_rotate(2,conjqm[ibody],quat[ibody],inertia[ibody],dtq); - no_squish_rotate(1,conjqm[ibody],quat[ibody],inertia[ibody],dtv); - no_squish_rotate(2,conjqm[ibody],quat[ibody],inertia[ibody],dtq); - no_squish_rotate(3,conjqm[ibody],quat[ibody],inertia[ibody],dtq); - - // update the exyz_space from new quaternion - // transform p back to angmom - // update angular velocity - - MathExtra::q_to_exyz(quat[ibody],ex_space[ibody],ey_space[ibody], - ez_space[ibody]); - MathExtra::invquatvec(quat[ibody],conjqm[ibody],mbody); - MathExtra::matvec(ex_space[ibody],ey_space[ibody],ez_space[ibody], - mbody,angmom[ibody]); - - angmom[ibody][0] *= 0.5; - angmom[ibody][1] *= 0.5; - angmom[ibody][2] *= 0.5; - - MathExtra::angmom_to_omega(angmom[ibody],ex_space[ibody],ey_space[ibody], - ez_space[ibody],inertia[ibody],omega[ibody]); - - akin_r += angmom[ibody][0]*omega[ibody][0] + - angmom[ibody][1]*omega[ibody][1] + angmom[ibody][2]*omega[ibody][2]; - } - - // update thermostat chains - - update_nhcp(akin_t,akin_r); - - // virial setup before call to set_xv - - if (vflag) v_setup(vflag); - else evflag = 0; - - // set coords/orient and velocity/rotation of atoms in rigid bodies - // from quarternion and omega - - set_xv(); -} - -/* ---------------------------------------------------------------------- */ - -void FixRigidNVT::final_integrate() -{ - int i,ibody; - double tmp,scale_t,scale_r; - double dtfm,xy,xz,yz; - - // compute velocity scales for translation and rotation - - tmp = -1.0 * dtq * eta_dot_t[0]; - scale_t = exp(tmp); - tmp = -1.0 * dtq * eta_dot_r[0]; - scale_r = exp(tmp); - - // sum over atoms to get force and torque on rigid body - - tagint *image = atom->image; - double **x = atom->x; - double **f = atom->f; - int nlocal = atom->nlocal; - - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; - if (triclinic) { - xy = domain->xy; - xz = domain->xz; - yz = domain->yz; - } - - int xbox,ybox,zbox; - double xunwrap,yunwrap,zunwrap,dx,dy,dz; - for (ibody = 0; ibody < nbody; ibody++) - for (i = 0; i < 6; i++) sum[ibody][i] = 0.0; - - for (i = 0; i < nlocal; i++) { - if (body[i] < 0) continue; - ibody = body[i]; - - sum[ibody][0] += f[i][0]; - sum[ibody][1] += f[i][1]; - sum[ibody][2] += f[i][2]; - - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; - - if (triclinic == 0) { - xunwrap = x[i][0] + xbox*xprd; - yunwrap = x[i][1] + ybox*yprd; - zunwrap = x[i][2] + zbox*zprd; - } else { - xunwrap = x[i][0] + xbox*xprd + ybox*xy + zbox*xz; - yunwrap = x[i][1] + ybox*yprd + zbox*yz; - zunwrap = x[i][2] + zbox*zprd; - } - - dx = xunwrap - xcm[ibody][0]; - dy = yunwrap - xcm[ibody][1]; - dz = zunwrap - xcm[ibody][2]; - - sum[ibody][3] += dy*f[i][2] - dz*f[i][1]; - sum[ibody][4] += dz*f[i][0] - dx*f[i][2]; - sum[ibody][5] += dx*f[i][1] - dy*f[i][0]; - } - - // extended particles add their torque to torque of body - - if (extended) { - double **torque_one = atom->torque; - - for (i = 0; i < nlocal; i++) { - if (body[i] < 0) continue; - ibody = body[i]; - - if (eflags[i] & TORQUE) { - sum[ibody][3] += torque_one[i][0]; - sum[ibody][4] += torque_one[i][1]; - sum[ibody][5] += torque_one[i][2]; - } - } - } - - MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world); - - double mbody[3],tbody[3],fquat[4]; - double dtf2 = dtf * 2.0; - - for (ibody = 0; ibody < nbody; ibody++) { - fcm[ibody][0] = all[ibody][0]; - fcm[ibody][1] = all[ibody][1]; - fcm[ibody][2] = all[ibody][2]; - torque[ibody][0] = all[ibody][3]; - torque[ibody][1] = all[ibody][4]; - torque[ibody][2] = all[ibody][5]; - - // 2.5-2.6 update vcm by 1/2 step - - dtfm = dtf / masstotal[ibody]; - vcm[ibody][0] *= scale_t; - vcm[ibody][1] *= scale_t; - vcm[ibody][2] *= scale_t; - vcm[ibody][0] += dtfm * fcm[ibody][0] * fflag[ibody][0]; - vcm[ibody][1] += dtfm * fcm[ibody][1] * fflag[ibody][1]; - vcm[ibody][2] += dtfm * fcm[ibody][2] * fflag[ibody][2]; - - // 2.1-2.4 update conjqm, angular momentum and angular velocity - // apply body torque flags - - torque[ibody][0] *= tflag[ibody][0]; - torque[ibody][1] *= tflag[ibody][1]; - torque[ibody][2] *= tflag[ibody][2]; - - // convert torque to the body frame - - MathExtra::transpose_matvec(ex_space[ibody],ey_space[ibody],ez_space[ibody], - torque[ibody],tbody); - - // compute "force" for quaternion - - MathExtra::quatvec(quat[ibody],tbody,fquat); - - // update the conjugate quaternion momentum (conjqm) - - conjqm[ibody][0] = scale_r * conjqm[ibody][0] + dtf2 * fquat[0]; - conjqm[ibody][1] = scale_r * conjqm[ibody][1] + dtf2 * fquat[1]; - conjqm[ibody][2] = scale_r * conjqm[ibody][2] + dtf2 * fquat[2]; - conjqm[ibody][3] = scale_r * conjqm[ibody][3] + dtf2 * fquat[3]; - - // compute angular momentum in the body frame - // then convert to the space-fixed frame - - MathExtra::invquatvec(quat[ibody],conjqm[ibody],mbody); - MathExtra::matvec(ex_space[ibody],ey_space[ibody],ez_space[ibody], - mbody,angmom[ibody]); - - angmom[ibody][0] *= 0.5; - angmom[ibody][1] *= 0.5; - angmom[ibody][2] *= 0.5; - - // compute new angular velocity - - MathExtra::angmom_to_omega(angmom[ibody],ex_space[ibody],ey_space[ibody], - ez_space[ibody],inertia[ibody],omega[ibody]); - } - - // set velocity/rotation of atoms in rigid bodies - // virial is already setup from initial_integrate - - set_v(); -} - -/* ---------------------------------------------------------------------- */ - -void FixRigidNVT::update_nhcp(double akin_t, double akin_r) -{ - 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; - - akin_t *= force->mvv2e; - akin_r *= force->mvv2e; - - // 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 order/timestep dependent coefficients - - for (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; - } - - // update force of thermostats coupled to particles - - f_eta_t[0] = (akin_t - gfkt_t) / q_t[0]; - f_eta_r[0] = (akin_r - 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]; - } - } -} - -/* ---------------------------------------------------------------------- - pack entire state of Fix into one write -------------------------------------------------------------------------- */ - -void FixRigidNVT::write_restart(FILE *fp) -{ - int n = 0; - double *list = new double[1+6*t_chain]; - 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++] = f_eta_t[i]; - list[n++] = f_eta_r[i]; - } - - if (comm->me == 0) { - int size = (1 + 6*t_chain)*sizeof(double); - fwrite(&size,sizeof(int),1,fp); - fwrite(list,sizeof(double),1+6*t_chain,fp); - } - - delete list; -} - -/* ---------------------------------------------------------------------- - compute kinetic energy in the extended Hamiltonian - conserved quantity = sum of returned energy and potential energy ------------------------------------------------------------------------*/ - -double FixRigidNVT::compute_scalar() -{ - int i,k,ibody; - double kt = boltz * t_target; - double energy,ke_t,ke_q,tmp,Pkq[4]; - - // compute the kinetic parts of H_NVE in Kameraj et al (JCP 2005, pp 224114) - - // translational kinetic energy - - ke_t = 0.0; - for (ibody = 0; ibody < nbody; ibody++) - ke_t += 0.5 * masstotal[ibody] * (vcm[ibody][0]*vcm[ibody][0] + - vcm[ibody][1]*vcm[ibody][1] + - vcm[ibody][2]*vcm[ibody][2]); - - // rotational kinetic energy - - ke_q = 0.0; - for (ibody = 0; ibody < nbody; ibody++) { - for (k = 1; k < 4; k++) { - if (k == 1) { - Pkq[0] = -quat[ibody][1]; - Pkq[1] = quat[ibody][0]; - Pkq[2] = quat[ibody][3]; - Pkq[3] = -quat[ibody][2]; - } else if (k == 2) { - Pkq[0] = -quat[ibody][2]; - Pkq[1] = -quat[ibody][3]; - Pkq[2] = quat[ibody][0]; - Pkq[3] = quat[ibody][1]; - } else if (k == 3) { - Pkq[0] = -quat[ibody][3]; - Pkq[1] = quat[ibody][2]; - Pkq[2] = -quat[ibody][1]; - Pkq[3] = quat[ibody][0]; - } - - tmp = conjqm[ibody][0]*Pkq[0] + conjqm[ibody][1]*Pkq[1] + - conjqm[ibody][2]*Pkq[2] + conjqm[ibody][3]*Pkq[3]; - tmp *= tmp; - - if (fabs(inertia[ibody][k-1]) < 1e-6) tmp = 0.0; - else tmp /= (8.0 * inertia[ibody][k-1]); - ke_q += tmp; - } - } - - energy = ke_t + ke_q; - - // 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]); - } - - return energy; -} - -/* ---------------------------------------------------------------------- - use state info from restart file to restart the Fix -------------------------------------------------------------------------- */ - -void FixRigidNVT::restart(char *buf) -{ - int n = 0; - double *list = (double *) buf; - - int t_chain_prev = static_cast (list[n++]); - if (t_chain_prev != t_chain) - error->all(FLERR,"Cannot restart fix rigid/nvt with different # of chains"); - - 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++]; - f_eta_t[i] = list[n++]; - f_eta_r[i] = list[n++]; - } -} - -/* ---------------------------------------------------------------------- */ - -void FixRigidNVT::reset_target(double t_new) -{ - t_start = t_stop = t_new; -} - -/* ---------------------------------------------------------------------- */ - -void FixRigidNVT::allocate_chain() -{ - 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]; -} - -/* ---------------------------------------------------------------------- */ - -void FixRigidNVT::allocate_order() -{ - w = new double[t_order]; - wdti1 = new double[t_order]; - wdti2 = new double[t_order]; - wdti4 = new double[t_order]; -} - -/* ---------------------------------------------------------------------- */ - -void FixRigidNVT::deallocate_chain() -{ - 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; -} - -/* ---------------------------------------------------------------------- */ - -void FixRigidNVT::deallocate_order() -{ - delete [] w; - delete [] wdti1; - delete [] wdti2; - delete [] wdti4; + if (t_order != 3 && t_order != 5) + error->all(FLERR,"Fix_modify order must be 3 or 5"); } diff --git a/src/fix_rigid_nvt.h b/src/fix_rigid_nvt.h index 9286e0fa02..337de6d9ed 100644 --- a/src/fix_rigid_nvt.h +++ b/src/fix_rigid_nvt.h @@ -5,7 +5,7 @@ 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 + 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. @@ -20,85 +20,17 @@ FixStyle(rigid/nvt,FixRigidNVT) #ifndef LMP_FIX_RIGID_NVT_H #define LMP_FIX_RIGID_NVT_H -#include "fix_rigid.h" +#include "fix_rigid_nh.h" namespace LAMMPS_NS { -class FixRigidNVT : public FixRigid { +class FixRigidNVT : public FixRigidNH { public: FixRigidNVT(class LAMMPS *, int, char **); - ~FixRigidNVT(); - int setmask(); - void init(); - void setup(int); - void initial_integrate(int); - void final_integrate(); - double compute_scalar(); - void write_restart(FILE *); - void restart(char *); - void reset_target(double); - - private: - double **conjqm; // conjugate quaternion momentum - double boltz; // boltzman constant - double t_target; - - int nf_t,nf_r; // 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 - - void update_nhcp(double, double); - inline double maclaurin_series(double); - - void allocate_chain(); - void allocate_order(); - void deallocate_chain(); - void deallocate_order(); + ~FixRigidNVT() {} }; -inline double FixRigidNVT::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 #endif - -/* ERROR/WARNING messages: - -E: Did not set temp for fix rigid/nvt - -The temp keyword must be used. - -E: Target temperature for fix rigid/nvt cannot be 0.0 - -Self-explanatory. - -E: Fix rigid/nvt period must be > 0.0 - -Self-explanatory - -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: Fix_modify order must be 3 or 5 - -Self-explanatory. - -E: Cannot restart fix rigid/nvt with different # of chains - -This is because the restart file contains per-chain info. - -*/