diff --git a/src/RIGID/fix_rigid.cpp b/src/RIGID/fix_rigid.cpp index 1b0b7402a..d1ce77a29 100644 --- a/src/RIGID/fix_rigid.cpp +++ b/src/RIGID/fix_rigid.cpp @@ -1,2563 +1,2513 @@ /* ---------------------------------------------------------------------- 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. ------------------------------------------------------------------------- */ #include "math.h" #include "stdio.h" #include "stdlib.h" #include "string.h" #include "fix_rigid.h" #include "math_extra.h" #include "atom.h" #include "atom_vec_ellipsoid.h" #include "atom_vec_line.h" #include "atom_vec_tri.h" #include "domain.h" #include "update.h" #include "respa.h" #include "modify.h" #include "group.h" #include "comm.h" #include "random_mars.h" #include "force.h" #include "output.h" #include "math_const.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; 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 #define ATTRIBUTE_PERBODY 11 #define TOLERANCE 1.0e-6 #define EPSILON 1.0e-7 #define SINERTIA 0.4 // moment of inertia prefactor for sphere #define EINERTIA 0.4 // moment of inertia prefactor for ellipsoid #define LINERTIA (1.0/12.0) // moment of inertia prefactor for line segment /* ---------------------------------------------------------------------- */ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) { int i,ibody; scalar_flag = 1; extscalar = 0; time_integrate = 1; rigid_flag = 1; virial_flag = 1; create_attribute = 1; MPI_Comm_rank(world,&me); MPI_Comm_size(world,&nprocs); // perform initial allocation of atom-based arrays // register with Atom class extended = orientflag = dorientflag = 0; body = NULL; displace = NULL; eflags = NULL; orient = NULL; dorient = NULL; grow_arrays(atom->nmax); atom->add_callback(0); // parse args for rigid body specification // set nbody and body[i] for each atom if (narg < 4) error->all(FLERR,"Illegal fix rigid command"); int iarg; mol2body = NULL; body2mol = NULL; // single rigid body // nbody = 1 // all atoms in fix group are part of body if (strcmp(arg[3],"single") == 0) { rstyle = SINGLE; iarg = 4; nbody = 1; int *mask = atom->mask; int nlocal = atom->nlocal; for (i = 0; i < nlocal; i++) { body[i] = -1; if (mask[i] & groupbit) body[i] = 0; } // each molecule in fix group is a rigid body // maxmol = largest molecule ID // ncount = # of atoms in each molecule (have to sum across procs) // nbody = # of non-zero ncount values // use nall as incremented ptr to set body[] values for each atom } else if (strcmp(arg[3],"molecule") == 0) { rstyle = MOLECULE; iarg = 4; if (atom->molecule_flag == 0) error->all(FLERR,"Fix rigid molecule requires atom attribute molecule"); int *mask = atom->mask; tagint *molecule = atom->molecule; int nlocal = atom->nlocal; tagint maxmol_tag = -1; for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) maxmol_tag = MAX(maxmol_tag,molecule[i]); tagint itmp; MPI_Allreduce(&maxmol_tag,&itmp,1,MPI_LMP_TAGINT,MPI_MAX,world); if (itmp+1 > MAXSMALLINT) error->all(FLERR,"Too many molecules for fix rigid"); maxmol = (int) itmp; int *ncount; memory->create(ncount,maxmol+1,"rigid:ncount"); for (i = 0; i <= maxmol; i++) ncount[i] = 0; for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) ncount[molecule[i]]++; memory->create(mol2body,maxmol+1,"rigid:mol2body"); MPI_Allreduce(ncount,mol2body,maxmol+1,MPI_INT,MPI_SUM,world); nbody = 0; for (i = 0; i <= maxmol; i++) if (mol2body[i]) mol2body[i] = nbody++; else mol2body[i] = -1; memory->create(body2mol,nbody,"rigid:body2mol"); nbody = 0; for (i = 0; i <= maxmol; i++) if (mol2body[i] >= 0) body2mol[nbody++] = i; for (i = 0; i < nlocal; i++) { body[i] = -1; if (mask[i] & groupbit) body[i] = mol2body[molecule[i]]; } memory->destroy(ncount); // each listed group is a rigid body // check if all listed groups exist // an atom must belong to fix group and listed group to be in rigid body // error if atom belongs to more than 1 rigid body } else if (strcmp(arg[3],"group") == 0) { if (narg < 5) error->all(FLERR,"Illegal fix rigid command"); rstyle = GROUP; nbody = force->inumeric(FLERR,arg[4]); if (nbody <= 0) error->all(FLERR,"Illegal fix rigid command"); if (narg < 5+nbody) error->all(FLERR,"Illegal fix rigid command"); iarg = 5+nbody; int *igroups = new int[nbody]; for (ibody = 0; ibody < nbody; ibody++) { igroups[ibody] = group->find(arg[5+ibody]); if (igroups[ibody] == -1) error->all(FLERR,"Could not find fix rigid group ID"); } int *mask = atom->mask; int nlocal = atom->nlocal; int flag = 0; for (i = 0; i < nlocal; i++) { body[i] = -1; if (mask[i] & groupbit) for (ibody = 0; ibody < nbody; ibody++) if (mask[i] & group->bitmask[igroups[ibody]]) { if (body[i] >= 0) flag = 1; body[i] = ibody; } } int flagall; MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world); if (flagall) error->all(FLERR,"One or more atoms belong to multiple rigid bodies"); delete [] igroups; } else error->all(FLERR,"Illegal fix rigid command"); // error check on nbody if (nbody == 0) error->all(FLERR,"No rigid bodies defined"); // create all nbody-length arrays memory->create(nrigid,nbody,"rigid:nrigid"); memory->create(masstotal,nbody,"rigid:masstotal"); memory->create(xcm,nbody,3,"rigid:xcm"); memory->create(vcm,nbody,3,"rigid:vcm"); memory->create(fcm,nbody,3,"rigid:fcm"); memory->create(inertia,nbody,3,"rigid:inertia"); memory->create(ex_space,nbody,3,"rigid:ex_space"); memory->create(ey_space,nbody,3,"rigid:ey_space"); memory->create(ez_space,nbody,3,"rigid:ez_space"); memory->create(angmom,nbody,3,"rigid:angmom"); memory->create(omega,nbody,3,"rigid:omega"); memory->create(torque,nbody,3,"rigid:torque"); memory->create(quat,nbody,4,"rigid:quat"); memory->create(imagebody,nbody,"rigid:imagebody"); memory->create(fflag,nbody,3,"rigid:fflag"); memory->create(tflag,nbody,3,"rigid:tflag"); memory->create(langextra,nbody,6,"rigid:langextra"); memory->create(sum,nbody,6,"rigid:sum"); memory->create(all,nbody,6,"rigid:all"); memory->create(remapflag,nbody,4,"rigid:remapflag"); // initialize force/torque flags to default = 1.0 // for 2d: fz, tx, ty = 0.0 array_flag = 1; size_array_rows = nbody; size_array_cols = 15; global_freq = 1; extarray = 0; for (i = 0; i < nbody; i++) { fflag[i][0] = fflag[i][1] = fflag[i][2] = 1.0; tflag[i][0] = tflag[i][1] = tflag[i][2] = 1.0; if (domain->dimension == 2) fflag[i][2] = tflag[i][0] = tflag[i][1] = 0.0; } // parse optional args int seed; langflag = 0; 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"); int mlo,mhi; force->bounds(arg[iarg+1],nbody,mlo,mhi); double xflag,yflag,zflag; if (strcmp(arg[iarg+2],"off") == 0) xflag = 0.0; else if (strcmp(arg[iarg+2],"on") == 0) xflag = 1.0; else error->all(FLERR,"Illegal fix rigid command"); if (strcmp(arg[iarg+3],"off") == 0) yflag = 0.0; else if (strcmp(arg[iarg+3],"on") == 0) yflag = 1.0; else error->all(FLERR,"Illegal fix rigid command"); if (strcmp(arg[iarg+4],"off") == 0) zflag = 0.0; else if (strcmp(arg[iarg+4],"on") == 0) zflag = 1.0; else error->all(FLERR,"Illegal fix rigid command"); if (domain->dimension == 2 && zflag == 1.0) error->all(FLERR,"Fix rigid z force cannot be on for 2d simulation"); int count = 0; for (int m = mlo; m <= mhi; m++) { fflag[m-1][0] = xflag; fflag[m-1][1] = yflag; fflag[m-1][2] = zflag; count++; } if (count == 0) error->all(FLERR,"Illegal fix rigid command"); iarg += 5; } else if (strcmp(arg[iarg],"torque") == 0) { if (iarg+5 > narg) error->all(FLERR,"Illegal fix rigid command"); int mlo,mhi; force->bounds(arg[iarg+1],nbody,mlo,mhi); double xflag,yflag,zflag; if (strcmp(arg[iarg+2],"off") == 0) xflag = 0.0; else if (strcmp(arg[iarg+2],"on") == 0) xflag = 1.0; else error->all(FLERR,"Illegal fix rigid command"); if (strcmp(arg[iarg+3],"off") == 0) yflag = 0.0; else if (strcmp(arg[iarg+3],"on") == 0) yflag = 1.0; else error->all(FLERR,"Illegal fix rigid command"); if (strcmp(arg[iarg+4],"off") == 0) zflag = 0.0; else if (strcmp(arg[iarg+4],"on") == 0) zflag = 1.0; else error->all(FLERR,"Illegal fix rigid command"); if (domain->dimension == 2 && (xflag == 1.0 || yflag == 1.0)) error->all(FLERR,"Fix rigid xy torque cannot be on for 2d simulation"); int count = 0; for (int m = mlo; m <= mhi; m++) { tflag[m-1][0] = xflag; tflag[m-1][1] = yflag; tflag[m-1][2] = zflag; count++; } if (count == 0) error->all(FLERR,"Illegal fix rigid command"); iarg += 5; } else if (strcmp(arg[iarg],"langevin") == 0) { if (iarg+5 > narg) error->all(FLERR,"Illegal fix rigid command"); if (strcmp(style,"rigid") != 0 && strcmp(style,"rigid/nve") != 0) error->all(FLERR,"Illegal fix rigid command"); langflag = 1; t_start = force->numeric(FLERR,arg[iarg+1]); t_stop = force->numeric(FLERR,arg[iarg+2]); t_period = force->numeric(FLERR,arg[iarg+3]); seed = force->inumeric(FLERR,arg[iarg+4]); if (t_period <= 0.0) error->all(FLERR,"Fix rigid langevin period must be > 0.0"); if (seed <= 0) error->all(FLERR,"Illegal fix rigid command"); iarg += 5; } else if (strcmp(arg[iarg],"temp") == 0) { 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"); tstat_flag = 1; t_start = force->numeric(FLERR,arg[iarg+1]); t_stop = force->numeric(FLERR,arg[iarg+2]); t_period = force->numeric(FLERR,arg[iarg+3]); iarg += 4; } else if (strcmp(arg[iarg],"iso") == 0) { if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid command"); 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] = force->numeric(FLERR,arg[iarg+1]); p_stop[0] = p_stop[1] = p_stop[2] = force->numeric(FLERR,arg[iarg+2]); p_period[0] = p_period[1] = p_period[2] = force->numeric(FLERR,arg[iarg+3]); p_flag[0] = p_flag[1] = p_flag[2] = 1; if (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] = force->numeric(FLERR,arg[iarg+1]); p_stop[0] = p_stop[1] = p_stop[2] = force->numeric(FLERR,arg[iarg+2]); p_period[0] = p_period[1] = p_period[2] = force->numeric(FLERR,arg[iarg+3]); p_flag[0] = p_flag[1] = p_flag[2] = 1; if (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] = force->numeric(FLERR,arg[iarg+1]); p_stop[0] = force->numeric(FLERR,arg[iarg+2]); p_period[0] = force->numeric(FLERR,arg[iarg+3]); p_flag[0] = 1; iarg += 4; } else if (strcmp(arg[iarg],"y") == 0) { if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid command"); p_start[1] = force->numeric(FLERR,arg[iarg+1]); p_stop[1] = force->numeric(FLERR,arg[iarg+2]); p_period[1] = force->numeric(FLERR,arg[iarg+3]); p_flag[1] = 1; iarg += 4; } else if (strcmp(arg[iarg],"z") == 0) { if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid command"); p_start[2] = force->numeric(FLERR,arg[iarg+1]); p_stop[2] = force->numeric(FLERR,arg[iarg+2]); p_period[2] = force->numeric(FLERR,arg[iarg+3]); p_flag[2] = 1; iarg += 4; } else if (strcmp(arg[iarg],"couple") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix rigid 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 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 && strcmp(style,"rigid/npt") != 0) error->all(FLERR,"Illegal fix rigid command"); t_chain = force->inumeric(FLERR,arg[iarg+1]); t_iter = force->inumeric(FLERR,arg[iarg+2]); t_order = force->inumeric(FLERR,arg[iarg+3]); iarg += 4; } else if (strcmp(arg[iarg],"pchain") == 0) { if (iarg+2 > 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_chain = force->inumeric(FLERR,arg[iarg+1]); iarg += 2; } else if (strcmp(arg[iarg],"infile") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix rigid command"); delete [] infile; int n = strlen(arg[iarg+1]) + 1; infile = new char[n]; strcpy(infile,arg[iarg+1]); restart_file = 1; iarg += 2; } else error->all(FLERR,"Illegal fix rigid command"); } // 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); else random = NULL; // initialize vector output quantities in case accessed before run for (i = 0; i < nbody; i++) { xcm[i][0] = xcm[i][1] = xcm[i][2] = 0.0; vcm[i][0] = vcm[i][1] = vcm[i][2] = 0.0; fcm[i][0] = fcm[i][1] = fcm[i][2] = 0.0; torque[i][0] = torque[i][1] = torque[i][2] = 0.0; } // nrigid[n] = # of atoms in Nth rigid body // error if one or zero atoms int *ncount = new int[nbody]; for (ibody = 0; ibody < nbody; ibody++) ncount[ibody] = 0; int nlocal = atom->nlocal; for (i = 0; i < nlocal; i++) if (body[i] >= 0) ncount[body[i]]++; MPI_Allreduce(ncount,nrigid,nbody,MPI_INT,MPI_SUM,world); delete [] ncount; for (ibody = 0; ibody < nbody; ibody++) if (nrigid[ibody] <= 1) error->all(FLERR,"One or zero atoms in rigid body"); // bitmasks for properties of extended particles POINT = 1; SPHERE = 2; ELLIPSOID = 4; LINE = 8; TRIANGLE = 16; DIPOLE = 32; OMEGA = 64; ANGMOM = 128; TORQUE = 256; MINUSPI = -MY_PI; TWOPI = 2.0*MY_PI; // if infile set, only setup rigid bodies once, using info from file // this means users cannot change atom properties like mass between runs // if they do, bodies will not reflect the changes staticflag = 0; if (infile) setup_bodies_static(); // print statistics int nsum = 0; for (ibody = 0; ibody < nbody; ibody++) nsum += nrigid[ibody]; if (me == 0) { if (screen) fprintf(screen,"%d rigid bodies with %d atoms\n",nbody,nsum); if (logfile) fprintf(logfile,"%d rigid bodies with %d atoms\n",nbody,nsum); } } /* ---------------------------------------------------------------------- */ FixRigid::~FixRigid() { // unregister callbacks to this fix from Atom class atom->delete_callback(id,0); delete random; delete [] infile; memory->destroy(mol2body); memory->destroy(body2mol); // delete locally stored arrays memory->destroy(body); memory->destroy(displace); memory->destroy(eflags); memory->destroy(orient); memory->destroy(dorient); // delete nbody-length arrays memory->destroy(nrigid); memory->destroy(masstotal); memory->destroy(xcm); memory->destroy(vcm); memory->destroy(fcm); memory->destroy(inertia); memory->destroy(ex_space); memory->destroy(ey_space); memory->destroy(ez_space); memory->destroy(angmom); memory->destroy(omega); memory->destroy(torque); memory->destroy(quat); memory->destroy(imagebody); memory->destroy(fflag); memory->destroy(tflag); memory->destroy(langextra); memory->destroy(sum); memory->destroy(all); memory->destroy(remapflag); } /* ---------------------------------------------------------------------- */ int FixRigid::setmask() { int mask = 0; mask |= INITIAL_INTEGRATE; mask |= FINAL_INTEGRATE; if (langflag) mask |= POST_FORCE; mask |= PRE_NEIGHBOR; mask |= INITIAL_INTEGRATE_RESPA; mask |= FINAL_INTEGRATE_RESPA; return mask; } /* ---------------------------------------------------------------------- */ void FixRigid::init() { int i,ibody; triclinic = domain->triclinic; // atom style pointers to particles that store extra info avec_ellipsoid = (AtomVecEllipsoid *) atom->style_match("ellipsoid"); avec_line = (AtomVecLine *) atom->style_match("line"); avec_tri = (AtomVecTri *) atom->style_match("tri"); // warn if more than one rigid fix int count = 0; for (i = 0; i < modify->nfix; i++) if (strcmp(modify->fix[i]->style,"rigid") == 0) count++; if (count > 1 && me == 0) error->warning(FLERR,"More than one fix rigid"); // error if npt,nph fix comes before rigid fix for (i = 0; i < modify->nfix; i++) { if (strcmp(modify->fix[i]->style,"npt") == 0) break; if (strcmp(modify->fix[i]->style,"nph") == 0) break; } if (i < modify->nfix) { for (int j = i; j < modify->nfix; j++) if (strcmp(modify->fix[j]->style,"rigid") == 0) error->all(FLERR,"Rigid fix must come before NPT/NPH fix"); } // timestep info dtv = update->dt; dtf = 0.5 * update->dt * force->ftm2v; dtq = 0.5 * update->dt; if (strstr(update->integrate_style,"respa")) step_respa = ((Respa *) update->integrate)->step; // setup rigid bodies, using current atom info // allows resetting of atom properties like mass between runs // only do this if not using an infile, else was called in constructor if (!infile) setup_bodies_static(); // temperature scale factor double ndof = 0.0; for (ibody = 0; ibody < nbody; ibody++) { ndof += fflag[ibody][0] + fflag[ibody][1] + fflag[ibody][2]; ndof += tflag[ibody][0] + tflag[ibody][1] + tflag[ibody][2]; } if (ndof > 0.0) tfactor = force->mvv2e / (ndof * force->boltz); else tfactor = 0.0; } /* ---------------------------------------------------------------------- compute initial fcm and torque on bodies, also initial virial reset all particle velocities to be consistent with vcm and omega ------------------------------------------------------------------------- */ void FixRigid::setup(int vflag) { int i,n,ibody; // setup_bodies_dynamic sets vcm and angmom // so angmom_to_omega() and set_v() below will set per-atom vels correctly // re-calling it every run allows reset of body/atom velocities between runs setup_bodies_dynamic(); // fcm = force on center-of-mass of each rigid body double **f = atom->f; int nlocal = atom->nlocal; 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]; } MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world); 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 = torque on each rigid body imageint *image = atom->image; double **x = atom->x; double dx,dy,dz; double unwrap[3]; 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]; domain->unmap(x[i],image[i],unwrap); dx = unwrap[0] - xcm[ibody][0]; dy = unwrap[1] - xcm[ibody][1]; dz = unwrap[2] - xcm[ibody][2]; sum[ibody][0] += dy * f[i][2] - dz * f[i][1]; sum[ibody][1] += dz * f[i][0] - dx * f[i][2]; sum[ibody][2] += 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][0] += torque_one[i][0]; sum[ibody][1] += torque_one[i][1]; sum[ibody][2] += torque_one[i][2]; } } } MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world); for (ibody = 0; ibody < nbody; ibody++) { torque[ibody][0] = all[ibody][0]; torque[ibody][1] = all[ibody][1]; torque[ibody][2] = all[ibody][2]; } // zero langextra in case Langevin thermostat not used // no point to calling post_force() here since langextra // is only added to fcm/torque in final_integrate() for (ibody = 0; ibody < nbody; ibody++) for (i = 0; i < 6; i++) langextra[ibody][i] = 0.0; // virial setup before call to set_v if (vflag) v_setup(vflag); else evflag = 0; // set velocities from angmom & omega for (ibody = 0; ibody < nbody; ibody++) MathExtra::angmom_to_omega(angmom[ibody],ex_space[ibody],ey_space[ibody], ez_space[ibody],inertia[ibody],omega[ibody]); set_v(); // guesstimate virial as 2x the set_v contribution if (vflag_global) for (n = 0; n < 6; n++) virial[n] *= 2.0; if (vflag_atom) { for (i = 0; i < nlocal; i++) for (n = 0; n < 6; n++) vatom[i][n] *= 2.0; } } /* ---------------------------------------------------------------------- */ void FixRigid::initial_integrate(int vflag) { double dtfm; for (int ibody = 0; ibody < nbody; ibody++) { // 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 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]; // update angular momentum by 1/2 step angmom[ibody][0] += dtf * torque[ibody][0] * tflag[ibody][0]; angmom[ibody][1] += dtf * torque[ibody][1] * tflag[ibody][1]; angmom[ibody][2] += dtf * torque[ibody][2] * tflag[ibody][2]; // compute omega at 1/2 step from angmom at 1/2 step and current q // update quaternion a full step via Richardson iteration // returns new normalized quaternion, also updated omega at 1/2 step // update ex,ey,ez to reflect new quaternion MathExtra::angmom_to_omega(angmom[ibody],ex_space[ibody],ey_space[ibody], ez_space[ibody],inertia[ibody],omega[ibody]); MathExtra::richardson(quat[ibody],angmom[ibody],omega[ibody], inertia[ibody],dtq); MathExtra::q_to_exyz(quat[ibody], ex_space[ibody],ey_space[ibody],ez_space[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(); } /* ---------------------------------------------------------------------- apply Langevin thermostat to all 6 DOF of rigid bodies computed by proc 0, broadcast to other procs unlike fix langevin, this stores extra force in extra arrays, which are added in when final_integrate() calculates a new fcm/torque ------------------------------------------------------------------------- */ void FixRigid::post_force(int vflag) { if (me == 0) { double gamma1,gamma2; double delta = update->ntimestep - update->beginstep; if (delta != 0.0) delta /= update->endstep - update->beginstep; t_target = t_start + delta * (t_stop-t_start); double tsqrt = sqrt(t_target); double boltz = force->boltz; double dt = update->dt; double mvv2e = force->mvv2e; double ftm2v = force->ftm2v; for (int i = 0; i < nbody; i++) { gamma1 = -masstotal[i] / t_period / ftm2v; gamma2 = sqrt(masstotal[i]) * tsqrt * sqrt(24.0*boltz/t_period/dt/mvv2e) / ftm2v; langextra[i][0] = gamma1*vcm[i][0] + gamma2*(random->uniform()-0.5); langextra[i][1] = gamma1*vcm[i][1] + gamma2*(random->uniform()-0.5); langextra[i][2] = gamma1*vcm[i][2] + gamma2*(random->uniform()-0.5); gamma1 = -1.0 / t_period / ftm2v; gamma2 = tsqrt * sqrt(24.0*boltz/t_period/dt/mvv2e) / ftm2v; langextra[i][3] = inertia[i][0]*gamma1*omega[i][0] + sqrt(inertia[i][0])*gamma2*(random->uniform()-0.5); langextra[i][4] = inertia[i][1]*gamma1*omega[i][1] + sqrt(inertia[i][1])*gamma2*(random->uniform()-0.5); langextra[i][5] = inertia[i][2]*gamma1*omega[i][2] + sqrt(inertia[i][2])*gamma2*(random->uniform()-0.5); } } MPI_Bcast(&langextra[0][0],6*nbody,MPI_DOUBLE,0,world); } /* ---------------------------------------------------------------------- */ void FixRigid::final_integrate() { int i,ibody; double dtfm; // sum over atoms to get force and torque on rigid body imageint *image = atom->image; double **x = atom->x; double **f = atom->f; int nlocal = atom->nlocal; double dx,dy,dz; double unwrap[3]; 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]; domain->unmap(x[i],image[i],unwrap); dx = unwrap[0] - xcm[ibody][0]; dy = unwrap[1] - xcm[ibody][1]; dz = unwrap[2] - 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]; 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 angular momentum by 1/2 step angmom[ibody][0] += dtf * torque[ibody][0] * tflag[ibody][0]; angmom[ibody][1] += dtf * torque[ibody][1] * tflag[ibody][1]; angmom[ibody][2] += dtf * torque[ibody][2] * tflag[ibody][2]; 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(); } -/* ---------------------------------------------------------------------- - apply evolution operators to quat, quat momentum - see Miller paper cited in fix rigid/nvt and fix rigid/npt -------------------------------------------------------------------------- */ - -void FixRigid::no_squish_rotate(int k, double *p, double *q, - double *inertia, double dt) const -{ - double phi,c_phi,s_phi,kp[4],kq[4]; - - // apply permuation operator on p and q, get kp and kq - - if (k == 1) { - kq[0] = -q[1]; kp[0] = -p[1]; - kq[1] = q[0]; kp[1] = p[0]; - kq[2] = q[3]; kp[2] = p[3]; - kq[3] = -q[2]; kp[3] = -p[2]; - } else if (k == 2) { - kq[0] = -q[2]; kp[0] = -p[2]; - kq[1] = -q[3]; kp[1] = -p[3]; - kq[2] = q[0]; kp[2] = p[0]; - kq[3] = q[1]; kp[3] = p[1]; - } else if (k == 3) { - kq[0] = -q[3]; kp[0] = -p[3]; - kq[1] = q[2]; kp[1] = p[2]; - kq[2] = -q[1]; kp[2] = -p[1]; - kq[3] = q[0]; kp[3] = p[0]; - } - - // obtain phi, cosines and sines - - phi = p[0]*kq[0] + p[1]*kq[1] + p[2]*kq[2] + p[3]*kq[3]; - if (fabs(inertia[k-1]) < 1e-6) phi *= 0.0; - else phi /= 4.0 * inertia[k-1]; - c_phi = cos(dt * phi); - s_phi = sin(dt * phi); - - // advance p and q - - p[0] = c_phi*p[0] + s_phi*kp[0]; - p[1] = c_phi*p[1] + s_phi*kp[1]; - p[2] = c_phi*p[2] + s_phi*kp[2]; - p[3] = c_phi*p[3] + s_phi*kp[3]; - - q[0] = c_phi*q[0] + s_phi*kq[0]; - q[1] = c_phi*q[1] + s_phi*kq[1]; - q[2] = c_phi*q[2] + s_phi*kq[2]; - q[3] = c_phi*q[3] + s_phi*kq[3]; -} - /* ---------------------------------------------------------------------- */ void FixRigid::initial_integrate_respa(int vflag, int ilevel, int iloop) { dtv = step_respa[ilevel]; dtf = 0.5 * step_respa[ilevel] * force->ftm2v; dtq = 0.5 * step_respa[ilevel]; if (ilevel == 0) initial_integrate(vflag); else final_integrate(); } /* ---------------------------------------------------------------------- */ void FixRigid::final_integrate_respa(int ilevel, int iloop) { dtf = 0.5 * step_respa[ilevel] * force->ftm2v; final_integrate(); } /* ---------------------------------------------------------------------- remap xcm of each rigid body back into periodic simulation box done during pre_neighbor so will be after call to pbc() and after fix_deform::pre_exchange() may have flipped box use domain->remap() in case xcm is far away from box due to 1st definition of rigid body or due to box flip if don't do this, then atoms of a body which drifts far away from a triclinic box will be remapped back into box with huge displacements when the box tilt changes via set_x() adjust image flag of body and image flags of all atoms in body ------------------------------------------------------------------------- */ void FixRigid::pre_neighbor() { imageint original,oldimage,newimage; for (int ibody = 0; ibody < nbody; ibody++) { original = imagebody[ibody]; domain->remap(xcm[ibody],imagebody[ibody]); if (original == imagebody[ibody]) remapflag[ibody][3] = 0; else { oldimage = original & IMGMASK; newimage = imagebody[ibody] & IMGMASK; remapflag[ibody][0] = newimage - oldimage; oldimage = (original >> IMGBITS) & IMGMASK; newimage = (imagebody[ibody] >> IMGBITS) & IMGMASK; remapflag[ibody][1] = newimage - oldimage; oldimage = original >> IMG2BITS; newimage = imagebody[ibody] >> IMG2BITS; remapflag[ibody][2] = newimage - oldimage; remapflag[ibody][3] = 1; } } // adjust image flags of any atom in a rigid body whose xcm was remapped // subtracting remapflag = new-old keeps ix,iy,iz near 0 // so body is always in central simulation box imageint *image = atom->image; int nlocal = atom->nlocal; int ibody; imageint idim,otherdims; for (int i = 0; i < nlocal; i++) { if (body[i] == -1) continue; if (remapflag[body[i]][3] == 0) continue; ibody = body[i]; if (remapflag[ibody][0]) { idim = image[i] & IMGMASK; otherdims = image[i] ^ idim; idim -= remapflag[ibody][0]; idim &= IMGMASK; image[i] = otherdims | idim; } if (remapflag[ibody][1]) { idim = (image[i] >> IMGBITS) & IMGMASK; otherdims = image[i] ^ (idim << IMGBITS); idim -= remapflag[ibody][1]; idim &= IMGMASK; image[i] = otherdims | (idim << IMGBITS); } if (remapflag[ibody][2]) { idim = image[i] >> IMG2BITS; otherdims = image[i] ^ (idim << IMG2BITS); idim -= remapflag[ibody][2]; idim &= IMGMASK; image[i] = otherdims | (idim << IMG2BITS); } } } /* ---------------------------------------------------------------------- count # of DOF removed by rigid bodies for atoms in igroup return total count of DOF ------------------------------------------------------------------------- */ int FixRigid::dof(int tgroup) { // cannot count DOF correctly unless setup_bodies_static() has been called if (!staticflag) { if (comm->me == 0) error->warning(FLERR,"Cannot count rigid body degrees-of-freedom " "before bodies are initialized"); return 0; } int tgroupbit = group->bitmask[tgroup]; // nall = # of point particles in each rigid body // mall = # of finite-size particles in each rigid body // particles must also be in temperature group int *mask = atom->mask; int nlocal = atom->nlocal; int *ncount = new int[nbody]; int *mcount = new int[nbody]; for (int ibody = 0; ibody < nbody; ibody++) ncount[ibody] = mcount[ibody] = 0; for (int i = 0; i < nlocal; i++) if (body[i] >= 0 && mask[i] & tgroupbit) { // do not count point particles or point dipoles as extended particles // a spheroid dipole will be counted as extended if (extended && (eflags[i] & ~(POINT | DIPOLE))) mcount[body[i]]++; else ncount[body[i]]++; } int *nall = new int[nbody]; int *mall = new int[nbody]; MPI_Allreduce(ncount,nall,nbody,MPI_INT,MPI_SUM,world); MPI_Allreduce(mcount,mall,nbody,MPI_INT,MPI_SUM,world); // warn if nall+mall != nrigid for any body included in temperature group int flag = 0; for (int ibody = 0; ibody < nbody; ibody++) { if (nall[ibody]+mall[ibody] > 0 && nall[ibody]+mall[ibody] != nrigid[ibody]) flag = 1; } if (flag && me == 0) error->warning(FLERR,"Computing temperature of portions of rigid bodies"); // remove appropriate DOFs for each rigid body wholly in temperature group // N = # of point particles in body // M = # of finite-size particles in body // 3d body has 3N + 6M dof to start with // 2d body has 2N + 3M dof to start with // 3d point-particle body with all non-zero I should have 6 dof, remove 3N-6 // 3d point-particle body (linear) with a 0 I should have 5 dof, remove 3N-5 // 2d point-particle body should have 3 dof, remove 2N-3 // 3d body with any finite-size M should have 6 dof, remove (3N+6M) - 6 // 2d body with any finite-size M should have 3 dof, remove (2N+3M) - 3 int n = 0; if (domain->dimension == 3) { for (int ibody = 0; ibody < nbody; ibody++) if (nall[ibody]+mall[ibody] == nrigid[ibody]) { n += 3*nall[ibody] + 6*mall[ibody] - 6; if (inertia[ibody][0] == 0.0 || inertia[ibody][1] == 0.0 || inertia[ibody][2] == 0.0) n++; } } else if (domain->dimension == 2) { for (int ibody = 0; ibody < nbody; ibody++) if (nall[ibody]+mall[ibody] == nrigid[ibody]) n += 2*nall[ibody] + 3*mall[ibody] - 3; } delete [] ncount; delete [] mcount; delete [] nall; delete [] mall; return n; } /* ---------------------------------------------------------------------- adjust xcm of each rigid body due to box deformation called by various fixes that change box size/shape flag = 0/1 means map from box to lamda coords or vice versa ------------------------------------------------------------------------- */ void FixRigid::deform(int flag) { if (flag == 0) for (int ibody = 0; ibody < nbody; ibody++) domain->x2lamda(xcm[ibody],xcm[ibody]); else for (int ibody = 0; ibody < nbody; ibody++) domain->lamda2x(xcm[ibody],xcm[ibody]); } /* ---------------------------------------------------------------------- set space-frame coords and velocity of each atom in each rigid body set orientation and rotation of extended particles x = Q displace + Xcm, mapped back to periodic box v = Vcm + (W cross (x - Xcm)) ------------------------------------------------------------------------- */ void FixRigid::set_xv() { int ibody; int xbox,ybox,zbox; double x0,x1,x2,v0,v1,v2,fc0,fc1,fc2,massone; double xy,xz,yz; double ione[3],exone[3],eyone[3],ezone[3],vr[6],p[3][3]; imageint *image = atom->image; double **x = atom->x; double **v = atom->v; double **f = atom->f; double *rmass = atom->rmass; double *mass = atom->mass; int *type = atom->type; 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; } // set x and v of each atom for (int i = 0; i < nlocal; i++) { if (body[i] < 0) continue; ibody = body[i]; xbox = (image[i] & IMGMASK) - IMGMAX; ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; zbox = (image[i] >> IMG2BITS) - IMGMAX; // save old positions and velocities for virial if (evflag) { if (triclinic == 0) { x0 = x[i][0] + xbox*xprd; x1 = x[i][1] + ybox*yprd; x2 = x[i][2] + zbox*zprd; } else { x0 = x[i][0] + xbox*xprd + ybox*xy + zbox*xz; x1 = x[i][1] + ybox*yprd + zbox*yz; x2 = x[i][2] + zbox*zprd; } v0 = v[i][0]; v1 = v[i][1]; v2 = v[i][2]; } // x = displacement from center-of-mass, based on body orientation // v = vcm + omega around center-of-mass MathExtra::matvec(ex_space[ibody],ey_space[ibody], ez_space[ibody],displace[i],x[i]); v[i][0] = omega[ibody][1]*x[i][2] - omega[ibody][2]*x[i][1] + vcm[ibody][0]; v[i][1] = omega[ibody][2]*x[i][0] - omega[ibody][0]*x[i][2] + vcm[ibody][1]; v[i][2] = omega[ibody][0]*x[i][1] - omega[ibody][1]*x[i][0] + vcm[ibody][2]; // add center of mass to displacement // map back into periodic box via xbox,ybox,zbox // for triclinic, add in box tilt factors as well if (triclinic == 0) { x[i][0] += xcm[ibody][0] - xbox*xprd; x[i][1] += xcm[ibody][1] - ybox*yprd; x[i][2] += xcm[ibody][2] - zbox*zprd; } else { x[i][0] += xcm[ibody][0] - xbox*xprd - ybox*xy - zbox*xz; x[i][1] += xcm[ibody][1] - ybox*yprd - zbox*yz; x[i][2] += xcm[ibody][2] - zbox*zprd; } // virial = unwrapped coords dotted into body constraint force // body constraint force = implied force due to v change minus f external // assume f does not include forces internal to body // 1/2 factor b/c final_integrate contributes other half // assume per-atom contribution is due to constraint force on that atom if (evflag) { if (rmass) massone = rmass[i]; else massone = mass[type[i]]; fc0 = massone*(v[i][0] - v0)/dtf - f[i][0]; fc1 = massone*(v[i][1] - v1)/dtf - f[i][1]; fc2 = massone*(v[i][2] - v2)/dtf - f[i][2]; vr[0] = 0.5*x0*fc0; vr[1] = 0.5*x1*fc1; vr[2] = 0.5*x2*fc2; vr[3] = 0.5*x0*fc1; vr[4] = 0.5*x0*fc2; vr[5] = 0.5*x1*fc2; v_tally(1,&i,1.0,vr); } } // set orientation, omega, angmom of each extended particle if (extended) { double theta_body,theta; double *shape,*quatatom,*inertiaatom; AtomVecEllipsoid::Bonus *ebonus; if (avec_ellipsoid) ebonus = avec_ellipsoid->bonus; AtomVecLine::Bonus *lbonus; if (avec_line) lbonus = avec_line->bonus; AtomVecTri::Bonus *tbonus; if (avec_tri) tbonus = avec_tri->bonus; double **omega_one = atom->omega; double **angmom_one = atom->angmom; double **mu = atom->mu; int *ellipsoid = atom->ellipsoid; int *line = atom->line; int *tri = atom->tri; for (int i = 0; i < nlocal; i++) { if (body[i] < 0) continue; ibody = body[i]; if (eflags[i] & SPHERE) { omega_one[i][0] = omega[ibody][0]; omega_one[i][1] = omega[ibody][1]; omega_one[i][2] = omega[ibody][2]; } else if (eflags[i] & ELLIPSOID) { shape = ebonus[ellipsoid[i]].shape; quatatom = ebonus[ellipsoid[i]].quat; MathExtra::quatquat(quat[ibody],orient[i],quatatom); MathExtra::qnormalize(quatatom); ione[0] = EINERTIA*rmass[i] * (shape[1]*shape[1] + shape[2]*shape[2]); ione[1] = EINERTIA*rmass[i] * (shape[0]*shape[0] + shape[2]*shape[2]); ione[2] = EINERTIA*rmass[i] * (shape[0]*shape[0] + shape[1]*shape[1]); MathExtra::q_to_exyz(quatatom,exone,eyone,ezone); MathExtra::omega_to_angmom(omega[ibody],exone,eyone,ezone,ione, angmom_one[i]); } else if (eflags[i] & LINE) { if (quat[ibody][3] >= 0.0) theta_body = 2.0*acos(quat[ibody][0]); else theta_body = -2.0*acos(quat[ibody][0]); theta = orient[i][0] + theta_body; while (theta <= MINUSPI) theta += TWOPI; while (theta > MY_PI) theta -= TWOPI; lbonus[line[i]].theta = theta; omega_one[i][0] = omega[ibody][0]; omega_one[i][1] = omega[ibody][1]; omega_one[i][2] = omega[ibody][2]; } else if (eflags[i] & TRIANGLE) { inertiaatom = tbonus[tri[i]].inertia; quatatom = tbonus[tri[i]].quat; MathExtra::quatquat(quat[ibody],orient[i],quatatom); MathExtra::qnormalize(quatatom); MathExtra::q_to_exyz(quatatom,exone,eyone,ezone); MathExtra::omega_to_angmom(omega[ibody],exone,eyone,ezone, inertiaatom,angmom_one[i]); } if (eflags[i] & DIPOLE) { MathExtra::quat_to_mat(quat[ibody],p); MathExtra::matvec(p,dorient[i],mu[i]); MathExtra::snormalize3(mu[i][3],mu[i],mu[i]); } } } } /* ---------------------------------------------------------------------- set space-frame velocity of each atom in a rigid body set omega and angmom of extended particles v = Vcm + (W cross (x - Xcm)) ------------------------------------------------------------------------- */ void FixRigid::set_v() { int xbox,ybox,zbox; double x0,x1,x2,v0,v1,v2,fc0,fc1,fc2,massone; double xy,xz,yz; double ione[3],exone[3],eyone[3],ezone[3],delta[3],vr[6]; double **x = atom->x; double **v = atom->v; double **f = atom->f; double *rmass = atom->rmass; double *mass = atom->mass; int *type = atom->type; imageint *image = atom->image; 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; } // set v of each atom for (int i = 0; i < nlocal; i++) { if (body[i] < 0) continue; const int ibody = body[i]; MathExtra::matvec(ex_space[ibody],ey_space[ibody], ez_space[ibody],displace[i],delta); // save old velocities for virial if (evflag) { v0 = v[i][0]; v1 = v[i][1]; v2 = v[i][2]; } v[i][0] = omega[ibody][1]*delta[2] - omega[ibody][2]*delta[1] + vcm[ibody][0]; v[i][1] = omega[ibody][2]*delta[0] - omega[ibody][0]*delta[2] + vcm[ibody][1]; v[i][2] = omega[ibody][0]*delta[1] - omega[ibody][1]*delta[0] + vcm[ibody][2]; // virial = unwrapped coords dotted into body constraint force // body constraint force = implied force due to v change minus f external // assume f does not include forces internal to body // 1/2 factor b/c initial_integrate contributes other half // assume per-atom contribution is due to constraint force on that atom if (evflag) { if (rmass) massone = rmass[i]; else massone = mass[type[i]]; fc0 = massone*(v[i][0] - v0)/dtf - f[i][0]; fc1 = massone*(v[i][1] - v1)/dtf - f[i][1]; fc2 = massone*(v[i][2] - v2)/dtf - f[i][2]; xbox = (image[i] & IMGMASK) - IMGMAX; ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; zbox = (image[i] >> IMG2BITS) - IMGMAX; if (triclinic == 0) { x0 = x[i][0] + xbox*xprd; x1 = x[i][1] + ybox*yprd; x2 = x[i][2] + zbox*zprd; } else { x0 = x[i][0] + xbox*xprd + ybox*xy + zbox*xz; x1 = x[i][1] + ybox*yprd + zbox*yz; x2 = x[i][2] + zbox*zprd; } vr[0] = 0.5*x0*fc0; vr[1] = 0.5*x1*fc1; vr[2] = 0.5*x2*fc2; vr[3] = 0.5*x0*fc1; vr[4] = 0.5*x0*fc2; vr[5] = 0.5*x1*fc2; v_tally(1,&i,1.0,vr); } } // set omega, angmom of each extended particle if (extended) { double *shape,*quatatom,*inertiaatom; AtomVecEllipsoid::Bonus *ebonus; if (avec_ellipsoid) ebonus = avec_ellipsoid->bonus; AtomVecTri::Bonus *tbonus; if (avec_tri) tbonus = avec_tri->bonus; double **omega_one = atom->omega; double **angmom_one = atom->angmom; int *ellipsoid = atom->ellipsoid; int *tri = atom->tri; for (int i = 0; i < nlocal; i++) { if (body[i] < 0) continue; const int ibody = body[i]; if (eflags[i] & SPHERE) { omega_one[i][0] = omega[ibody][0]; omega_one[i][1] = omega[ibody][1]; omega_one[i][2] = omega[ibody][2]; } else if (eflags[i] & ELLIPSOID) { shape = ebonus[ellipsoid[i]].shape; quatatom = ebonus[ellipsoid[i]].quat; ione[0] = EINERTIA*rmass[i] * (shape[1]*shape[1] + shape[2]*shape[2]); ione[1] = EINERTIA*rmass[i] * (shape[0]*shape[0] + shape[2]*shape[2]); ione[2] = EINERTIA*rmass[i] * (shape[0]*shape[0] + shape[1]*shape[1]); MathExtra::q_to_exyz(quatatom,exone,eyone,ezone); MathExtra::omega_to_angmom(omega[ibody],exone,eyone,ezone,ione, angmom_one[i]); } else if (eflags[i] & LINE) { omega_one[i][0] = omega[ibody][0]; omega_one[i][1] = omega[ibody][1]; omega_one[i][2] = omega[ibody][2]; } else if (eflags[i] & TRIANGLE) { inertiaatom = tbonus[tri[i]].inertia; quatatom = tbonus[tri[i]].quat; MathExtra::q_to_exyz(quatatom,exone,eyone,ezone); MathExtra::omega_to_angmom(omega[ibody],exone,eyone,ezone, inertiaatom,angmom_one[i]); } } } } /* ---------------------------------------------------------------------- initialization of static rigid body attributes called from init() so body/atom properties can be changed between runs unless reading from infile, in which case called once from constructor sets extended flags, masstotal, center-of-mass sets Cartesian and diagonalized inertia tensor sets body image flags, but only on first call ------------------------------------------------------------------------- */ void FixRigid::setup_bodies_static() { int i,ibody; // extended = 1 if any particle in a rigid body is finite size // or has a dipole moment extended = orientflag = dorientflag = 0; AtomVecEllipsoid::Bonus *ebonus; if (avec_ellipsoid) ebonus = avec_ellipsoid->bonus; AtomVecLine::Bonus *lbonus; if (avec_line) lbonus = avec_line->bonus; AtomVecTri::Bonus *tbonus; if (avec_tri) tbonus = avec_tri->bonus; double **mu = atom->mu; double *radius = atom->radius; double *rmass = atom->rmass; double *mass = atom->mass; int *ellipsoid = atom->ellipsoid; int *line = atom->line; int *tri = atom->tri; int *type = atom->type; int nlocal = atom->nlocal; if (atom->radius_flag || atom->ellipsoid_flag || atom->line_flag || atom->tri_flag || atom->mu_flag) { int flag = 0; for (i = 0; i < nlocal; i++) { if (body[i] < 0) continue; if (radius && radius[i] > 0.0) flag = 1; if (ellipsoid && ellipsoid[i] >= 0) flag = 1; if (line && line[i] >= 0) flag = 1; if (tri && tri[i] >= 0) flag = 1; if (mu && mu[i][3] > 0.0) flag = 1; } MPI_Allreduce(&flag,&extended,1,MPI_INT,MPI_MAX,world); } // grow extended arrays and set extended flags for each particle // orientflag = 4 if any particle stores ellipsoid or tri orientation // orientflag = 1 if any particle stores line orientation // dorientflag = 1 if any particle stores dipole orientation if (extended) { if (atom->ellipsoid_flag) orientflag = 4; if (atom->line_flag) orientflag = 1; if (atom->tri_flag) orientflag = 4; if (atom->mu_flag) dorientflag = 1; grow_arrays(atom->nmax); for (i = 0; i < nlocal; i++) { eflags[i] = 0; if (body[i] < 0) continue; // set to POINT or SPHERE or ELLIPSOID or LINE if (radius && radius[i] > 0.0) { eflags[i] |= SPHERE; eflags[i] |= OMEGA; eflags[i] |= TORQUE; } else if (ellipsoid && ellipsoid[i] >= 0) { eflags[i] |= ELLIPSOID; eflags[i] |= ANGMOM; eflags[i] |= TORQUE; } else if (line && line[i] >= 0) { eflags[i] |= LINE; eflags[i] |= OMEGA; eflags[i] |= TORQUE; } else if (tri && tri[i] >= 0) { eflags[i] |= TRIANGLE; eflags[i] |= ANGMOM; eflags[i] |= TORQUE; } else eflags[i] |= POINT; // set DIPOLE if atom->mu and mu[3] > 0.0 if (atom->mu_flag && mu[i][3] > 0.0) eflags[i] |= DIPOLE; } } // compute masstotal & center-of-mass of each rigid body // error if image flag is not 0 in a non-periodic dim double **x = atom->x; imageint *image = atom->image; int *periodicity = domain->periodicity; double xprd = domain->xprd; double yprd = domain->yprd; double zprd = domain->zprd; double xy = domain->xy; double xz = domain->xz; double yz = domain->yz; for (ibody = 0; ibody < nbody; ibody++) for (i = 0; i < 6; i++) sum[ibody][i] = 0.0; int xbox,ybox,zbox; double massone,xunwrap,yunwrap,zunwrap; for (i = 0; i < nlocal; i++) { if (body[i] < 0) continue; ibody = body[i]; xbox = (image[i] & IMGMASK) - IMGMAX; ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; zbox = (image[i] >> IMG2BITS) - IMGMAX; if (rmass) massone = rmass[i]; else massone = mass[type[i]]; if ((xbox && !periodicity[0]) || (ybox && !periodicity[1]) || (zbox && !periodicity[2])) error->one(FLERR,"Fix rigid atom has non-zero image flag " "in a non-periodic dimension"); 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; } sum[ibody][0] += xunwrap * massone; sum[ibody][1] += yunwrap * massone; sum[ibody][2] += zunwrap * massone; sum[ibody][3] += massone; } MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world); for (ibody = 0; ibody < nbody; ibody++) { masstotal[ibody] = all[ibody][3]; xcm[ibody][0] = all[ibody][0]/masstotal[ibody]; xcm[ibody][1] = all[ibody][1]/masstotal[ibody]; xcm[ibody][2] = all[ibody][2]/masstotal[ibody]; } // overwrite masstotal and center-of-mass with file values // inbody[i] = 0/1 if Ith rigid body is initialized by file int *inbody; if (infile) { memory->create(inbody,nbody,"rigid:inbody"); for (ibody = 0; ibody < nbody; ibody++) inbody[ibody] = 0; readfile(0,masstotal,xcm,inbody); } // set image flags for each rigid body to default values // then remap the xcm of each body back into simulation box if needed // staticflag check insures this in only done once, not on successive runs if (!staticflag) { for (ibody = 0; ibody < nbody; ibody++) imagebody[ibody] = ((imageint) IMGMAX << IMG2BITS) | ((imageint) IMGMAX << IMGBITS) | IMGMAX; } pre_neighbor(); // compute 6 moments of inertia of each body in Cartesian reference frame // dx,dy,dz = coords relative to center-of-mass // symmetric 3x3 inertia tensor stored in Voigt notation as 6-vector double 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]; 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]; if (rmass) massone = rmass[i]; else massone = mass[type[i]]; sum[ibody][0] += massone * (dy*dy + dz*dz); sum[ibody][1] += massone * (dx*dx + dz*dz); sum[ibody][2] += massone * (dx*dx + dy*dy); sum[ibody][3] -= massone * dy*dz; sum[ibody][4] -= massone * dx*dz; sum[ibody][5] -= massone * dx*dy; } // extended particles may contribute extra terms to moments of inertia if (extended) { double ivec[6]; double *shape,*quatatom,*inertiaatom; double length,theta; for (i = 0; i < nlocal; i++) { if (body[i] < 0) continue; ibody = body[i]; if (rmass) massone = rmass[i]; else massone = mass[type[i]]; if (eflags[i] & SPHERE) { sum[ibody][0] += SINERTIA*massone * radius[i]*radius[i]; sum[ibody][1] += SINERTIA*massone * radius[i]*radius[i]; sum[ibody][2] += SINERTIA*massone * radius[i]*radius[i]; } else if (eflags[i] & ELLIPSOID) { shape = ebonus[ellipsoid[i]].shape; quatatom = ebonus[ellipsoid[i]].quat; MathExtra::inertia_ellipsoid(shape,quatatom,massone,ivec); sum[ibody][0] += ivec[0]; sum[ibody][1] += ivec[1]; sum[ibody][2] += ivec[2]; sum[ibody][3] += ivec[3]; sum[ibody][4] += ivec[4]; sum[ibody][5] += ivec[5]; } else if (eflags[i] & LINE) { length = lbonus[line[i]].length; theta = lbonus[line[i]].theta; MathExtra::inertia_line(length,theta,massone,ivec); sum[ibody][0] += ivec[0]; sum[ibody][1] += ivec[1]; sum[ibody][2] += ivec[2]; sum[ibody][3] += ivec[3]; sum[ibody][4] += ivec[4]; sum[ibody][5] += ivec[5]; } else if (eflags[i] & TRIANGLE) { inertiaatom = tbonus[tri[i]].inertia; quatatom = tbonus[tri[i]].quat; MathExtra::inertia_triangle(inertiaatom,quatatom,massone,ivec); sum[ibody][0] += ivec[0]; sum[ibody][1] += ivec[1]; sum[ibody][2] += ivec[2]; sum[ibody][3] += ivec[3]; sum[ibody][4] += ivec[4]; sum[ibody][5] += ivec[5]; } } } MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world); // overwrite Cartesian inertia tensor with file values if (infile) readfile(1,NULL,all,inbody); // diagonalize inertia tensor for each body via Jacobi rotations // inertia = 3 eigenvalues = principal moments of inertia // evectors and exzy_space = 3 evectors = principal axes of rigid body int ierror; double cross[3]; double tensor[3][3],evectors[3][3]; for (ibody = 0; ibody < nbody; ibody++) { tensor[0][0] = all[ibody][0]; tensor[1][1] = all[ibody][1]; tensor[2][2] = all[ibody][2]; tensor[1][2] = tensor[2][1] = all[ibody][3]; tensor[0][2] = tensor[2][0] = all[ibody][4]; tensor[0][1] = tensor[1][0] = all[ibody][5]; ierror = MathExtra::jacobi(tensor,inertia[ibody],evectors); if (ierror) error->all(FLERR, "Insufficient Jacobi rotations for rigid body"); ex_space[ibody][0] = evectors[0][0]; ex_space[ibody][1] = evectors[1][0]; ex_space[ibody][2] = evectors[2][0]; ey_space[ibody][0] = evectors[0][1]; ey_space[ibody][1] = evectors[1][1]; ey_space[ibody][2] = evectors[2][1]; ez_space[ibody][0] = evectors[0][2]; ez_space[ibody][1] = evectors[1][2]; ez_space[ibody][2] = evectors[2][2]; // if any principal moment < scaled EPSILON, set to 0.0 double max; max = MAX(inertia[ibody][0],inertia[ibody][1]); max = MAX(max,inertia[ibody][2]); if (inertia[ibody][0] < EPSILON*max) inertia[ibody][0] = 0.0; if (inertia[ibody][1] < EPSILON*max) inertia[ibody][1] = 0.0; if (inertia[ibody][2] < EPSILON*max) inertia[ibody][2] = 0.0; // enforce 3 evectors as a right-handed coordinate system // flip 3rd vector if needed MathExtra::cross3(ex_space[ibody],ey_space[ibody],cross); if (MathExtra::dot3(cross,ez_space[ibody]) < 0.0) MathExtra::negate3(ez_space[ibody]); // create initial quaternion MathExtra::exyz_to_q(ex_space[ibody],ey_space[ibody],ez_space[ibody], quat[ibody]); } // displace = initial atom coords in basis of principal axes // set displace = 0.0 for atoms not in any rigid body // for extended particles, set their orientation wrt to rigid body double qc[4],delta[3]; double *quatatom; double theta_body; for (i = 0; i < nlocal; i++) { if (body[i] < 0) { displace[i][0] = displace[i][1] = displace[i][2] = 0.0; continue; } ibody = body[i]; 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; } delta[0] = xunwrap - xcm[ibody][0]; delta[1] = yunwrap - xcm[ibody][1]; delta[2] = zunwrap - xcm[ibody][2]; MathExtra::transpose_matvec(ex_space[ibody],ey_space[ibody], ez_space[ibody],delta,displace[i]); if (extended) { if (eflags[i] & ELLIPSOID) { quatatom = ebonus[ellipsoid[i]].quat; MathExtra::qconjugate(quat[ibody],qc); MathExtra::quatquat(qc,quatatom,orient[i]); MathExtra::qnormalize(orient[i]); } else if (eflags[i] & LINE) { if (quat[ibody][3] >= 0.0) theta_body = 2.0*acos(quat[ibody][0]); else theta_body = -2.0*acos(quat[ibody][0]); orient[i][0] = lbonus[line[i]].theta - theta_body; while (orient[i][0] <= MINUSPI) orient[i][0] += TWOPI; while (orient[i][0] > MY_PI) orient[i][0] -= TWOPI; if (orientflag == 4) orient[i][1] = orient[i][2] = orient[i][3] = 0.0; } else if (eflags[i] & TRIANGLE) { quatatom = tbonus[tri[i]].quat; MathExtra::qconjugate(quat[ibody],qc); MathExtra::quatquat(qc,quatatom,orient[i]); MathExtra::qnormalize(orient[i]); } else if (orientflag == 4) { orient[i][0] = orient[i][1] = orient[i][2] = orient[i][3] = 0.0; } else if (orientflag == 1) orient[i][0] = 0.0; if (eflags[i] & DIPOLE) { MathExtra::transpose_matvec(ex_space[ibody],ey_space[ibody], ez_space[ibody],mu[i],dorient[i]); MathExtra::snormalize3(mu[i][3],dorient[i],dorient[i]); } else if (dorientflag) dorient[i][0] = dorient[i][1] = dorient[i][2] = 0.0; } } // test for valid principal moments & axes // recompute moments of inertia around new axes // 3 diagonal moments should equal principal moments // 3 off-diagonal moments should be 0.0 // extended particles may contribute extra terms to moments of inertia 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]; if (rmass) massone = rmass[i]; else massone = mass[type[i]]; sum[ibody][0] += massone * (displace[i][1]*displace[i][1] + displace[i][2]*displace[i][2]); sum[ibody][1] += massone * (displace[i][0]*displace[i][0] + displace[i][2]*displace[i][2]); sum[ibody][2] += massone * (displace[i][0]*displace[i][0] + displace[i][1]*displace[i][1]); sum[ibody][3] -= massone * displace[i][1]*displace[i][2]; sum[ibody][4] -= massone * displace[i][0]*displace[i][2]; sum[ibody][5] -= massone * displace[i][0]*displace[i][1]; } if (extended) { double ivec[6]; double *shape,*inertiaatom; double length; for (i = 0; i < nlocal; i++) { if (body[i] < 0) continue; ibody = body[i]; if (rmass) massone = rmass[i]; else massone = mass[type[i]]; if (eflags[i] & SPHERE) { sum[ibody][0] += SINERTIA*massone * radius[i]*radius[i]; sum[ibody][1] += SINERTIA*massone * radius[i]*radius[i]; sum[ibody][2] += SINERTIA*massone * radius[i]*radius[i]; } else if (eflags[i] & ELLIPSOID) { shape = ebonus[ellipsoid[i]].shape; MathExtra::inertia_ellipsoid(shape,orient[i],massone,ivec); sum[ibody][0] += ivec[0]; sum[ibody][1] += ivec[1]; sum[ibody][2] += ivec[2]; sum[ibody][3] += ivec[3]; sum[ibody][4] += ivec[4]; sum[ibody][5] += ivec[5]; } else if (eflags[i] & LINE) { length = lbonus[line[i]].length; MathExtra::inertia_line(length,orient[i][0],massone,ivec); sum[ibody][0] += ivec[0]; sum[ibody][1] += ivec[1]; sum[ibody][2] += ivec[2]; sum[ibody][3] += ivec[3]; sum[ibody][4] += ivec[4]; sum[ibody][5] += ivec[5]; } else if (eflags[i] & TRIANGLE) { inertiaatom = tbonus[tri[i]].inertia; MathExtra::inertia_triangle(inertiaatom,orient[i],massone,ivec); sum[ibody][0] += ivec[0]; sum[ibody][1] += ivec[1]; sum[ibody][2] += ivec[2]; sum[ibody][3] += ivec[3]; sum[ibody][4] += ivec[4]; sum[ibody][5] += ivec[5]; } } } MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world); // error check that re-computed momemts of inertia match diagonalized ones // do not do test for bodies with params read from infile double norm; for (ibody = 0; ibody < nbody; ibody++) { if (infile && inbody[ibody]) continue; if (inertia[ibody][0] == 0.0) { if (fabs(all[ibody][0]) > TOLERANCE) error->all(FLERR,"Fix rigid: Bad principal moments"); } else { if (fabs((all[ibody][0]-inertia[ibody][0])/inertia[ibody][0]) > TOLERANCE) error->all(FLERR,"Fix rigid: Bad principal moments"); } if (inertia[ibody][1] == 0.0) { if (fabs(all[ibody][1]) > TOLERANCE) error->all(FLERR,"Fix rigid: Bad principal moments"); } else { if (fabs((all[ibody][1]-inertia[ibody][1])/inertia[ibody][1]) > TOLERANCE) error->all(FLERR,"Fix rigid: Bad principal moments"); } if (inertia[ibody][2] == 0.0) { if (fabs(all[ibody][2]) > TOLERANCE) error->all(FLERR,"Fix rigid: Bad principal moments"); } else { if (fabs((all[ibody][2]-inertia[ibody][2])/inertia[ibody][2]) > TOLERANCE) error->all(FLERR,"Fix rigid: Bad principal moments"); } norm = (inertia[ibody][0] + inertia[ibody][1] + inertia[ibody][2]) / 3.0; if (fabs(all[ibody][3]/norm) > TOLERANCE || fabs(all[ibody][4]/norm) > TOLERANCE || fabs(all[ibody][5]/norm) > TOLERANCE) error->all(FLERR,"Fix rigid: Bad principal moments"); } if (infile) memory->destroy(inbody); // static properties have now been initialized once // used to prevent re-initialization which would re-read infile staticflag = 1; } /* ---------------------------------------------------------------------- initialization of dynamic rigid body attributes set vcm and angmom, computed explicitly from constituent particles OK if wrong for overlapping particles, since is just setting vcm/angmom at start of run, which can be estimated value ------------------------------------------------------------------------- */ void FixRigid::setup_bodies_dynamic() { int i,ibody; double massone,radone; // vcm = velocity of center-of-mass of each rigid body // angmom = angular momentum of each rigid body double **x = atom->x; double **v = atom->v; double *rmass = atom->rmass; double *mass = atom->mass; imageint *image = atom->image; int *type = atom->type; int nlocal = atom->nlocal; double dx,dy,dz; double unwrap[3]; 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]; if (rmass) massone = rmass[i]; else massone = mass[type[i]]; sum[ibody][0] += v[i][0] * massone; sum[ibody][1] += v[i][1] * massone; sum[ibody][2] += v[i][2] * massone; domain->unmap(x[i],image[i],unwrap); dx = unwrap[0] - xcm[ibody][0]; dy = unwrap[1] - xcm[ibody][1]; dz = unwrap[2] - xcm[ibody][2]; sum[ibody][3] += dy * massone*v[i][2] - dz * massone*v[i][1]; sum[ibody][4] += dz * massone*v[i][0] - dx * massone*v[i][2]; sum[ibody][5] += dx * massone*v[i][1] - dy * massone*v[i][0]; } // extended particles add their rotation to angmom of body if (extended) { AtomVecLine::Bonus *lbonus; if (avec_line) lbonus = avec_line->bonus; double **omega_one = atom->omega; double **angmom_one = atom->angmom; double *radius = atom->radius; int *line = atom->line; for (i = 0; i < nlocal; i++) { if (body[i] < 0) continue; ibody = body[i]; if (eflags[i] & OMEGA) { if (eflags[i] & SPHERE) { radone = radius[i]; sum[ibody][3] += SINERTIA*rmass[i] * radone*radone * omega_one[i][0]; sum[ibody][4] += SINERTIA*rmass[i] * radone*radone * omega_one[i][1]; sum[ibody][5] += SINERTIA*rmass[i] * radone*radone * omega_one[i][2]; } else if (eflags[i] & LINE) { radone = lbonus[line[i]].length; sum[ibody][5] += LINERTIA*rmass[i] * radone*radone * omega_one[i][2]; } } if (eflags[i] & ANGMOM) { sum[ibody][3] += angmom_one[i][0]; sum[ibody][4] += angmom_one[i][1]; sum[ibody][5] += angmom_one[i][2]; } } } MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world); // normalize velocity of COM for (ibody = 0; ibody < nbody; ibody++) { vcm[ibody][0] = all[ibody][0]/masstotal[ibody]; vcm[ibody][1] = all[ibody][1]/masstotal[ibody]; vcm[ibody][2] = all[ibody][2]/masstotal[ibody]; angmom[ibody][0] = all[ibody][3]; angmom[ibody][1] = all[ibody][4]; angmom[ibody][2] = all[ibody][5]; } } /* ---------------------------------------------------------------------- read per rigid body info from user-provided file which = 0 to read total mass and center-of-mass, store in vec and array which = 1 to read 6 moments of inertia, store in array flag inbody = 0 for bodies whose info is read from file nlines = # of lines of rigid body info one line = rigid-ID mass xcm ycm zcm ixx iyy izz ixy ixz iyz ------------------------------------------------------------------------- */ void FixRigid::readfile(int which, double *vec, double **array, int *inbody) { int j,nchunk,id,eofflag; int nlines; FILE *fp; char *eof,*start,*next,*buf; char line[MAXLINE]; if (me == 0) { fp = fopen(infile,"r"); if (fp == NULL) { char str[128]; sprintf(str,"Cannot open fix rigid infile %s",infile); error->one(FLERR,str); } while (1) { eof = fgets(line,MAXLINE,fp); if (eof == NULL) error->one(FLERR,"Unexpected end of fix rigid file"); start = &line[strspn(line," \t\n\v\f\r")]; if (*start != '\0' && *start != '#') break; } sscanf(line,"%d",&nlines); } MPI_Bcast(&nlines,1,MPI_INT,0,world); if (nlines == 0) error->all(FLERR,"Fix rigid file has no lines"); char *buffer = new char[CHUNK*MAXLINE]; char **values = new char*[ATTRIBUTE_PERBODY]; int nread = 0; while (nread < nlines) { nchunk = MIN(nlines-nread,CHUNK); eofflag = comm->read_lines_from_file(fp,nchunk,MAXLINE,buffer); if (eofflag) error->all(FLERR,"Unexpected end of fix rigid file"); buf = buffer; next = strchr(buf,'\n'); *next = '\0'; int nwords = atom->count_words(buf); *next = '\n'; 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 // use ID as-is for SINGLE, as mol-ID for MOLECULE, as-is for GROUP // for which = 0, store mass/com in vec/array // for which = 1, store interia tensor array, invert 3,4,5 values to Voigt 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) error->all(FLERR,"Invalid rigid body ID in fix rigid file"); id = mol2body[id]; } else id--; if (id < 0 || id >= nbody) error->all(FLERR,"Invalid rigid body ID in fix rigid file"); inbody[id] = 1; if (which == 0) { vec[id] = atof(values[1]); array[id][0] = atof(values[2]); array[id][1] = atof(values[3]); array[id][2] = atof(values[4]); } else { array[id][0] = atof(values[5]); array[id][1] = atof(values[6]); array[id][2] = atof(values[7]); array[id][3] = atof(values[10]); array[id][4] = atof(values[9]); array[id][5] = atof(values[8]); } buf = next + 1; } nread += nchunk; } if (me == 0) fclose(fp); delete [] buffer; delete [] values; } /* ---------------------------------------------------------------------- write out restart info for mass, COM, inertia tensor to file identical format to infile option, so info can be read in when restarting only proc 0 writes list of global bodies to file ------------------------------------------------------------------------- */ void FixRigid::write_restart_file(char *file) { if (me) return; char outfile[128]; sprintf(outfile,"%s.rigid",file); FILE *fp = fopen(outfile,"w"); if (fp == NULL) { char str[128]; sprintf(str,"Cannot open fix rigid restart file %s",outfile); error->one(FLERR,str); } fprintf(fp,"# fix rigid mass, COM, inertia tensor info for " "%d bodies on timestep " BIGINT_FORMAT "\n\n", nbody,update->ntimestep); fprintf(fp,"%d\n",nbody); // compute I tensor against xyz axes from diagonalized I and current quat // Ispace = P Idiag P_transpose // P is stored column-wise in exyz_space double p[3][3],pdiag[3][3],ispace[3][3]; int id; for (int i = 0; i < nbody; i++) { if (rstyle == SINGLE || rstyle == GROUP) id = i; else id = body2mol[i]; MathExtra::col2mat(ex_space[i],ey_space[i],ez_space[i],p); MathExtra::times3_diag(p,inertia[i],pdiag); MathExtra::times3_transpose(pdiag,p,ispace); fprintf(fp,"%d %-1.16e %-1.16e %-1.16e %-1.16e " "%-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e\n", id,masstotal[i],xcm[i][0],xcm[i][1],xcm[i][2], ispace[0][0],ispace[1][1],ispace[2][2], ispace[0][1],ispace[0][2],ispace[1][2]); } fclose(fp); } /* ---------------------------------------------------------------------- memory usage of local atom-based arrays ------------------------------------------------------------------------- */ double FixRigid::memory_usage() { int nmax = atom->nmax; double bytes = nmax * sizeof(int); bytes += nmax*3 * sizeof(double); bytes += maxvatom*6 * sizeof(double); // vatom if (extended) { bytes += nmax * sizeof(int); if (orientflag) bytes = nmax*orientflag * sizeof(double); if (dorientflag) bytes = nmax*3 * sizeof(double); } return bytes; } /* ---------------------------------------------------------------------- allocate local atom-based arrays ------------------------------------------------------------------------- */ void FixRigid::grow_arrays(int nmax) { memory->grow(body,nmax,"rigid:body"); memory->grow(displace,nmax,3,"rigid:displace"); if (extended) { memory->grow(eflags,nmax,"rigid:eflags"); if (orientflag) memory->grow(orient,nmax,orientflag,"rigid:orient"); if (dorientflag) memory->grow(dorient,nmax,3,"rigid:dorient"); } } /* ---------------------------------------------------------------------- copy values within local atom-based arrays ------------------------------------------------------------------------- */ void FixRigid::copy_arrays(int i, int j, int delflag) { body[j] = body[i]; displace[j][0] = displace[i][0]; displace[j][1] = displace[i][1]; displace[j][2] = displace[i][2]; if (extended) { eflags[j] = eflags[i]; for (int k = 0; k < orientflag; k++) orient[j][k] = orient[i][k]; if (dorientflag) { dorient[j][0] = dorient[i][0]; dorient[j][1] = dorient[i][1]; dorient[j][2] = dorient[i][2]; } } } /* ---------------------------------------------------------------------- initialize one atom's array values, called when atom is created ------------------------------------------------------------------------- */ void FixRigid::set_arrays(int i) { body[i] = -1; displace[i][0] = 0.0; displace[i][1] = 0.0; displace[i][2] = 0.0; } /* ---------------------------------------------------------------------- pack values in local atom-based arrays for exchange with another proc ------------------------------------------------------------------------- */ int FixRigid::pack_exchange(int i, double *buf) { buf[0] = body[i]; buf[1] = displace[i][0]; buf[2] = displace[i][1]; buf[3] = displace[i][2]; if (!extended) return 4; int m = 4; buf[m++] = eflags[i]; for (int j = 0; j < orientflag; j++) buf[m++] = orient[i][j]; if (dorientflag) { buf[m++] = dorient[i][0]; buf[m++] = dorient[i][1]; buf[m++] = dorient[i][2]; } return m; } /* ---------------------------------------------------------------------- unpack values in local atom-based arrays from exchange with another proc ------------------------------------------------------------------------- */ int FixRigid::unpack_exchange(int nlocal, double *buf) { body[nlocal] = static_cast (buf[0]); displace[nlocal][0] = buf[1]; displace[nlocal][1] = buf[2]; displace[nlocal][2] = buf[3]; if (!extended) return 4; int m = 4; eflags[nlocal] = static_cast (buf[m++]); for (int j = 0; j < orientflag; j++) orient[nlocal][j] = buf[m++]; if (dorientflag) { dorient[nlocal][0] = buf[m++]; dorient[nlocal][1] = buf[m++]; dorient[nlocal][2] = buf[m++]; } return m; } /* ---------------------------------------------------------------------- */ void FixRigid::reset_dt() { dtv = update->dt; dtf = 0.5 * update->dt * force->ftm2v; dtq = 0.5 * update->dt; } /* ---------------------------------------------------------------------- zero linear momentum of each rigid body set Vcm to 0.0, then reset velocities of particles via set_v() ------------------------------------------------------------------------- */ void FixRigid::zero_momentum() { for (int ibody = 0; ibody < nbody; ibody++) vcm[ibody][0] = vcm[ibody][1] = vcm[ibody][2] = 0.0; evflag = 0; set_v(); } /* ---------------------------------------------------------------------- zero angular momentum of each rigid body set angmom/omega to 0.0, then reset velocities of particles via set_v() ------------------------------------------------------------------------- */ void FixRigid::zero_rotation() { for (int ibody = 0; ibody < nbody; ibody++) { angmom[ibody][0] = angmom[ibody][1] = angmom[ibody][2] = 0.0; omega[ibody][0] = omega[ibody][1] = omega[ibody][2] = 0.0; } evflag = 0; set_v(); } /* ---------------------------------------------------------------------- return temperature of collection of rigid bodies non-active DOF are removed by fflag/tflag and in tfactor ------------------------------------------------------------------------- */ double FixRigid::compute_scalar() { double wbody[3],rot[3][3]; double t = 0.0; for (int i = 0; i < nbody; i++) { t += masstotal[i] * (fflag[i][0]*vcm[i][0]*vcm[i][0] + fflag[i][1]*vcm[i][1]*vcm[i][1] + fflag[i][2]*vcm[i][2]*vcm[i][2]); // wbody = angular velocity in body frame MathExtra::quat_to_mat(quat[i],rot); MathExtra::transpose_matvec(rot,angmom[i],wbody); if (inertia[i][0] == 0.0) wbody[0] = 0.0; else wbody[0] /= inertia[i][0]; if (inertia[i][1] == 0.0) wbody[1] = 0.0; else wbody[1] /= inertia[i][1]; if (inertia[i][2] == 0.0) wbody[2] = 0.0; else wbody[2] /= inertia[i][2]; t += tflag[i][0]*inertia[i][0]*wbody[0]*wbody[0] + tflag[i][1]*inertia[i][1]*wbody[1]*wbody[1] + tflag[i][2]*inertia[i][2]*wbody[2]*wbody[2]; } t *= tfactor; return t; } /* ---------------------------------------------------------------------- */ void *FixRigid::extract(const char *str, int &dim) { if (strcmp(str,"body") == 0) { dim = 1; return body; } if (strcmp(str,"masstotal") == 0) { dim = 1; return masstotal; } if (strcmp(str,"t_target") == 0) { dim = 0; return &t_target; } return NULL; } /* ---------------------------------------------------------------------- return translational KE for all rigid bodies KE = 1/2 M Vcm^2 ------------------------------------------------------------------------- */ double FixRigid::extract_ke() { double ke = 0.0; for (int i = 0; i < nbody; i++) ke += masstotal[i] * (vcm[i][0]*vcm[i][0] + vcm[i][1]*vcm[i][1] + vcm[i][2]*vcm[i][2]); return 0.5*ke; } /* ---------------------------------------------------------------------- return rotational KE for all rigid bodies Erotational = 1/2 I wbody^2 ------------------------------------------------------------------------- */ double FixRigid::extract_erotational() { double wbody[3],rot[3][3]; double erotate = 0.0; for (int i = 0; i < nbody; i++) { // wbody = angular velocity in body frame MathExtra::quat_to_mat(quat[i],rot); MathExtra::transpose_matvec(rot,angmom[i],wbody); if (inertia[i][0] == 0.0) wbody[0] = 0.0; else wbody[0] /= inertia[i][0]; if (inertia[i][1] == 0.0) wbody[1] = 0.0; else wbody[1] /= inertia[i][1]; if (inertia[i][2] == 0.0) wbody[2] = 0.0; else wbody[2] /= inertia[i][2]; erotate += inertia[i][0]*wbody[0]*wbody[0] + inertia[i][1]*wbody[1]*wbody[1] + inertia[i][2]*wbody[2]*wbody[2]; } return 0.5*erotate; } /* ---------------------------------------------------------------------- return attributes of a rigid body 15 values per body xcm = 0,1,2; vcm = 3,4,5; fcm = 6,7,8; torque = 9,10,11; image = 12,13,14 ------------------------------------------------------------------------- */ double FixRigid::compute_array(int i, int j) { if (j < 3) return xcm[i][j]; if (j < 6) return vcm[i][j-3]; if (j < 9) return fcm[i][j-6]; if (j < 12) return torque[i][j-9]; if (j == 12) return (imagebody[i] & IMGMASK) - IMGMAX; if (j == 13) return (imagebody[i] >> IMGBITS & IMGMASK) - IMGMAX; return (imagebody[i] >> IMG2BITS) - IMGMAX; } diff --git a/src/RIGID/fix_rigid.h b/src/RIGID/fix_rigid.h index 5467afd81..33ca2df34 100644 --- a/src/RIGID/fix_rigid.h +++ b/src/RIGID/fix_rigid.h @@ -1,261 +1,260 @@ /* -*- c++ -*- ---------------------------------------------------------- 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,FixRigid) #else #ifndef LMP_FIX_RIGID_H #define LMP_FIX_RIGID_H #include "fix.h" namespace LAMMPS_NS { class FixRigid : public Fix { public: FixRigid(class LAMMPS *, int, char **); virtual ~FixRigid(); virtual int setmask(); virtual void init(); virtual void setup(int); virtual void initial_integrate(int); void post_force(int); virtual void final_integrate(); void initial_integrate_respa(int, int, int); void final_integrate_respa(int, int); void write_restart_file(char *); virtual double compute_scalar(); virtual int modify_param(int, char **) {return 0;} double memory_usage(); void grow_arrays(int); void copy_arrays(int, int, int); void set_arrays(int); int pack_exchange(int, double *); int unpack_exchange(int, double *); void pre_neighbor(); int dof(int); void deform(int); void reset_dt(); void zero_momentum(); void zero_rotation(); virtual void *extract(const char*, int &); double extract_ke(); double extract_erotational(); double compute_array(int, int); protected: int me,nprocs; double dtv,dtf,dtq; double *step_respa; int triclinic; double MINUSPI,TWOPI; char *infile; // file to read rigid body attributes from int rstyle; // SINGLE,MOLECULE,GROUP int staticflag; // 1 if static body properties are setup, else 0 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 int *body2mol; // convert rigid body index to mol-ID int maxmol; // size of mol2body = max mol-ID int *body; // which body each atom is part of (-1 if none) double **displace; // displacement of each atom in body coords double *masstotal; // total mass of each rigid body double **xcm; // coords of center-of-mass of each rigid body double **vcm; // velocity of center-of-mass of each double **fcm; // force on center-of-mass of each double **inertia; // 3 principal components of inertia of each double **ex_space,**ey_space,**ez_space; // principal axes of each in space coords double **angmom; // angular momentum of each in space coords double **omega; // angular velocity of each in space coords double **torque; // torque on each rigid body in space coords double **quat; // quaternion of each rigid body imageint *imagebody; // image flags of xcm of each rigid body double **fflag; // flag for on/off of center-of-mass force double **tflag; // flag for on/off of center-of-mass torque double **langextra; // Langevin thermostat forces and torques double **sum,**all; // work vectors for each rigid body int **remapflag; // PBC remap flags for each rigid body int extended; // 1 if any particles have extended attributes int orientflag; // 1 if particles store spatial orientation int dorientflag; // 1 if particles store dipole orientation int *eflags; // flags for extended particles double **orient; // orientation vector of particle wrt rigid body double **dorient; // orientation of dipole mu wrt rigid body double tfactor; // scale factor on temperature of rigid bodies int langflag; // 0/1 = no/yes Langevin thermostat int tstat_flag; // NVT settings double t_start,t_stop,t_target; double t_period,t_freq; int t_chain,t_iter,t_order; 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; class AtomVecTri *avec_tri; int POINT,SPHERE,ELLIPSOID,LINE,TRIANGLE,DIPOLE; // bitmasks for eflags int OMEGA,ANGMOM,TORQUE; - void no_squish_rotate(int, double *, double *, double *, double) const; void set_xv(); void set_v(); void setup_bodies_static(); void setup_bodies_dynamic(); void readfile(int, double *, double **, int *); }; } #endif #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: Fix rigid molecule requires atom attribute molecule Self-explanatory. E: Too many molecules for fix rigid The limit is 2^31 = ~2 billion molecules. E: Could not find fix rigid group ID A group ID used in the fix rigid command does not exist. E: One or more atoms belong to multiple rigid bodies Two or more rigid bodies defined by the fix rigid command cannot contain the same atom. E: No rigid bodies defined The fix specification did not end up defining any rigid bodies. E: Fix rigid z force cannot be on for 2d simulation Self-explanatory. E: Fix rigid xy torque cannot be on for 2d simulation Self-explanatory. E: Fix rigid langevin period must be > 0.0 Self-explanatory. E: Fix rigid npt/nph dilate group ID does not exist Self-explanatory. E: One or zero atoms in rigid body Any rigid body defined by the fix rigid command must contain 2 or more atoms. W: More than one fix rigid It is not efficient to use fix rigid more than once. E: Rigid fix must come before NPT/NPH fix NPT/NPH fix must be defined in input script after all rigid fixes, else the rigid fix contribution to the pressure virial is incorrect. W: Cannot count rigid body degrees-of-freedom before bodies are fully initialized This means the temperature associated with the rigid bodies may be incorrect on this timestep. W: Computing temperature of portions of rigid bodies The group defined by the temperature compute does not encompass all the atoms in one or more rigid bodies, so the change in degrees-of-freedom for the atoms in those partial rigid bodies will not be accounted for. E: Fix rigid atom has non-zero image flag in a non-periodic dimension Image flags for non-periodic dimensions should not be set. E: Insufficient Jacobi rotations for rigid body Eigensolve for rigid body was not sufficiently accurate. E: Fix rigid: Bad principal moments The principal moments of inertia computed for a rigid body are not within the required tolerances. E: Cannot open fix rigid infile %s The specified file cannot be opened. Check that the path and name are correct. E: Unexpected end of fix rigid file A read operation from the file failed. E: Fix rigid file has no lines Self-explanatory. E: Incorrect rigid body format in fix rigid file The number of fields per line is not what expected. E: Invalid rigid body ID in fix rigid file The ID does not match the number of an existing ID of rigid bodies that are defined by the fix rigid command. E: Cannot open fix rigid restart file %s The specified file cannot be opened. Check that the path and name are correct. */ diff --git a/src/RIGID/fix_rigid_nh.cpp b/src/RIGID/fix_rigid_nh.cpp index 4c503ab4a..64cd9a09c 100644 --- a/src/RIGID/fix_rigid_nh.cpp +++ b/src/RIGID/fix_rigid_nh.cpp @@ -1,1388 +1,1391 @@ /* ---------------------------------------------------------------------- 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; +using namespace MathExtra; 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 (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--; - + + nf_t = dimension * nbody; + if (dimension == 3) { + 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--; + } else if (dimension == 2) { + nf_r = nbody; + for (int ibody = 0; ibody < nbody; ibody++) + if (fabs(inertia[ibody][2]) < EPSILON) nf_r--; + } + + g_f = nf_t + nf_r; + onednft = 1.0 + (double)(dimension) / (double)g_f; + onednfr = (double) (dimension) / (double)g_f; + // 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,"Temperature ID for fix rigid nvt/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,"Pressure 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(-dtq * eta_dot_t[0]); scale_t[0] = scale_t[1] = scale_t[2] = tmp; tmp = exp(-dtq * eta_dot_r[0]); scale_r = tmp; } if (pstat_flag) { 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]; } } + + // compute target temperature + // update thermostat chains using akin_t and akin_r + // refer to update_nhcp() in Kamberaj et al. + + if (tstat_flag) { + compute_temp_target(); + nhc_temp_integrate(); + } + // update thermostat chains coupled with barostat + // refer to update_nhcb() in Kamberaj et al. + + if (pstat_flag) { + nhc_press_integrate(); + } + // 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)); + // reset akin_t and akin_r, to be accumulated for use in nh_epsilon_dot() + akin_t = akin_r = 0.0; } // sum over atoms to get force and torque on rigid body imageint *image = atom->image; double **x = atom->x; double **f = atom->f; int nlocal = atom->nlocal; double 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 - + + // compute current temperature if (tcomputeflag) t_current = temperature->compute_scalar(); + + // compute current and target pressures + // update epsilon dot using akin_t and akin_r + if (pstat_flag) { - if (pstyle == ISO) pressure->compute_scalar(); - else pressure->compute_vector(); + if (pstyle == ISO) { + temperature->compute_scalar(); + pressure->compute_scalar(); + } else { + temperature->compute_vector(); + 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(); + compute_press_target(); + + nh_epsilon_dot(); + } } /* ---------------------------------------------------------------------- */ 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]; + 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; + int i,j,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; + 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-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]; } + kecurrent /= pdim; f_eta_b[0] = (kecurrent - lkt_press) / q_b[0]; + + // multiple timestep iteration - // update thermostat velocities a half step + for (i = 0; i < t_iter; i++) { + for (j = 0; j < t_order; j++) { + + // update thermostat velocities a half step - eta_dot_b[p_chain-1] += 0.5 * dtq * f_eta_b[p_chain-1]; + eta_dot_b[p_chain-1] += wdti2[j] * 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; - } + for (k = 1; k < p_chain; k++) { + tmp = wdti4[j] * eta_dot_b[p_chain-k]; + ms = maclaurin_series(tmp); + s = exp(-0.5 * tmp); + s2 = s * s; + eta_dot_b[p_chain-k-1] = eta_dot_b[p_chain-k-1] * s2 + + wdti2[j] * f_eta_b[p_chain-k-1] * s * ms; + } - // update thermostat positions + // update thermostat positions - for (k = 0; k < p_chain; k++) - eta_b[k] += dtv * eta_dot_b[k]; + for (k = 0; k < p_chain; k++) + eta_b[k] += wdti1[j] * eta_dot_b[k]; - // update epsilon dot + // update thermostat forces - 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]; + for (k = 1; k < p_chain; k++) { + f_eta_b[k] = q_b[k-1] * eta_dot_b[k-1] * eta_dot_b[k-1] - kt; + f_eta_b[k] /= q_b[k]; + } + + // update thermostat velocites a full step - // update thermostat velocites a full step + for (k = 0; k < p_chain-1; k++) { + tmp = wdti4[j] * 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 + wdti2[j] * 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]; + } - 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] += wdti2[j] * f_eta_b[p_chain-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() { const double kt = boltz * t_target; double energy; int i; energy = FixRigid::compute_scalar(); 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 + double e = 0.0; for (i = 0; i < 3; i++) - energy += 0.5 * epsilon_mass[i] * epsilon_dot[i] * epsilon_dot[i]; + if (p_flag[i]) + e += epsilon_mass[i] * epsilon_dot[i] * epsilon_dot[i]; + energy += e*(0.5/pdim); 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) + if (nrigidfix) for (i = 0; i < nrigidfix; i++) modify->fix[rfix[i]]->deform(0); // reset global and local box to new size/shape for (i = 0; i < 3; i++) { if (p_flag[i]) { oldlo = domain->boxlo[i]; oldhi = domain->boxhi[i]; ctr = 0.5 * (oldlo + oldhi); expfac = exp(dtq * epsilon_dot[i]); domain->boxlo[i] = (oldlo-ctr)*expfac + ctr; domain->boxhi[i] = (oldhi-ctr)*expfac + ctr; } } domain->set_global_box(); domain->set_local_box(); // convert pertinent atoms and rigid bodies back to box coords if (allremap) domain->lamda2x(nlocal); else { for (i = 0; i < nlocal; i++) if (mask[i] & dilate_group_bit) domain->lamda2x(x[i],x[i]); } - if (nrigid) + if (nrigidfix) 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 (delta != 0.0) delta /= update->endstep - update->beginstep; 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 (delta != 0.0) delta /= update->endstep - update->beginstep; p_hydro = 0.0; for (int i = 0; i < 3; i++) if (p_flag[i]) { p_target[i] = p_start[i] + delta * (p_stop[i]-p_start[i]); p_hydro += p_target[i]; } p_hydro /= pdim; } /* ---------------------------------------------------------------------- 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/RIGID/fix_rigid_nh_small.cpp b/src/RIGID/fix_rigid_nh_small.cpp index c050f2de9..de9b5a881 100644 --- a/src/RIGID/fix_rigid_nh_small.cpp +++ b/src/RIGID/fix_rigid_nh_small.cpp @@ -1,1532 +1,1484 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- Contributing author: Trung Dac Nguyen (ORNL) references: Kamberaj et al., J. Chem. Phys. 122, 224114 (2005) Miller et al., J Chem Phys. 116, 8649-8659 (2002) ------------------------------------------------------------------------- */ #include "math.h" #include "stdio.h" #include "string.h" #include "fix_rigid_nh_small.h" #include "math_extra.h" #include "atom.h" #include "compute.h" #include "domain.h" #include "update.h" #include "modify.h" #include "fix_deform.h" #include "group.h" #include "comm.h" #include "force.h" #include "kspace.h" #include "output.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; using namespace FixConst; +using namespace MathExtra; enum{NONE,XYZ,XY,YZ,XZ}; // same as in FixRigid enum{ISO,ANISO,TRICLINIC}; // same as in FixRigid #define EPSILON 1.0e-7 enum{FULL_BODY,INITIAL,FINAL,FORCE_TORQUE,VCM_ANGMOM,XCM_MASS,ITENSOR,DOF}; /* ---------------------------------------------------------------------- */ FixRigidNHSmall::FixRigidNHSmall(LAMMPS *lmp, int narg, char **arg) : FixRigidSmall(lmp, narg, arg) { // error checks if ((p_flag[0] == 1 && p_period[0] <= 0.0) || (p_flag[1] == 1 && p_period[1] <= 0.0) || (p_flag[2] == 1 && p_period[2] <= 0.0)) error->all(FLERR,"Fix rigid/small npt/nph period must be > 0.0"); if (domain->dimension == 2 && p_flag[2]) error->all(FLERR,"Invalid fix rigid/small npt/nph command for a 2d simulation"); if (domain->dimension == 2 && (pcouple == YZ || pcouple == XZ)) error->all(FLERR,"Invalid fix rigid/small npt/nph command for a 2d simulation"); if (pcouple == XYZ && (p_flag[0] == 0 || p_flag[1] == 0)) error->all(FLERR,"Invalid fix rigid/small npt/nph command pressure settings"); if (pcouple == XYZ && domain->dimension == 3 && p_flag[2] == 0) error->all(FLERR,"Invalid fix rigid/small npt/nph command pressure settings"); if (pcouple == XY && (p_flag[0] == 0 || p_flag[1] == 0)) error->all(FLERR,"Invalid fix rigid/small npt/nph command pressure settings"); if (pcouple == YZ && (p_flag[1] == 0 || p_flag[2] == 0)) error->all(FLERR,"Invalid fix rigid/small npt/nph command pressure settings"); if (pcouple == XZ && (p_flag[0] == 0 || p_flag[2] == 0)) error->all(FLERR,"Invalid fix rigid/small npt/nph command pressure settings"); // require periodicity in tensile dimension if (p_flag[0] && domain->xperiodic == 0) error->all(FLERR, "Cannot use fix rigid/small npt/nph on a non-periodic dimension"); if (p_flag[1] && domain->yperiodic == 0) error->all(FLERR, "Cannot use fix rigid/small npt/nph on a non-periodic dimension"); if (p_flag[2] && domain->zperiodic == 0) error->all(FLERR, "Cannot use fix rigid/small npt/nph on a non-periodic dimension"); if (pcouple == XYZ && domain->dimension == 3 && (p_start[0] != p_start[1] || p_start[0] != p_start[2] || p_stop[0] != p_stop[1] || p_stop[0] != p_stop[2] || p_period[0] != p_period[1] || p_period[0] != p_period[2])) error->all(FLERR,"Invalid fix rigid/small npt/nph pressure settings"); if (pcouple == XYZ && domain->dimension == 2 && (p_start[0] != p_start[1] || p_stop[0] != p_stop[1] || p_period[0] != p_period[1])) error->all(FLERR,"Invalid fix rigid/small npt/nph pressure settings"); if (pcouple == XY && (p_start[0] != p_start[1] || p_stop[0] != p_stop[1] || p_period[0] != p_period[1])) error->all(FLERR,"Invalid fix rigid/small npt/nph pressure settings"); if (pcouple == YZ && (p_start[1] != p_start[2] || p_stop[1] != p_stop[2] || p_period[1] != p_period[2])) error->all(FLERR,"Invalid fix rigid/small npt/nph pressure settings"); if (pcouple == XZ && (p_start[0] != p_start[2] || p_stop[0] != p_stop[2] || p_period[0] != p_period[2])) error->all(FLERR,"Invalid fix rigid/small npt/nph pressure settings"); if ((tstat_flag && t_period <= 0.0) || (p_flag[0] && p_period[0] <= 0.0) || (p_flag[1] && p_period[1] <= 0.0) || (p_flag[2] && p_period[2] <= 0.0)) error->all(FLERR,"Fix rigid/small nvt/npt/nph damping parameters must be > 0.0"); // memory allocation and initialization if (tstat_flag || pstat_flag) { allocate_chain(); allocate_order(); } if (tstat_flag) { eta_t[0] = eta_r[0] = 0.0; eta_dot_t[0] = eta_dot_r[0] = 0.0; f_eta_t[0] = f_eta_r[0] = 0.0; for (int i = 1; i < t_chain; i++) { eta_t[i] = eta_r[i] = 0.0; eta_dot_t[i] = eta_dot_r[i] = 0.0; } } if (pstat_flag) { epsilon_dot[0] = epsilon_dot[1] = epsilon_dot[2] = 0.0; eta_b[0] = eta_dot_b[0] = f_eta_b[0] = 0.0; for (int i = 1; i < p_chain; i++) eta_b[i] = eta_dot_b[i] = 0.0; } // rigid body pointers nrigidfix = 0; rfix = NULL; vol0 = 0.0; t0 = 1.0; tcomputeflag = 0; pcomputeflag = 0; } /* ---------------------------------------------------------------------- */ FixRigidNHSmall::~FixRigidNHSmall() { if (tstat_flag || pstat_flag) { deallocate_chain(); deallocate_order(); } if (rfix) delete [] rfix; if (tcomputeflag) { modify->delete_compute(id_temp); delete [] id_temp; } // delete pressure if fix created it if (pstat_flag) { if (pcomputeflag) modify->delete_compute(id_press); delete [] id_press; } } /* ---------------------------------------------------------------------- */ int FixRigidNHSmall::setmask() { int mask = 0; mask = FixRigidSmall::setmask(); if (tstat_flag || pstat_flag) mask |= THERMO_ENERGY; return mask; } /* ---------------------------------------------------------------------- */ void FixRigidNHSmall::init() { FixRigidSmall::init(); // recheck that dilate group has not been deleted if (allremap == 0) { int idilate = group->find(id_dilate); if (idilate == -1) error->all(FLERR,"Fix rigid npt/nph dilate group ID does not exist"); dilate_group_bit = group->bitmask[idilate]; } // initialize thermostats // set timesteps, constants // store Yoshida-Suzuki integrator parameters dtv = update->dt; dtf = 0.5 * update->dt * force->ftm2v; dtq = 0.5 * update->dt; boltz = force->boltz; nktv2p = force->nktv2p; mvv2e = force->mvv2e; dimension = domain->dimension; if (force->kspace) kspace_flag = 1; else kspace_flag = 0; - + // see Table 1 in Kamberaj et al if (tstat_flag || pstat_flag) { if (t_order == 3) { w[0] = 1.0 / (2.0 - pow(2.0, 1.0/3.0)); w[1] = 1.0 - 2.0*w[0]; w[2] = w[0]; } else if (t_order == 5) { w[0] = 1.0 / (4.0 - pow(4.0, 1.0/3.0)); w[1] = w[0]; w[2] = 1.0 - 4.0 * w[0]; w[3] = w[0]; w[4] = w[0]; } - } + } int icompute; if (tcomputeflag) { icompute = modify->find_compute(id_temp); if (icompute < 0) - error->all(FLERR,"Temp ID for fix rigid npt/nph does not exist"); + error->all(FLERR,"Temperature ID for fix rigid nvt/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"); + 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"); + error->all(FLERR,"Pressure ID for fix rigid npt/nph does not exist"); pressure = modify->compute[icompute]; // detect if any rigid fixes exist so rigid bodies move on remap // rfix[] = indices to each fix rigid // this will include self if (rfix) delete [] rfix; nrigidfix = 0; rfix = NULL; for (int i = 0; i < modify->nfix; i++) if (modify->fix[i]->rigid_flag) nrigidfix++; if (nrigidfix) { rfix = new int[nrigidfix]; nrigidfix = 0; for (int i = 0; i < modify->nfix; i++) if (modify->fix[i]->rigid_flag) rfix[nrigidfix++] = i; } } } /* ---------------------------------------------------------------------- */ void FixRigidNHSmall::setup(int vflag) { FixRigidSmall::setup(vflag); // total translational and rotational degrees of freedom - int k,ibody; - - nf_t = nf_r = dimension * nlocal_body; - for (ibody = 0; ibody < nlocal_body; ibody++) { - for (k = 0; k < domain->dimension; k++) - if (fabs(body[ibody].inertia[k]) < EPSILON) nf_r--; + nf_t = dimension * nlocal_body; + if (dimension == 3) { + nf_r = dimension * nlocal_body; + for (int ibody = 0; ibody < nlocal_body; ibody++) { + Body *b = &body[ibody]; + for (int k = 0; k < domain->dimension; k++) + if (fabs(b->inertia[k]) < EPSILON) nf_r--; + } + } else if (dimension == 2) { + nf_r = nlocal_body; + for (int ibody = 0; ibody < nlocal_body; ibody++) { + Body *b = &body[ibody]; + if (fabs(b->inertia[2]) < EPSILON) nf_r--; + } } double nf[2], nfall[2]; nf[0] = nf_t; nf[1] = nf_r; MPI_Allreduce(nf,nfall,2,MPI_DOUBLE,MPI_SUM,world); nf_t = nfall[0]; nf_r = nfall[1]; g_f = nf_t + nf_r; onednft = 1.0 + (double)(dimension) / (double)g_f; onednfr = (double) (dimension) / (double)g_f; double mbody[3]; akin_t = akin_r = 0.0; for (int ibody = 0; ibody < nlocal_body; ibody++) { Body *b = &body[ibody]; MathExtra::transpose_matvec(b->ex_space,b->ey_space,b->ez_space, b->angmom,mbody); MathExtra::quatvec(b->quat,mbody,b->conjqm); b->conjqm[0] *= 2.0; b->conjqm[1] *= 2.0; b->conjqm[2] *= 2.0; b->conjqm[3] *= 2.0; if (tstat_flag || pstat_flag) { akin_t += b->mass*(b->vcm[0]*b->vcm[0] + b->vcm[1]*b->vcm[1] + b->vcm[2]*b->vcm[2]); akin_r += b->angmom[0]*b->omega[0] + b->angmom[1]*b->omega[1] + b->angmom[2]*b->omega[2]; } } // accumulate translational and rotational kinetic energies if (tstat_flag || pstat_flag) { double ke[2],keall[2]; ke[0] = akin_t; ke[1] = akin_r; MPI_Allreduce(ke,keall,2,MPI_DOUBLE,MPI_SUM,world); akin_t = keall[0]; akin_r = keall[1]; } // compute target temperature if (tstat_flag) compute_temp_target(); else if (pstat_flag) { t0 = temperature->compute_scalar(); if (t0 == 0.0) { if (strcmp(update->unit_style,"lj") == 0) t0 = 1.0; else t0 = 300.0; } t_target = t0; } // compute target pressure // compute current pressure // trigger virial computation on next timestep if (pstat_flag) { compute_press_target(); temperature->compute_scalar(); if (pstyle == ISO) pressure->compute_scalar(); else pressure->compute_vector(); couple(); pressure->addstep(update->ntimestep+1); } // initialize thermostat/barostat settings double kt, t_mass, tb_mass; kt = boltz * t_target; if (tstat_flag) { t_mass = kt / (t_freq*t_freq); q_t[0] = nf_t * t_mass; q_r[0] = nf_r * t_mass; for (int i = 1; i < t_chain; i++) q_t[i] = q_r[i] = t_mass; for (int i = 1; i < t_chain; i++) { f_eta_t[i] = (q_t[i-1] * eta_dot_t[i-1] * eta_dot_t[i-1] - kt)/q_t[i]; f_eta_r[i] = (q_r[i-1] * eta_dot_r[i-1] * eta_dot_r[i-1] - kt)/q_r[i]; } } // initial forces on barostat thermostat variables if (pstat_flag) { for (int i = 0; i < 3; i++) if (p_flag[i]) { epsilon_mass[i] = (g_f + dimension) * kt / (p_freq[i]*p_freq[i]); epsilon[i] = log(vol0)/dimension; } tb_mass = kt / (p_freq_max * p_freq_max); q_b[0] = dimension * dimension * tb_mass; for (int i = 1; i < p_chain; i++) { q_b[i] = tb_mass; f_eta_b[i] = (q_b[i] * eta_dot_b[i-1] * eta_dot_b[i-1] - kt)/q_b[i]; } } // update order/timestep dependent coefficients if (tstat_flag || pstat_flag) { for (int i = 0; i < t_order; i++) { wdti1[i] = w[i] * dtv / t_iter; wdti2[i] = wdti1[i] / 2.0; wdti4[i] = wdti1[i] / 4.0; } } } /* ---------------------------------------------------------------------- perform preforce velocity Verlet integration see Kamberaj paper for step references ------------------------------------------------------------------------- */ void FixRigidNHSmall::initial_integrate(int vflag) { double tmp,scale_r,scale_t[3],scale_v[3]; double dtfm,mbody[3],tbody[3],fquat[4]; double dtf2 = dtf * 2.0; - // compute target temperature - // update thermostat chains coupled to particles - - if (tstat_flag) { - compute_temp_target(); - nhc_temp_integrate(); - } - - // compute target pressure - // update epsilon dot - // update thermostat coupled to barostat - - if (pstat_flag) { - nhc_press_integrate(); - - if (pstyle == ISO) { - temperature->compute_scalar(); - pressure->compute_scalar(); - } else { - temperature->compute_vector(); - pressure->compute_vector(); - } - couple(); - pressure->addstep(update->ntimestep+1); - - compute_press_target(); - nh_epsilon_dot(); - } - // compute scale variables scale_t[0] = scale_t[1] = scale_t[2] = 1.0; scale_v[0] = scale_v[1] = scale_v[2] = 1.0; scale_r = 1.0; if (tstat_flag) { tmp = exp(-dtq * eta_dot_t[0]); scale_t[0] = scale_t[1] = scale_t[2] = tmp; tmp = exp(-dtq * eta_dot_r[0]); scale_r = tmp; } if (pstat_flag) { scale_t[0] *= exp(-dtq * (epsilon_dot[0] + mtk_term2)); scale_t[1] *= exp(-dtq * (epsilon_dot[1] + mtk_term2)); scale_t[2] *= exp(-dtq * (epsilon_dot[2] + mtk_term2)); scale_r *= exp(-dtq * (pdim * mtk_term2)); tmp = dtq * epsilon_dot[0]; scale_v[0] = dtv * exp(tmp) * maclaurin_series(tmp); tmp = dtq * epsilon_dot[1]; scale_v[1] = dtv * exp(tmp) * maclaurin_series(tmp); tmp = dtq * epsilon_dot[2]; scale_v[2] = dtv * exp(tmp) * maclaurin_series(tmp); } // update xcm, vcm, quat, conjqm and angmom for (int ibody = 0; ibody < nlocal_body; ibody++) { Body *b = &body[ibody]; // step 1.1 - update vcm by 1/2 step dtfm = dtf / b->mass; b->vcm[0] += dtfm * b->fcm[0]; b->vcm[1] += dtfm * b->fcm[1]; b->vcm[2] += dtfm * b->fcm[2]; if (tstat_flag || pstat_flag) { b->vcm[0] *= scale_t[0]; b->vcm[1] *= scale_t[1]; b->vcm[2] *= scale_t[2]; } // step 1.2 - update xcm by full step if (!pstat_flag) { b->xcm[0] += dtv * b->vcm[0]; b->xcm[1] += dtv * b->vcm[1]; b->xcm[2] += dtv * b->vcm[2]; } else { b->xcm[0] += scale_v[0] * b->vcm[0]; b->xcm[1] += scale_v[1] * b->vcm[1]; b->xcm[2] += scale_v[2] * b->vcm[2]; } // step 1.3 - apply torque (body coords) to quaternion momentum MathExtra::transpose_matvec(b->ex_space,b->ey_space,b->ez_space, b->torque,tbody); MathExtra::quatvec(b->quat,tbody,fquat); b->conjqm[0] += dtf2 * fquat[0]; b->conjqm[1] += dtf2 * fquat[1]; b->conjqm[2] += dtf2 * fquat[2]; b->conjqm[3] += dtf2 * fquat[3]; if (tstat_flag || pstat_flag) { b->conjqm[0] *= scale_r; b->conjqm[1] *= scale_r; b->conjqm[2] *= scale_r; b->conjqm[3] *= scale_r; } // step 1.4 to 1.13 - use no_squish rotate to update p and q no_squish_rotate(3,b->conjqm,b->quat,b->inertia,dtq); no_squish_rotate(2,b->conjqm,b->quat,b->inertia,dtq); no_squish_rotate(1,b->conjqm,b->quat,b->inertia,dtv); no_squish_rotate(2,b->conjqm,b->quat,b->inertia,dtq); no_squish_rotate(3,b->conjqm,b->quat,b->inertia,dtq); // update exyz_space // transform p back to angmom // update angular velocity MathExtra::q_to_exyz(b->quat,b->ex_space,b->ey_space, b->ez_space); MathExtra::invquatvec(b->quat,b->conjqm,mbody); MathExtra::matvec(b->ex_space,b->ey_space,b->ez_space, mbody,b->angmom); b->angmom[0] *= 0.5; b->angmom[1] *= 0.5; b->angmom[2] *= 0.5; MathExtra::angmom_to_omega(b->angmom,b->ex_space,b->ey_space, b->ez_space,b->inertia,b->omega); } - // virial setup before call to set_xv - - if (vflag) v_setup(vflag); - else evflag = 0; - // forward communicate updated info of all bodies commflag = INITIAL; comm->forward_comm_fix(this,26); // accumulate translational and rotational kinetic energies if (tstat_flag || pstat_flag) { akin_t = akin_r = 0.0; for (int ibody = 0; ibody < nlocal_body; ibody++) { Body *b = &body[ibody]; akin_t += b->mass*(b->vcm[0]*b->vcm[0] + b->vcm[1]*b->vcm[1] + b->vcm[2]*b->vcm[2]); akin_r += b->angmom[0]*b->omega[0] + b->angmom[1]*b->omega[1] + b->angmom[2]*b->omega[2]; } double ke[2],keall[2]; ke[0] = akin_t; ke[1] = akin_r; MPI_Allreduce(ke,keall,2,MPI_DOUBLE,MPI_SUM,world); akin_t = keall[0]; akin_r = keall[1]; } + // compute target temperature + // update thermostat chains using akin_t and akin_r + // refer to update_nhcp() in Kamberaj et al. + + if (tstat_flag) { + compute_temp_target(); + nhc_temp_integrate(); + } + + // update thermostat chains coupled with barostat + // refer to update_nhcb() in Kamberaj et al. + + if (pstat_flag) { + nhc_press_integrate(); + } + + // 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 FixRigidNHSmall::final_integrate() { int i,ibody; double tmp,scale_t[3],scale_r; double dtfm; double mbody[3],tbody[3],fquat[4]; double dtf2 = dtf * 2.0; // compute scale variables scale_t[0] = scale_t[1] = scale_t[2] = 1.0; scale_r = 1.0; if (tstat_flag) { tmp = exp(-1.0 * dtq * eta_dot_t[0]); scale_t[0] = scale_t[1] = scale_t[2] = tmp; scale_r = exp(-1.0 * dtq * eta_dot_r[0]); } if (pstat_flag) { scale_t[0] *= exp(-dtq * (epsilon_dot[0] + mtk_term2)); scale_t[1] *= exp(-dtq * (epsilon_dot[1] + mtk_term2)); scale_t[2] *= exp(-dtq * (epsilon_dot[2] + mtk_term2)); scale_r *= exp(-dtq * (pdim * mtk_term2)); } // sum over atoms to get force and torque on rigid body imageint *image = atom->image; double **x = atom->x; double **f = atom->f; int nlocal = atom->nlocal; double dx,dy,dz; double unwrap[3]; double *xcm,*fcm,*tcm; for (ibody = 0; ibody < nlocal_body+nghost_body; ibody++) { fcm = body[ibody].fcm; fcm[0] = fcm[1] = fcm[2] = 0.0; tcm = body[ibody].torque; tcm[0] = tcm[1] = tcm[2] = 0.0; } for (i = 0; i < nlocal; i++) { if (atom2body[i] < 0) continue; Body *b = &body[atom2body[i]]; fcm = b->fcm; fcm[0] += f[i][0]; fcm[1] += f[i][1]; fcm[2] += f[i][2]; domain->unmap(x[i],image[i],unwrap); xcm = b->xcm; dx = unwrap[0] - xcm[0]; dy = unwrap[1] - xcm[1]; dz = unwrap[2] - xcm[2]; tcm = b->torque; tcm[0] += dy*f[i][2] - dz*f[i][1]; tcm[1] += dz*f[i][0] - dx*f[i][2]; tcm[2] += dx*f[i][1] - dy*f[i][0]; } // extended particles add their torque to torque of body if (extended) { double **torque = atom->torque; for (i = 0; i < nlocal; i++) { if (atom2body[i] < 0) continue; if (eflags[i] & TORQUE) { tcm = body[atom2body[i]].torque; tcm[0] += torque[i][0]; tcm[1] += torque[i][1]; tcm[2] += torque[i][2]; } } } // reverse communicate fcm, torque of all bodies commflag = FORCE_TORQUE; comm->reverse_comm_fix(this,6); // include Langevin thermostat forces and torques if (langflag) { for (int ibody = 0; ibody < nlocal_body; ibody++) { fcm = body[ibody].fcm; fcm[0] += langextra[ibody][0]; fcm[1] += langextra[ibody][1]; fcm[2] += langextra[ibody][2]; tcm = body[ibody].torque; tcm[0] += langextra[ibody][3]; tcm[1] += langextra[ibody][4]; tcm[2] += langextra[ibody][5]; } } // update vcm and angmom // include Langevin thermostat forces // fflag,tflag = 0 for some dimensions in 2d for (ibody = 0; ibody < nbody; ibody++) { Body *b = &body[ibody]; // update vcm by 1/2 step dtfm = dtf / b->mass; if (tstat_flag || pstat_flag) { b->vcm[0] *= scale_t[0]; b->vcm[1] *= scale_t[1]; b->vcm[2] *= scale_t[2]; } b->vcm[0] += dtfm * b->fcm[0]; b->vcm[1] += dtfm * b->fcm[1]; b->vcm[2] += dtfm * b->fcm[2]; // update conjqm, then transform to angmom, set velocity again // virial is already setup from initial_integrate MathExtra::transpose_matvec(b->ex_space,b->ey_space, b->ez_space,b->torque,tbody); MathExtra::quatvec(b->quat,tbody,fquat); if (tstat_flag || pstat_flag) { b->conjqm[0] = scale_r * b->conjqm[0] + dtf2 * fquat[0]; b->conjqm[1] = scale_r * b->conjqm[1] + dtf2 * fquat[1]; b->conjqm[2] = scale_r * b->conjqm[2] + dtf2 * fquat[2]; b->conjqm[3] = scale_r * b->conjqm[3] + dtf2 * fquat[3]; } else { b->conjqm[0] += dtf2 * fquat[0]; b->conjqm[1] += dtf2 * fquat[1]; b->conjqm[2] += dtf2 * fquat[2]; b->conjqm[3] += dtf2 * fquat[3]; } MathExtra::invquatvec(b->quat,b->conjqm,mbody); MathExtra::matvec(b->ex_space,b->ey_space,b->ez_space,mbody,b->angmom); b->angmom[0] *= 0.5; b->angmom[1] *= 0.5; b->angmom[2] *= 0.5; MathExtra::angmom_to_omega(b->angmom,b->ex_space,b->ey_space, b->ez_space,b->inertia,b->omega); } // forward communicate updated info of all bodies commflag = FINAL; comm->forward_comm_fix(this,10); // accumulate translational and rotational kinetic energies if (pstat_flag) { akin_t = akin_r = 0.0; for (int ibody = 0; ibody < nlocal_body; ibody++) { Body *b = &body[ibody]; akin_t += b->mass*(b->vcm[0]*b->vcm[0] + b->vcm[1]*b->vcm[1] + b->vcm[2]*b->vcm[2]); akin_r += b->angmom[0]*b->omega[0] + b->angmom[1]*b->omega[1] + b->angmom[2]*b->omega[2]; } double ke[2],keall[2]; ke[0] = akin_t; ke[1] = akin_r; MPI_Allreduce(ke,keall,2,MPI_DOUBLE,MPI_SUM,world); akin_t = keall[0]; akin_r = keall[1]; } // set velocity/rotation of atoms in rigid bodies // virial is already setup from initial_integrate set_v(); - // compute temperature and pressure tensor - // couple to compute current pressure components - // trigger virial computation on next timestep - + // compute current temperature if (tcomputeflag) t_current = temperature->compute_scalar(); + + // compute current and target pressures + // update epsilon dot using akin_t and akin_r + if (pstat_flag) { - if (pstyle == ISO) pressure->compute_scalar(); - else pressure->compute_vector(); + if (pstyle == ISO) { + temperature->compute_scalar(); + pressure->compute_scalar(); + } else { + temperature->compute_vector(); + 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(); + compute_press_target(); + + nh_epsilon_dot(); + } } /* ---------------------------------------------------------------------- */ void FixRigidNHSmall::nhc_temp_integrate() { int i,j,k; double kt,gfkt_t,gfkt_r,tmp,ms,s,s2; kt = boltz * t_target; gfkt_t = nf_t * kt; gfkt_r = nf_r * kt; // update thermostat masses double t_mass = boltz * t_target / (t_freq * t_freq); q_t[0] = nf_t * t_mass; q_r[0] = nf_r * t_mass; for (i = 1; i < t_chain; i++) q_t[i] = q_r[i] = t_mass; // update force of thermostats coupled to particles f_eta_t[0] = (akin_t * mvv2e - gfkt_t) / q_t[0]; f_eta_r[0] = (akin_r * mvv2e - gfkt_r) / q_r[0]; // multiple timestep iteration for (i = 0; i < t_iter; i++) { for (j = 0; j < t_order; j++) { // update thermostat velocities half step eta_dot_t[t_chain-1] += wdti2[j] * f_eta_t[t_chain-1]; eta_dot_r[t_chain-1] += wdti2[j] * f_eta_r[t_chain-1]; for (k = 1; k < t_chain; k++) { tmp = wdti4[j] * eta_dot_t[t_chain-k]; ms = maclaurin_series(tmp); s = exp(-1.0 * tmp); s2 = s * s; eta_dot_t[t_chain-k-1] = eta_dot_t[t_chain-k-1] * s2 + wdti2[j] * f_eta_t[t_chain-k-1] * s * ms; tmp = wdti4[j] * eta_dot_r[t_chain-k]; ms = maclaurin_series(tmp); s = exp(-1.0 * tmp); s2 = s * s; eta_dot_r[t_chain-k-1] = eta_dot_r[t_chain-k-1] * s2 + wdti2[j] * f_eta_r[t_chain-k-1] * s * ms; } // update thermostat positions a full step for (k = 0; k < t_chain; k++) { eta_t[k] += wdti1[j] * eta_dot_t[k]; eta_r[k] += wdti1[j] * eta_dot_r[k]; } // update thermostat forces for (k = 1; k < t_chain; k++) { f_eta_t[k] = q_t[k-1] * eta_dot_t[k-1] * eta_dot_t[k-1] - kt; f_eta_t[k] /= q_t[k]; f_eta_r[k] = q_r[k-1] * eta_dot_r[k-1] * eta_dot_r[k-1] - kt; f_eta_r[k] /= q_r[k]; } // update thermostat velocities a full step for (k = 0; k < t_chain-1; k++) { tmp = wdti4[j] * eta_dot_t[k+1]; ms = maclaurin_series(tmp); s = exp(-1.0 * tmp); s2 = s * s; eta_dot_t[k] = eta_dot_t[k] * s2 + wdti2[j] * f_eta_t[k] * s * ms; tmp = q_t[k] * eta_dot_t[k] * eta_dot_t[k] - kt; f_eta_t[k+1] = tmp / q_t[k+1]; tmp = wdti4[j] * eta_dot_r[k+1]; ms = maclaurin_series(tmp); s = exp(-1.0 * tmp); s2 = s * s; eta_dot_r[k] = eta_dot_r[k] * s2 + wdti2[j] * f_eta_r[k] * s * ms; tmp = q_r[k] * eta_dot_r[k] * eta_dot_r[k] - kt; - f_eta_r[k+1] = tmp / q_r[k+1]; + f_eta_r[k+1] = tmp / q_r[k+1]; } eta_dot_t[t_chain-1] += wdti2[j] * f_eta_t[t_chain-1]; eta_dot_r[t_chain-1] += wdti2[j] * f_eta_r[t_chain-1]; } } } /* ---------------------------------------------------------------------- */ void FixRigidNHSmall::nhc_press_integrate() { - int i,k; + int i,j,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; + 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-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]; } + kecurrent /= pdim; f_eta_b[0] = (kecurrent - lkt_press) / q_b[0]; + + // multiple timestep iteration - // update thermostat velocities a half step + for (i = 0; i < t_iter; i++) { + for (j = 0; j < t_order; j++) { + + // update thermostat velocities a half step - eta_dot_b[p_chain-1] += 0.5 * dtq * f_eta_b[p_chain-1]; + eta_dot_b[p_chain-1] += wdti2[j] * 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; - } + for (k = 1; k < p_chain; k++) { + tmp = wdti4[j] * eta_dot_b[p_chain-k]; + ms = maclaurin_series(tmp); + s = exp(-0.5 * tmp); + s2 = s * s; + eta_dot_b[p_chain-k-1] = eta_dot_b[p_chain-k-1] * s2 + + wdti2[j] * f_eta_b[p_chain-k-1] * s * ms; + } - // update thermostat positions + // update thermostat positions - for (k = 0; k < p_chain; k++) - eta_b[k] += dtv * eta_dot_b[k]; + for (k = 0; k < p_chain; k++) + eta_b[k] += wdti1[j] * eta_dot_b[k]; - // update epsilon dot + // update thermostat forces - 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]; + for (k = 1; k < p_chain; k++) { + f_eta_b[k] = q_b[k-1] * eta_dot_b[k-1] * eta_dot_b[k-1] - kt; + f_eta_b[k] /= q_b[k]; + } + + // update thermostat velocites a full step - // update thermostat velocites a full step + for (k = 0; k < p_chain-1; k++) { + tmp = wdti4[j] * 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 + wdti2[j] * 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]; + } - 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] += wdti2[j] * f_eta_b[p_chain-1]; + } } - - eta_dot_b[p_chain-1] += 0.5 * dtq * f_eta_b[p_chain-1]; - } /* ---------------------------------------------------------------------- compute kinetic energy in the extended Hamiltonian conserved quantity = sum of returned energy and potential energy -----------------------------------------------------------------------*/ double FixRigidNHSmall::compute_scalar() { int i,k; double kt = boltz * t_target; double energy,ke_t,ke_q,tmp,Pkq[4]; double *vcm,*quat; // compute the kinetic parts of H_NVE in Kameraj et al (JCP 2005, pp 224114) // translational and rotational kinetic energies ke_t = 0.0; ke_q = 0.0; for (int i = 0; i < nlocal_body; i++) { vcm = body[i].vcm; quat = body[i].quat; ke_t += body[i].mass * (vcm[0]*vcm[0] + vcm[1]*vcm[1] + vcm[2]*vcm[2]); for (k = 1; k < 4; k++) { if (k == 1) { Pkq[0] = -quat[1]; Pkq[1] = quat[0]; Pkq[2] = quat[3]; Pkq[3] = -quat[2]; } else if (k == 2) { Pkq[0] = -quat[2]; Pkq[1] = -quat[3]; Pkq[2] = quat[0]; Pkq[3] = quat[1]; } else if (k == 3) { Pkq[0] = -quat[3]; Pkq[1] = quat[2]; Pkq[2] = -quat[1]; Pkq[3] = quat[0]; } tmp = body[i].conjqm[0]*Pkq[0] + body[i].conjqm[1]*Pkq[1] + body[i].conjqm[2]*Pkq[2] + body[i].conjqm[3]*Pkq[3]; tmp *= tmp; if (fabs(body[i].inertia[k-1]) < 1e-6) tmp = 0.0; else tmp /= (8.0 * body[i].inertia[k-1]); ke_q += tmp; } } double ke[2],keall[2]; ke[0] = ke_t; ke[1] = ke_q; MPI_Allreduce(ke,keall,2,MPI_DOUBLE,MPI_SUM,world); ke_t = keall[0]; ke_q = keall[1]; energy = (ke_t + ke_q) * mvv2e; if (tstat_flag) { // thermostat chain energy: from equation 12 in Kameraj et al (JCP 2005) energy += kt * (nf_t * eta_t[0] + nf_r * eta_r[0]); for (i = 1; i < t_chain; i++) energy += kt * (eta_t[i] + eta_r[i]); for (i = 0; i < t_chain; i++) { energy += 0.5 * q_t[i] * (eta_dot_t[i] * eta_dot_t[i]); energy += 0.5 * q_r[i] * (eta_dot_r[i] * eta_dot_r[i]); } } if (pstat_flag) { // using equation 22 in Kameraj et al for H_NPT + double e = 0.0; for (i = 0; i < 3; i++) - energy += 0.5 * epsilon_mass[i] * epsilon_dot[i] * epsilon_dot[i]; + if (p_flag[i]) + e += epsilon_mass[i] * epsilon_dot[i] * epsilon_dot[i]; + energy += e*(0.5/pdim); double vol; if (dimension == 2) vol = domain->xprd * domain->yprd; else vol = domain->xprd * domain->yprd * domain->zprd; double p0 = (p_target[0] + p_target[1] + p_target[2]) / 3.0; energy += p0 * vol / nktv2p; - + for (i = 0; i < p_chain; i++) { energy += kt * eta_b[i]; energy += 0.5 * q_b[i] * (eta_dot_b[i] * eta_dot_b[i]); } } return energy; } /* ---------------------------------------------------------------------- */ void FixRigidNHSmall::couple() { double *tensor = pressure->vector; if (pstyle == ISO) { p_current[0] = p_current[1] = p_current[2] = pressure->scalar; } else if (pcouple == XYZ) { double ave = 1.0/3.0 * (tensor[0] + tensor[1] + tensor[2]); p_current[0] = p_current[1] = p_current[2] = ave; } else if (pcouple == XY) { double ave = 0.5 * (tensor[0] + tensor[1]); p_current[0] = p_current[1] = ave; p_current[2] = tensor[2]; } else if (pcouple == YZ) { double ave = 0.5 * (tensor[1] + tensor[2]); p_current[1] = p_current[2] = ave; p_current[0] = tensor[0]; } else if (pcouple == XZ) { double ave = 0.5 * (tensor[0] + tensor[2]); p_current[0] = p_current[2] = ave; p_current[1] = tensor[1]; } else { p_current[0] = tensor[0]; p_current[1] = tensor[1]; p_current[2] = tensor[2]; } } /* ---------------------------------------------------------------------- */ void FixRigidNHSmall::remap() { int i; double oldlo,oldhi,ctr,expfac; double **x = atom->x; int *mask = atom->mask; int nlocal = atom->nlocal; // epsilon is not used, except for book-keeping for (i = 0; i < 3; i++) epsilon[i] += dtq * epsilon_dot[i]; // convert pertinent atoms and rigid bodies to lamda coords if (allremap) domain->x2lamda(nlocal); else { for (i = 0; i < nlocal; i++) if (mask[i] & dilate_group_bit) domain->x2lamda(x[i],x[i]); } if (nrigidfix) for (i = 0; i < nrigidfix; i++) modify->fix[rfix[i]]->deform(0); // reset global and local box to new size/shape for (i = 0; i < 3; i++) { if (p_flag[i]) { oldlo = domain->boxlo[i]; oldhi = domain->boxhi[i]; ctr = 0.5 * (oldlo + oldhi); expfac = exp(dtq * epsilon_dot[i]); domain->boxlo[i] = (oldlo-ctr)*expfac + ctr; domain->boxhi[i] = (oldhi-ctr)*expfac + ctr; } } domain->set_global_box(); domain->set_local_box(); // convert pertinent atoms and rigid bodies back to box coords if (allremap) domain->lamda2x(nlocal); else { for (i = 0; i < nlocal; i++) if (mask[i] & dilate_group_bit) domain->lamda2x(x[i],x[i]); } if (nrigidfix) for (i = 0; i< nrigidfix; i++) modify->fix[rfix[i]]->deform(1); } /* ---------------------------------------------------------------------- compute target temperature and kinetic energy -----------------------------------------------------------------------*/ void FixRigidNHSmall::compute_temp_target() { double delta = update->ntimestep - update->beginstep; if (delta != 0.0) delta /= update->endstep - update->beginstep; t_target = t_start + delta * (t_stop-t_start); } /* ---------------------------------------------------------------------- compute hydrostatic target pressure -----------------------------------------------------------------------*/ void FixRigidNHSmall::compute_press_target() { double delta = update->ntimestep - update->beginstep; if (delta != 0.0) delta /= update->endstep - update->beginstep; p_hydro = 0.0; for (int i = 0; i < 3; i++) if (p_flag[i]) { p_target[i] = p_start[i] + delta * (p_stop[i]-p_start[i]); p_hydro += p_target[i]; } p_hydro /= pdim; } -/* ---------------------------------------------------------------------- - apply evolution operators to quat, quat momentum - see Miller paper cited in fix rigid/nvt and fix rigid/npt -------------------------------------------------------------------------- */ - -void FixRigidNHSmall::no_squish_rotate(int k, double *p, double *q, - double *inertia, double dt) -{ - double phi,c_phi,s_phi,kp[4],kq[4]; - - // apply permuation operator on p and q, get kp and kq - - if (k == 1) { - kq[0] = -q[1]; kp[0] = -p[1]; - kq[1] = q[0]; kp[1] = p[0]; - kq[2] = q[3]; kp[2] = p[3]; - kq[3] = -q[2]; kp[3] = -p[2]; - } else if (k == 2) { - kq[0] = -q[2]; kp[0] = -p[2]; - kq[1] = -q[3]; kp[1] = -p[3]; - kq[2] = q[0]; kp[2] = p[0]; - kq[3] = q[1]; kp[3] = p[1]; - } else if (k == 3) { - kq[0] = -q[3]; kp[0] = -p[3]; - kq[1] = q[2]; kp[1] = p[2]; - kq[2] = -q[1]; kp[2] = -p[1]; - kq[3] = q[0]; kp[3] = p[0]; - } - - // obtain phi, cosines and sines - - phi = p[0]*kq[0] + p[1]*kq[1] + p[2]*kq[2] + p[3]*kq[3]; - if (fabs(inertia[k-1]) < 1e-6) phi *= 0.0; - else phi /= 4.0 * inertia[k-1]; - c_phi = cos(dt * phi); - s_phi = sin(dt * phi); - - // advance p and q - - p[0] = c_phi*p[0] + s_phi*kp[0]; - p[1] = c_phi*p[1] + s_phi*kp[1]; - p[2] = c_phi*p[2] + s_phi*kp[2]; - p[3] = c_phi*p[3] + s_phi*kp[3]; - - q[0] = c_phi*q[0] + s_phi*kq[0]; - q[1] = c_phi*q[1] + s_phi*kq[1]; - q[2] = c_phi*q[2] + s_phi*kq[2]; - q[3] = c_phi*q[3] + s_phi*kq[3]; -} - /* ---------------------------------------------------------------------- update epsilon_dot -----------------------------------------------------------------------*/ void FixRigidNHSmall::nh_epsilon_dot() { int i; double volume,scale,f_epsilon; if (dimension == 2) volume = domain->xprd*domain->yprd; else volume = domain->xprd*domain->yprd*domain->zprd; // MTK terms mtk_term1 = (akin_t + akin_r) * mvv2e / g_f; scale = exp(-1.0 * dtq * eta_dot_b[0]); for (i = 0; i < 3; i++) if (p_flag[i]) { f_epsilon = (p_current[i]-p_hydro)*volume / nktv2p + mtk_term1; f_epsilon /= epsilon_mass[i]; epsilon_dot[i] += dtq * f_epsilon; epsilon_dot[i] *= scale; } mtk_term2 = 0.0; for (i = 0; i < 3; i++) if (p_flag[i]) mtk_term2 += epsilon_dot[i]; mtk_term2 /= g_f; } /* ---------------------------------------------------------------------- pack entire state of Fix into one write ------------------------------------------------------------------------- */ void FixRigidNHSmall::write_restart(FILE *fp) { if (tstat_flag == 0 && pstat_flag == 0) return; int nsize = 2; // tstat_flag and pstat_flag if (tstat_flag) { nsize += 1; // t_chain nsize += 4*t_chain; // eta_t, eta_r, eta_dot_t, eta_dot_r } if (pstat_flag) { nsize += 7; // p_chain, epsilon(3) and epsilon_dot(3) nsize += 2*p_chain; } double *list; memory->create(list,nsize,"rigid_nh:list"); int n = 0; list[n++] = tstat_flag; if (tstat_flag) { list[n++] = t_chain; for (int i = 0; i < t_chain; i++) { list[n++] = eta_t[i]; list[n++] = eta_r[i]; list[n++] = eta_dot_t[i]; list[n++] = eta_dot_r[i]; } } list[n++] = pstat_flag; if (pstat_flag) { list[n++] = epsilon[0]; list[n++] = epsilon[1]; list[n++] = epsilon[2]; list[n++] = epsilon_dot[0]; list[n++] = epsilon_dot[1]; list[n++] = epsilon_dot[2]; list[n++] = p_chain; for (int i = 0; i < p_chain; i++) { list[n++] = eta_b[i]; list[n++] = eta_dot_b[i]; } } if (comm->me == 0) { int size = (nsize)*sizeof(double); fwrite(&size,sizeof(int),1,fp); fwrite(list,sizeof(double),nsize,fp); } memory->destroy(list); } /* ---------------------------------------------------------------------- use state info from restart file to restart the Fix ------------------------------------------------------------------------- */ void FixRigidNHSmall::restart(char *buf) { int n = 0; double *list = (double *) buf; int flag = static_cast (list[n++]); if (flag) { int m = static_cast (list[n++]); if (tstat_flag && m == t_chain) { for (int i = 0; i < t_chain; i++) { eta_t[i] = list[n++]; eta_r[i] = list[n++]; eta_dot_t[i] = list[n++]; eta_dot_r[i] = list[n++]; } } else n += 4*m; } flag = static_cast (list[n++]); if (flag) { epsilon[0] = list[n++]; epsilon[1] = list[n++]; epsilon[2] = list[n++]; epsilon_dot[0] = list[n++]; epsilon_dot[1] = list[n++]; epsilon_dot[2] = list[n++]; int m = static_cast (list[n++]); if (pstat_flag && m == p_chain) { for (int i = 0; i < p_chain; i++) { eta_b[i] = list[n++]; eta_dot_b[i] = list[n++]; } } else n += 2*m; } } /* ---------------------------------------------------------------------- */ int FixRigidNHSmall::modify_param(int narg, char **arg) { if (strcmp(arg[0],"temp") == 0) { if (narg < 2) error->all(FLERR,"Illegal fix_modify command"); if (!pstat_flag) error->all(FLERR,"Illegal fix_modify command"); if (tcomputeflag) { modify->delete_compute(id_temp); tcomputeflag = 0; } delete [] id_temp; int n = strlen(arg[1]) + 1; id_temp = new char[n]; strcpy(id_temp,arg[1]); int icompute = modify->find_compute(arg[1]); if (icompute < 0) error->all(FLERR,"Could not find fix_modify temperature ID"); temperature = modify->compute[icompute]; if (temperature->tempflag == 0) error->all(FLERR, "Fix_modify temperature ID does not compute temperature"); if (temperature->igroup != 0 && comm->me == 0) error->warning(FLERR,"Temperature for fix modify is not for group all"); // reset id_temp of pressure to new temperature ID if (pstat_flag) { icompute = modify->find_compute(id_press); if (icompute < 0) error->all(FLERR,"Pressure ID for fix modify does not exist"); modify->compute[icompute]->reset_extra_compute_fix(id_temp); } return 2; } else if (strcmp(arg[0],"press") == 0) { if (narg < 2) error->all(FLERR,"Illegal fix_modify command"); if (!pstat_flag) error->all(FLERR,"Illegal fix_modify command"); if (pcomputeflag) { modify->delete_compute(id_press); pcomputeflag = 0; } delete [] id_press; int n = strlen(arg[1]) + 1; id_press = new char[n]; strcpy(id_press,arg[1]); int icompute = modify->find_compute(arg[1]); if (icompute < 0) error->all(FLERR,"Could not find fix_modify pressure ID"); pressure = modify->compute[icompute]; if (pressure->pressflag == 0) error->all(FLERR,"Fix_modify pressure ID does not compute pressure"); return 2; } return 0; } /* ---------------------------------------------------------------------- */ void FixRigidNHSmall::allocate_chain() { if (tstat_flag) { q_t = new double[t_chain]; q_r = new double[t_chain]; eta_t = new double[t_chain]; eta_r = new double[t_chain]; eta_dot_t = new double[t_chain]; eta_dot_r = new double[t_chain]; f_eta_t = new double[t_chain]; f_eta_r = new double[t_chain]; } if (pstat_flag) { q_b = new double[p_chain]; eta_b = new double[p_chain]; eta_dot_b = new double[p_chain]; f_eta_b = new double[p_chain]; } } /* ---------------------------------------------------------------------- */ void FixRigidNHSmall::reset_target(double t_new) { t_start = t_stop = t_new; } /* ---------------------------------------------------------------------- */ void FixRigidNHSmall::allocate_order() { w = new double[t_order]; wdti1 = new double[t_order]; wdti2 = new double[t_order]; wdti4 = new double[t_order]; } /* ---------------------------------------------------------------------- */ void FixRigidNHSmall::deallocate_chain() { if (tstat_flag) { delete [] q_t; delete [] q_r; delete [] eta_t; delete [] eta_r; delete [] eta_dot_t; delete [] eta_dot_r; delete [] f_eta_t; delete [] f_eta_r; } if (pstat_flag) { delete [] q_b; delete [] eta_b; delete [] eta_dot_b; delete [] f_eta_b; } } /* ---------------------------------------------------------------------- */ void FixRigidNHSmall::deallocate_order() { delete [] w; delete [] wdti1; delete [] wdti2; delete [] wdti4; } diff --git a/src/RIGID/fix_rigid_nh_small.h b/src/RIGID/fix_rigid_nh_small.h index 68b47bbf8..cc148358c 100644 --- a/src/RIGID/fix_rigid_nh_small.h +++ b/src/RIGID/fix_rigid_nh_small.h @@ -1,195 +1,194 @@ /* -*- c++ -*- ---------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ #ifndef LMP_FIX_RIGID_NH_SMALL_H #define LMP_FIX_RIGID_NH_SMALL_H #include "fix_rigid_small.h" namespace LAMMPS_NS { class FixRigidNHSmall : public FixRigidSmall { public: FixRigidNHSmall(class LAMMPS *, int, char **); virtual ~FixRigidNHSmall(); virtual int setmask(); virtual void init(); virtual void setup(int); virtual void initial_integrate(int); virtual void final_integrate(); virtual double compute_scalar(); int modify_param(int, char **); void write_restart(FILE *); void restart(char *buf); void reset_target(double); protected: double boltz,nktv2p,mvv2e; // boltzman constant, conversion factors int dimension; // # of dimensions int nf_t,nf_r; // trans/rot degrees of freedom double onednft,onednfr; // factors 1 + dimension/trans(rot) degrees of freedom double *w,*wdti1,*wdti2,*wdti4; // Yoshida-Suzuki coefficients double *q_t,*q_r; // trans/rot thermostat masses double *eta_t,*eta_r; // trans/rot thermostat positions double *eta_dot_t,*eta_dot_r; // trans/rot thermostat velocities double *f_eta_t,*f_eta_r; // trans/rot thermostat forces double epsilon_mass[3], *q_b; // baro/thermo masses double epsilon[3],*eta_b; // baro/thermo positions double epsilon_dot[3],*eta_dot_b; // baro/thermo velocities double *f_eta_b; // thermo forces double akin_t,akin_r; // translational/rotational kinetic energies int kspace_flag; // 1 if KSpace invoked, 0 if not int nrigidfix; // number of rigid fixes int *rfix; // indicies of rigid fixes double vol0; // reference volume double t0; // reference temperature int pdim,g_f; // number of barostatted dims, total DoFs double p_hydro; // hydrostatic target pressure double p_freq_max; // maximum barostat frequency double mtk_term1,mtk_term2; // Martyna-Tobias-Klein corrections double t_target,t_current; double t_freq; char *id_temp,*id_press; class Compute *temperature,*pressure; int tcomputeflag,pcomputeflag; void couple(); void remap(); void nhc_temp_integrate(); void nhc_press_integrate(); virtual void compute_temp_target(); void compute_press_target(); void nh_epsilon_dot(); void allocate_chain(); void allocate_order(); void deallocate_chain(); void deallocate_order(); - void no_squish_rotate(int, double *, double *, double *, double); inline double maclaurin_series(double); }; inline double FixRigidNHSmall::maclaurin_series(double x) { double x2,x4; x2 = x * x; x4 = x2 * x2; return (1.0 + (1.0/6.0) * x2 + (1.0/120.0) * x4 + (1.0/5040.0) * x2 * x4 + (1.0/362880.0) * x4 * x4); } } #endif /* ERROR/WARNING messages: E: Fix rigid npt/nph period must be > 0.0 Self-explanatory. E: Invalid fix rigid npt/nph command for a 2d simulation Cannot control z dimension in a 2d model. E: Invalid fix rigid npt/nph command pressure settings If multiple dimensions are coupled, those dimensions must be specified. E: Cannot use fix rigid npt/nph on a non-periodic dimension When specifying a diagonal pressure component, the dimension must be periodic. E: Invalid fix rigid npt/nph pressure settings Settings for coupled dimensions must be the same. E: Fix rigid nvt/npt/nph damping parameters must be > 0.0 Self-explanatory. E: Fix rigid npt/nph dilate group ID does not exist Self-explanatory. E: Temp ID for fix rigid npt/nph does not exist Specified compute temperature must be valid. E: fix rigid npt/nph does not yet allow triclinic box Self-explanatory. E: Cannot use fix rigid npt/nph and fix deform on same component of stress tensor This would be changing the same box dimension twice. E: Press ID for fix rigid npt/nph does not exist Specified compute pressure must be valid. E: Illegal ... command Self-explanatory. Check the input script syntax and compare to the documentation for the command. You can use -echo screen as a command-line option when running LAMMPS to see the offending line. E: Could not find fix_modify temperature ID The compute ID for computing temperature does not exist. E: Fix_modify temperature ID does not compute temperature The compute ID assigned to the fix must compute temperature. W: Temperature for fix modify is not for group all The temperature compute is being used with a pressure calculation which does operate on group all, so this may be inconsistent. E: Pressure ID for fix modify does not exist Self-explanatory. E: Could not find fix_modify pressure ID The compute ID for computing pressure does not exist. E: Fix_modify pressure ID does not compute pressure The compute ID assigned to the fix must compute pressure. U: Target temperature for fix rigid nvt/npt cannot be 0.0 Self-explanatory. U: Temperature ID for fix rigid npt/nph does not exist Self-explanatory. U: Pressure ID for fix rigid npt/nph does not exist Self-explanatory. */