diff --git a/src/GRANULAR/pair_gran_hertz_history.cpp b/src/GRANULAR/pair_gran_hertz_history.cpp index 9960df2c8f..719bbae98c 100644 --- a/src/GRANULAR/pair_gran_hertz_history.cpp +++ b/src/GRANULAR/pair_gran_hertz_history.cpp @@ -23,7 +23,7 @@ #include "atom.h" #include "update.h" #include "force.h" -#include "fix_rigid.h" +#include "fix.h" #include "neighbor.h" #include "neigh_list.h" #include "comm.h" @@ -60,10 +60,12 @@ void PairGranHertzHistory::compute(int eflag, int vflag) int shearupdate = 1; if (update->setupflag) shearupdate = 0; - // update body ptr and values for ghost atoms if using FixRigid masses + // update rigid body ptrs and values for ghost atoms if using FixRigid masses if (fix_rigid && neighbor->ago == 0) { - body = fix_rigid->body; + int tmp; + body = (int *) fix_rigid->extract("body",tmp); + mass_rigid = (double *) fix_rigid->extract("masstotal",tmp); comm->forward_comm_pair(this); } @@ -162,8 +164,8 @@ void PairGranHertzHistory::compute(int eflag, int vflag) mj = mass[type[j]]; } if (fix_rigid) { - if (body[i] >= 0) mi = fix_rigid->masstotal[body[i]]; - if (body[j] >= 0) mj = fix_rigid->masstotal[body[j]]; + if (body[i] >= 0) mi = mass_rigid[body[i]]; + if (body[j] >= 0) mj = mass_rigid[body[j]]; } meff = mi*mj / (mi+mj); @@ -374,9 +376,11 @@ double PairGranHertzHistory::single(int i, int j, int itype, int jtype, if (fix_rigid) { // NOTE: need to make sure ghost atoms have updated body? // depends on where single() is called from - body = fix_rigid->body; - if (body[i] >= 0) mi = fix_rigid->masstotal[body[i]]; - if (body[j] >= 0) mj = fix_rigid->masstotal[body[j]]; + int tmp; + body = (int *) fix_rigid->extract("body",tmp); + mass_rigid = (double *) fix_rigid->extract("masstotal",tmp); + if (body[i] >= 0) mi = mass_rigid[body[i]]; + if (body[j] >= 0) mj = mass_rigid[body[j]]; } meff = mi*mj / (mi+mj); diff --git a/src/GRANULAR/pair_gran_hooke.cpp b/src/GRANULAR/pair_gran_hooke.cpp index 88075420b7..9272909870 100644 --- a/src/GRANULAR/pair_gran_hooke.cpp +++ b/src/GRANULAR/pair_gran_hooke.cpp @@ -21,7 +21,7 @@ #include "pair_gran_hooke.h" #include "atom.h" #include "force.h" -#include "fix_rigid.h" +#include "fix.h" #include "neighbor.h" #include "neigh_list.h" #include "comm.h" @@ -53,10 +53,12 @@ void PairGranHooke::compute(int eflag, int vflag) if (eflag || vflag) ev_setup(eflag,vflag); else evflag = vflag_fdotr = 0; - // update body ptr and values for ghost atoms if using FixRigid masses + // update rigid body ptrs and values for ghost atoms if using FixRigid masses if (fix_rigid && neighbor->ago == 0) { - body = fix_rigid->body; + int tmp; + body = (int *) fix_rigid->extract("body",tmp); + mass_rigid = (double *) fix_rigid->extract("masstotal",tmp); comm->forward_comm_pair(this); } @@ -142,8 +144,8 @@ void PairGranHooke::compute(int eflag, int vflag) mj = mass[type[j]]; } if (fix_rigid) { - if (body[i] >= 0) mi = fix_rigid->masstotal[body[i]]; - if (body[j] >= 0) mj = fix_rigid->masstotal[body[j]]; + if (body[i] >= 0) mi = mass_rigid[body[i]]; + if (body[j] >= 0) mj = mass_rigid[body[j]]; } meff = mi*mj / (mi+mj); @@ -291,9 +293,11 @@ double PairGranHooke::single(int i, int j, int itype, int jtype, double rsq, if (fix_rigid) { // NOTE: need to make sure ghost atoms have updated body? // depends on where single() is called from - body = fix_rigid->body; - if (body[i] >= 0) mi = fix_rigid->masstotal[body[i]]; - if (body[j] >= 0) mj = fix_rigid->masstotal[body[j]]; + int tmp; + body = (int *) fix_rigid->extract("body",tmp); + mass_rigid = (double *) fix_rigid->extract("masstotal",tmp); + if (body[i] >= 0) mi = mass_rigid[body[i]]; + if (body[j] >= 0) mj = mass_rigid[body[j]]; } meff = mi*mj / (mi+mj); diff --git a/src/GRANULAR/pair_gran_hooke_history.cpp b/src/GRANULAR/pair_gran_hooke_history.cpp index 876efc4cf4..a94ae78697 100644 --- a/src/GRANULAR/pair_gran_hooke_history.cpp +++ b/src/GRANULAR/pair_gran_hooke_history.cpp @@ -29,7 +29,6 @@ #include "fix.h" #include "fix_pour.h" #include "fix_shear_history.h" -#include "fix_rigid.h" #include "comm.h" #include "neighbor.h" #include "neigh_list.h" @@ -99,10 +98,12 @@ void PairGranHookeHistory::compute(int eflag, int vflag) int shearupdate = 1; if (update->setupflag) shearupdate = 0; - // update body ptr and values for ghost atoms if using FixRigid masses + // update rigid body ptrs and values for ghost atoms if using FixRigid masses if (fix_rigid && neighbor->ago == 0) { - body = fix_rigid->body; + int tmp; + body = (int *) fix_rigid->extract("body",tmp); + mass_rigid = (double *) fix_rigid->extract("masstotal",tmp); comm->forward_comm_pair(this); } @@ -201,8 +202,8 @@ void PairGranHookeHistory::compute(int eflag, int vflag) mj = mass[type[j]]; } if (fix_rigid) { - if (body[i] >= 0) mi = fix_rigid->masstotal[body[i]]; - if (body[j] >= 0) mj = fix_rigid->masstotal[body[j]]; + if (body[i] >= 0) mi = mass_rigid[body[i]]; + if (body[j] >= 0) mj = mass_rigid[body[j]]; } meff = mi*mj / (mi+mj); @@ -444,8 +445,8 @@ void PairGranHookeHistory::init_style() fix_rigid = NULL; for (i = 0; i < modify->nfix; i++) - if (strstr(modify->fix[i]->style,"rigid")) break; - if (i < modify->nfix) fix_rigid = (FixRigid *) modify->fix[i]; + if (modify->fix[i]->rigid_flag) break; + if (i < modify->nfix) fix_rigid = modify->fix[i]; // set maxrad_dynamic and maxrad_frozen for each type // include future Fix pour particles as dynamic @@ -654,9 +655,11 @@ double PairGranHookeHistory::single(int i, int j, int itype, int jtype, if (fix_rigid) { // NOTE: need to make sure ghost atoms have updated body? // depends on where single() is called from - body = fix_rigid->body; - if (body[i] >= 0) mi = fix_rigid->masstotal[body[i]]; - if (body[j] >= 0) mj = fix_rigid->masstotal[body[j]]; + int tmp; + body = (int *) fix_rigid->extract("body",tmp); + mass_rigid = (double *) fix_rigid->extract("masstotal",tmp); + if (body[i] >= 0) mi = mass_rigid[body[i]]; + if (body[j] >= 0) mj = mass_rigid[body[j]]; } meff = mi*mj / (mi+mj); diff --git a/src/GRANULAR/pair_gran_hooke_history.h b/src/GRANULAR/pair_gran_hooke_history.h index c0b8e9d6a2..4d1aca8b53 100644 --- a/src/GRANULAR/pair_gran_hooke_history.h +++ b/src/GRANULAR/pair_gran_hooke_history.h @@ -59,8 +59,9 @@ class PairGranHookeHistory : public Pair { double *maxrad_dynamic,*maxrad_frozen; class FixShearHistory *fix_history; - class FixRigid *fix_rigid; + class Fix *fix_rigid; int *body; + double *mass_rigid; void allocate(); }; diff --git a/src/Makefile b/src/Makefile index b0cf43b46b..5d3116f7ad 100755 --- a/src/Makefile +++ b/src/Makefile @@ -15,7 +15,7 @@ OBJ = $(SRC:.cpp=.o) PACKAGE = asphere class2 colloid dipole fld gpu granular kim \ kspace manybody mc meam molecule opt peri poems reax replica \ - shock srd xtc + rigid shock srd xtc PACKUSER = user-misc user-atc user-awpmd user-cg-cmm user-colvars \ user-cuda user-eff user-ewaldn user-omp user-molfile \ diff --git a/src/RIGID/Install.sh b/src/RIGID/Install.sh new file mode 100644 index 0000000000..3979831f8e --- /dev/null +++ b/src/RIGID/Install.sh @@ -0,0 +1,36 @@ +# Install/unInstall package files in LAMMPS + +if (test $1 = 1) then + + cp fix_rigid.cpp .. + cp fix_rigid_nh.cpp .. + cp fix_rigid_nph.cpp .. + cp fix_rigid_npt.cpp .. + cp fix_rigid_nve.cpp .. + cp fix_rigid_nvt.cpp .. + + cp fix_rigid.h .. + cp fix_rigid_nh.h .. + cp fix_rigid_nph.h .. + cp fix_rigid_npt.h .. + cp fix_rigid_nve.h .. + cp fix_rigid_nvt.h .. + +elif (test $1 = 0) then + + rm -f ../fix_rigid.cpp + rm -f ../fix_rigid_nh.cpp + rm -f ../fix_rigid_nph.cpp + rm -f ../fix_rigid_npt.cpp + rm -f ../fix_rigid_nve.cpp + rm -f ../fix_rigid_nvt.cpp + + rm -f ../fix_rigid.h + rm -f ../fix_rigid_nh.h + rm -f ../fix_rigid_nph.h + rm -f ../fix_rigid_npt.h + rm -f ../fix_rigid_nve.h + rm -f ../fix_rigid_nvt.h + +fi + diff --git a/src/RIGID/fix_rigid.cpp b/src/RIGID/fix_rigid.cpp new file mode 100644 index 0000000000..5514a48482 --- /dev/null +++ b/src/RIGID/fix_rigid.cpp @@ -0,0 +1,2386 @@ +/* ---------------------------------------------------------------------- + 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; + + // 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 # + // 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; + int *molecule = atom->molecule; + int nlocal = atom->nlocal; + + maxmol = -1; + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) maxmol = MAX(maxmol,molecule[i]); + + int itmp; + MPI_Allreduce(&maxmol,&itmp,1,MPI_INT,MPI_MAX,world); + maxmol = 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:ncount"); + 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; + + 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 = atoi(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 = atof(arg[iarg+1]); + t_stop = atof(arg[iarg+2]); + t_period = atof(arg[iarg+3]); + seed = atoi(arg[iarg+4]); + if (t_period <= 0.0) + error->all(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 = atof(arg[iarg+1]); + t_stop = atof(arg[iarg+2]); + t_period = atof(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] = atof(arg[iarg+1]); + p_stop[0] = p_stop[1] = p_stop[2] = atof(arg[iarg+2]); + p_period[0] = p_period[1] = p_period[2] = atof(arg[iarg+3]); + p_flag[0] = p_flag[1] = p_flag[2] = 1; + if (dimension == 2) { + p_start[2] = p_stop[2] = p_period[2] = 0.0; + p_flag[2] = 0; + } + iarg += 4; + + } else if (strcmp(arg[iarg],"aniso") == 0) { + if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid command"); + if (strcmp(style,"rigid/npt") != 0 && strcmp(style,"rigid/nph") != 0) + error->all(FLERR,"Illegal fix rigid command"); + p_start[0] = p_start[1] = p_start[2] = atof(arg[iarg+1]); + p_stop[0] = p_stop[1] = p_stop[2] = atof(arg[iarg+2]); + p_period[0] = p_period[1] = p_period[2] = atof(arg[iarg+3]); + p_flag[0] = p_flag[1] = p_flag[2] = 1; + if (dimension == 2) { + p_start[2] = p_stop[2] = p_period[2] = 0.0; + p_flag[2] = 0; + } + iarg += 4; + + } else if (strcmp(arg[iarg],"x") == 0) { + if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid command"); + p_start[0] = atof(arg[iarg+1]); + p_stop[0] = atof(arg[iarg+2]); + p_period[0] = atof(arg[iarg+3]); + p_flag[0] = 1; + iarg += 4; + + } else if (strcmp(arg[iarg],"y") == 0) { + if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid command"); + p_start[1] = atof(arg[iarg+1]); + p_stop[1] = atof(arg[iarg+2]); + p_period[1] = atof(arg[iarg+3]); + p_flag[1] = 1; + iarg += 4; + + } else if (strcmp(arg[iarg],"z") == 0) { + if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid command"); + p_start[2] = atof(arg[iarg+1]); + p_stop[2] = atof(arg[iarg+2]); + p_period[2] = atof(arg[iarg+3]); + p_flag[2] = 1; + iarg += 4; + + } else if (strcmp(arg[iarg],"couple") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix rigid command"); + if (strcmp(arg[iarg+1],"xyz") == 0) pcouple = XYZ; + else if (strcmp(arg[iarg+1],"xy") == 0) pcouple = XY; + else if (strcmp(arg[iarg+1],"yz") == 0) pcouple = YZ; + else if (strcmp(arg[iarg+1],"xz") == 0) pcouple = XZ; + else if (strcmp(arg[iarg+1],"none") == 0) pcouple = NONE; + else error->all(FLERR,"Illegal fix rigid command"); + iarg += 2; + + } else if (strcmp(arg[iarg],"dilate") == 0) { + if (iarg+2 > narg) + error->all(FLERR,"Illegal fix rigid nvt/npt/nph command"); + if (strcmp(arg[iarg+1],"all") == 0) allremap = 1; + else { + allremap = 0; + delete [] id_dilate; + int n = strlen(arg[iarg+1]) + 1; + id_dilate = new char[n]; + strcpy(id_dilate,arg[iarg+1]); + int idilate = group->find(id_dilate); + if (idilate == -1) + error->all(FLERR, + "Fix rigid nvt/npt/nph dilate group ID does not exist"); + } + iarg += 2; + + } else if (strcmp(arg[iarg],"tparam") == 0) { + if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid command"); + if (strcmp(style,"rigid/nvt") != 0 && strcmp(style,"rigid/npt") != 0) + error->all(FLERR,"Illegal fix rigid command"); + t_chain = atoi(arg[iarg+1]); + t_iter = atoi(arg[iarg+2]); + t_order = atoi(arg[iarg+3]); + iarg += 4; + + } else if (strcmp(arg[iarg],"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 = atoi(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]); + 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; + + // 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"); + + // 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); + } + + // firstflag = 1 triggers one-time initialization of rigid body attributes + + firstflag = 1; +} + +/* ---------------------------------------------------------------------- */ + +FixRigid::~FixRigid() +{ + // unregister callbacks to this fix from Atom class + + atom->delete_callback(id,0); + + delete random; + delete [] infile; + memory->destroy(mol2body); + + // 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; + + // 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; + + // one-time initialization of rigid body attributes + // extended flags, masstotal, COM, inertia tensor + + if (firstflag) setup_bodies(); + firstflag = 0; + + // 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; +} + +/* ---------------------------------------------------------------------- */ + +void FixRigid::setup(int vflag) +{ + int i,n,ibody; + double massone,radone; + + // vcm = velocity of center-of-mass of each rigid body + // fcm = force on center-of-mass of each rigid body + + double **v = atom->v; + double **f = atom->f; + double *rmass = atom->rmass; + double *mass = atom->mass; + int *type = atom->type; + 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]; + 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; + sum[ibody][3] += f[i][0]; + sum[ibody][4] += f[i][1]; + sum[ibody][5] += f[i][2]; + } + + MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world); + + 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]; + fcm[ibody][0] = all[ibody][3]; + fcm[ibody][1] = all[ibody][4]; + fcm[ibody][2] = all[ibody][5]; + } + + // angmom = angular momentum of each rigid body + // torque = torque on each rigid body + + tagint *image = atom->image; + double **x = atom->x; + + 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 xunwrap,yunwrap,zunwrap,dx,dy,dz; + + 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] += dy * massone*v[i][2] - dz * massone*v[i][1]; + sum[ibody][1] += dz * massone*v[i][0] - dx * massone*v[i][2]; + sum[ibody][2] += dx * massone*v[i][1] - dy * massone*v[i][0]; + 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 rotation/torque to angmom/torque 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 **torque_one = atom->torque; + 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][0] += SINERTIA*rmass[i] * radone*radone * omega_one[i][0]; + sum[ibody][1] += SINERTIA*rmass[i] * radone*radone * omega_one[i][1]; + sum[ibody][2] += SINERTIA*rmass[i] * radone*radone * omega_one[i][2]; + } else if (eflags[i] & LINE) { + radone = lbonus[line[i]].length; + sum[ibody][2] += LINERTIA*rmass[i] * radone*radone * omega_one[i][2]; + } + } + if (eflags[i] & ANGMOM) { + sum[ibody][0] += angmom_one[i][0]; + sum[ibody][1] += angmom_one[i][1]; + sum[ibody][2] += angmom_one[i][2]; + } + 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); + + for (ibody = 0; ibody < nbody; ibody++) { + angmom[ibody][0] = all[ibody][0]; + angmom[ibody][1] = all[ibody][1]; + angmom[ibody][2] = all[ibody][2]; + torque[ibody][0] = all[ibody][3]; + torque[ibody][1] = all[ibody][4]; + torque[ibody][2] = all[ibody][5]; + } + + // zero langextra in case Langevin thermostat not used + // no point to calling post_force() here since langextra + // is only added to fcm/torque in final_integrate() + + for (ibody = 0; ibody < nbody; ibody++) + for (i = 0; i < 6; i++) langextra[ibody][i] = 0.0; + + // virial setup before call to set_v + + if (vflag) v_setup(vflag); + 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; + 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,xy,xz,yz; + + // sum over atoms to get force and torque on rigid body + + tagint *image = atom->image; + double **x = atom->x; + double **f = atom->f; + int nlocal = atom->nlocal; + + double xprd = domain->xprd; + double yprd = domain->yprd; + double zprd = domain->zprd; + if (triclinic) { + xy = domain->xy; + xz = domain->xz; + yz = domain->yz; + } + + int xbox,ybox,zbox; + double xunwrap,yunwrap,zunwrap,dx,dy,dz; + for (ibody = 0; ibody < nbody; ibody++) + for (i = 0; i < 6; i++) sum[ibody][i] = 0.0; + + for (i = 0; i < nlocal; i++) { + if (body[i] < 0) continue; + ibody = body[i]; + + sum[ibody][0] += f[i][0]; + sum[ibody][1] += f[i][1]; + sum[ibody][2] += f[i][2]; + + xbox = (image[i] & IMGMASK) - IMGMAX; + ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; + zbox = (image[i] >> IMG2BITS) - IMGMAX; + + if (triclinic == 0) { + xunwrap = x[i][0] + xbox*xprd; + yunwrap = x[i][1] + ybox*yprd; + zunwrap = x[i][2] + zbox*zprd; + } else { + xunwrap = x[i][0] + xbox*xprd + ybox*xy + zbox*xz; + yunwrap = x[i][1] + ybox*yprd + zbox*yz; + zunwrap = x[i][2] + zbox*zprd; + } + + dx = xunwrap - xcm[ibody][0]; + dy = yunwrap - xcm[ibody][1]; + dz = zunwrap - xcm[ibody][2]; + + sum[ibody][3] += dy*f[i][2] - dz*f[i][1]; + sum[ibody][4] += dz*f[i][0] - dx*f[i][2]; + sum[ibody][5] += dx*f[i][1] - dy*f[i][0]; + } + + // extended particles add their torque to torque of body + + if (extended) { + double **torque_one = atom->torque; + + for (i = 0; i < nlocal; i++) { + if (body[i] < 0) continue; + ibody = body[i]; + + if (eflags[i] & TORQUE) { + sum[ibody][3] += torque_one[i][0]; + sum[ibody][4] += torque_one[i][1]; + sum[ibody][5] += torque_one[i][2]; + } + } + } + + MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world); + + // update vcm and angmom + // include Langevin thermostat forces + // fflag,tflag = 0 for some dimensions in 2d + + 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) +{ + 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() +{ + tagint 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 + + tagint *image = atom->image; + int nlocal = atom->nlocal; + + int ibody; + tagint 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 degrees-of-freedom removed by fix_rigid for atoms in igroup +------------------------------------------------------------------------- */ + +int FixRigid::dof(int igroup) +{ + int groupbit = group->bitmask[igroup]; + + // 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] & groupbit) { + if (extended && eflags[i]) 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,itype; + 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]; + + tagint *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 ibody,itype; + int xbox,ybox,zbox; + double dx,dy,dz; + 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; + tagint *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; + 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; + 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]); + } + } + } +} + +/* ---------------------------------------------------------------------- + one-time initialization of rigid body attributes + extended flags, masstotal, center-of-mass + Cartesian and diagonalized inertia tensor + read per-body attributes from infile if specified +------------------------------------------------------------------------- */ + +void FixRigid::setup_bodies() +{ + int i,itype,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; + tagint *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 + + for (ibody = 0; ibody < nbody; ibody++) + imagebody[ibody] = ((tagint) IMGMAX << IMG2BITS) | + ((tagint) 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,rad; + + 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); +} + +/* ---------------------------------------------------------------------- + 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 i,j,m,nchunk,id; + 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); + + char *buffer = new char[CHUNK*MAXLINE]; + char **values = new char*[ATTRIBUTE_PERBODY]; + + int nread = 0; + while (nread < nlines) { + if (nlines-nread > CHUNK) nchunk = CHUNK; + else nchunk = nlines-nread; + if (me == 0) { + char *eof; + m = 0; + for (i = 0; i < nchunk; i++) { + eof = fgets(&buffer[m],MAXLINE,fp); + if (eof == NULL) error->one(FLERR,"Unexpected end of fix rigid file"); + m += strlen(&buffer[m]); + } + m++; + } + MPI_Bcast(&m,1,MPI_INT,0,world); + MPI_Bcast(buffer,m,MPI_CHAR,0,world); + + 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][5] = atof(values[8]); + array[id][4] = atof(values[9]); + array[id][3] = atof(values[10]); + } + + buf = next + 1; + } + + nread += nchunk; + } + + if (me == 0) fclose(fp); + + delete [] buffer; + delete [] values; +} + +/* ---------------------------------------------------------------------- + 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); + 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) +{ + 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; +} + +/* ---------------------------------------------------------------------- + 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; +} + +/* ---------------------------------------------------------------------- + extract thermostat properties +------------------------------------------------------------------------- */ + +void *FixRigid::extract(const char *str, int &dim) +{ + dim=0; + if (strcmp(str,"t_target") == 0) { + return &t_target; + } + return NULL; +} + +/* ---------------------------------------------------------------------- + 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 new file mode 100644 index 0000000000..9c0dcf7d2e --- /dev/null +++ b/src/RIGID/fix_rigid.h @@ -0,0 +1,234 @@ +/* -*- 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: + // public so that granular pair styles can access them + int *body; // which body each atom is part of (-1 if none) + double *masstotal; // total mass of each rigid body + + 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); + virtual double compute_scalar(); + virtual int modify_param(int, char **) {return 0;} + + double memory_usage(); + void grow_arrays(int); + void copy_arrays(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(); + virtual void *extract(const char*,int &); + double compute_array(int, int); + + protected: + int me,nprocs; + double dtv,dtf,dtq; + double *step_respa; + int triclinic; + double MINUSPI,TWOPI; + + int rstyle; // SINGLE,MOLECULE,GROUP + int firstflag; // 1 for first-time setup of rigid bodies + char *infile; // file to read rigid body attributes from + + int dimension; // # of dimensions + int nbody; // # of rigid bodies + int *nrigid; // # of atoms in each rigid body + int *mol2body; // convert mol-ID to rigid body index + int maxmol; // size of mol2body = max mol-ID + + 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 + tagint *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 **displace; // displacement of each atom in body coords + + 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); + void set_xv(); + void set_v(); + void setup_bodies(); + 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: 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: 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: 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 + +You cannot set image flags for non-periodic dimensions. + +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: 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 or an existing ID of rigid bodies +that are defined by the fix rigid command. + +*/ diff --git a/src/RIGID/fix_rigid_nh.cpp b/src/RIGID/fix_rigid_nh.cpp new file mode 100644 index 0000000000..dbedfe1530 --- /dev/null +++ b/src/RIGID/fix_rigid_nh.cpp @@ -0,0 +1,1436 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Tony Sheh (U Michigan), Trung Dac Nguyen (U Michigan) + references: Kamberaj et al., J. Chem. Phys. 122, 224114 (2005) + Miller et al., J Chem Phys. 116, 8649-8659 (2002) +------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdio.h" +#include "string.h" +#include "fix_rigid_nh.h" +#include "math_extra.h" +#include "atom.h" +#include "compute.h" +#include "domain.h" +#include "update.h" +#include "modify.h" +#include "fix_deform.h" +#include "group.h" +#include "comm.h" +#include "force.h" +#include "kspace.h" +#include "output.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; +using namespace FixConst; + +enum{NONE,XYZ,XY,YZ,XZ}; // same as in FixRigid +enum{ISO,ANISO,TRICLINIC}; // same as in FixRigid + +#define EPSILON 1.0e-7 + +/* ---------------------------------------------------------------------- */ + +FixRigidNH::FixRigidNH(LAMMPS *lmp, int narg, char **arg) : + FixRigid(lmp, narg, arg) +{ + // error checks: could be moved up to FixRigid + + if ((p_flag[0] == 1 && p_period[0] <= 0.0) || + (p_flag[1] == 1 && p_period[1] <= 0.0) || + (p_flag[2] == 1 && p_period[2] <= 0.0)) + error->all(FLERR,"Fix rigid npt/nph period must be > 0.0"); + if (domain->nonperiodic) + error->all(FLERR, + "Cannot use fix rigid npt/nph with non-periodic dimension"); + + if (dimension == 2 && p_flag[2]) + error->all(FLERR,"Invalid fix rigid npt/nph command for a 2d simulation"); + if (dimension == 2 && (pcouple == YZ || pcouple == XZ)) + error->all(FLERR,"Invalid fix rigid npt/nph command for a 2d simulation"); + + if (pcouple == XYZ && (p_flag[0] == 0 || p_flag[1] == 0)) + error->all(FLERR,"Invalid fix rigid npt/nph command pressure settings"); + if (pcouple == XYZ && dimension == 3 && p_flag[2] == 0) + error->all(FLERR,"Invalid fix rigid npt/nph command pressure settings"); + if (pcouple == XY && (p_flag[0] == 0 || p_flag[1] == 0)) + error->all(FLERR,"Invalid fix rigid npt/nph command pressure settings"); + if (pcouple == YZ && (p_flag[1] == 0 || p_flag[2] == 0)) + error->all(FLERR,"Invalid fix rigid npt/nph command pressure settings"); + if (pcouple == XZ && (p_flag[0] == 0 || p_flag[2] == 0)) + error->all(FLERR,"Invalid fix rigid npt/nph command pressure settings"); + + // require periodicity in tensile dimension + + if (p_flag[0] && domain->xperiodic == 0) + error->all(FLERR, + "Cannot use fix rigid npt/nph on a non-periodic dimension"); + if (p_flag[1] && domain->yperiodic == 0) + error->all(FLERR, + "Cannot use fix rigid npt/nph on a non-periodic dimension"); + if (p_flag[2] && domain->zperiodic == 0) + error->all(FLERR, + "Cannot use fix rigid npt/nph on a non-periodic dimension"); + + if (pcouple == XYZ && dimension == 3 && + (p_start[0] != p_start[1] || p_start[0] != p_start[2] || + p_stop[0] != p_stop[1] || p_stop[0] != p_stop[2] || + p_period[0] != p_period[1] || p_period[0] != p_period[2])) + error->all(FLERR,"Invalid fix rigid npt/nph pressure settings"); + if (pcouple == XYZ && dimension == 2 && + (p_start[0] != p_start[1] || p_stop[0] != p_stop[1] || + p_period[0] != p_period[1])) + error->all(FLERR,"Invalid fix rigid npt/nph pressure settings"); + if (pcouple == XY && + (p_start[0] != p_start[1] || p_stop[0] != p_stop[1] || + p_period[0] != p_period[1])) + error->all(FLERR,"Invalid fix rigid npt/nph pressure settings"); + if (pcouple == YZ && + (p_start[1] != p_start[2] || p_stop[1] != p_stop[2] || + p_period[1] != p_period[2])) + error->all(FLERR,"Invalid fix rigid npt/nph pressure settings"); + if (pcouple == XZ && + (p_start[0] != p_start[2] || p_stop[0] != p_stop[2] || + p_period[0] != p_period[2])) + error->all(FLERR,"Invalid fix rigid npt/nph pressure settings"); + + if ((tstat_flag && t_period <= 0.0) || + (p_flag[0] && p_period[0] <= 0.0) || + (p_flag[1] && p_period[1] <= 0.0) || + (p_flag[2] && p_period[2] <= 0.0)) + error->all(FLERR,"Fix rigid nvt/npt/nph damping parameters must be > 0.0"); + + // memory allocation and initialization + + memory->create(conjqm,nbody,4,"rigid_nh:conjqm"); + if (tstat_flag || pstat_flag) { + allocate_chain(); + allocate_order(); + } + + if (tstat_flag) { + eta_t[0] = eta_r[0] = 0.0; + eta_dot_t[0] = eta_dot_r[0] = 0.0; + f_eta_t[0] = f_eta_r[0] = 0.0; + + for (int i = 1; i < t_chain; i++) { + eta_t[i] = eta_r[i] = 0.0; + eta_dot_t[i] = eta_dot_r[i] = 0.0; + } + } + + if (pstat_flag) { + epsilon_dot[0] = epsilon_dot[1] = epsilon_dot[2] = 0.0; + eta_b[0] = eta_dot_b[0] = f_eta_b[0] = 0.0; + for (int i = 1; i < p_chain; i++) + eta_b[i] = eta_dot_b[i] = 0.0; + } + + // rigid body pointers + + nrigidfix = 0; + rfix = NULL; + + vol0 = 0.0; + t0 = 1.0; + + tcomputeflag = 0; + pcomputeflag = 0; +} + +/* ---------------------------------------------------------------------- */ + +FixRigidNH::~FixRigidNH() +{ + memory->destroy(conjqm); + if (tstat_flag || pstat_flag) { + deallocate_chain(); + deallocate_order(); + } + + if (rfix) delete [] rfix; + + if (tcomputeflag) { + modify->delete_compute(id_temp); + delete [] id_temp; + } + + // delete pressure if fix created it + + if (pstat_flag) { + if (pcomputeflag) modify->delete_compute(id_press); + delete [] id_press; + } +} + +/* ---------------------------------------------------------------------- */ + +int FixRigidNH::setmask() +{ + int mask = 0; + mask = FixRigid::setmask(); + if (tstat_flag || pstat_flag) mask |= THERMO_ENERGY; + + return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixRigidNH::init() +{ + FixRigid::init(); + + // recheck that dilate group has not been deleted + + if (allremap == 0) { + int idilate = group->find(id_dilate); + if (idilate == -1) + error->all(FLERR,"Fix rigid npt/nph dilate group ID does not exist"); + dilate_group_bit = group->bitmask[idilate]; + } + + // initialize thermostats + // set timesteps, constants + // store Yoshida-Suzuki integrator parameters + + dtv = update->dt; + dtf = 0.5 * update->dt * force->ftm2v; + dtq = 0.5 * update->dt; + + boltz = force->boltz; + nktv2p = force->nktv2p; + mvv2e = force->mvv2e; + + if (force->kspace) kspace_flag = 1; + else kspace_flag = 0; + + nf_t = nf_r = dimension * nbody; + for (int ibody = 0; ibody < nbody; ibody++) + for (int k = 0; k < domain->dimension; k++) + if (fabs(inertia[ibody][k]) < EPSILON) nf_r--; + + // see Table 1 in Kamberaj et al + + if (tstat_flag || pstat_flag) { + if (t_order == 3) { + w[0] = 1.0 / (2.0 - pow(2.0, 1.0/3.0)); + w[1] = 1.0 - 2.0*w[0]; + w[2] = w[0]; + } else if (t_order == 5) { + w[0] = 1.0 / (4.0 - pow(4.0, 1.0/3.0)); + w[1] = w[0]; + w[2] = 1.0 - 4.0 * w[0]; + w[3] = w[0]; + w[4] = w[0]; + } + } + + g_f = nf_t + nf_r; + onednft = 1.0 + (double)(dimension) / (double)g_f; + onednfr = (double) (dimension) / (double)g_f; + + int icompute; + if (tcomputeflag) { + icompute = modify->find_compute(id_temp); + if (icompute < 0) + error->all(FLERR,"Temp ID for fix rigid npt/nph does not exist"); + temperature = modify->compute[icompute]; + } + + if (pstat_flag) { + if (domain->triclinic) + error->all(FLERR,"fix rigid npt/nph does not yet allow triclinic box"); + + // ensure no conflict with fix deform + + for (int i = 0; i < modify->nfix; i++) + if (strcmp(modify->fix[i]->style,"deform") == 0) { + int *dimflag = ((FixDeform *) modify->fix[i])->dimflag; + if ((p_flag[0] && dimflag[0]) || (p_flag[1] && dimflag[1]) || + (p_flag[2] && dimflag[2])) + error->all(FLERR,"Cannot use fix rigid npt/nph and fix deform on " + "same component of stress tensor"); + } + + // set frequency + + p_freq_max = 0.0; + p_freq_max = MAX(p_freq[0],p_freq[1]); + p_freq_max = MAX(p_freq_max,p_freq[2]); + + // tally the number of dimensions that are barostatted + // set initial volume and reference cell, if not already done + + pdim = p_flag[0] + p_flag[1] + p_flag[2]; + if (vol0 == 0.0) { + if (dimension == 2) vol0 = domain->xprd * domain->yprd; + else vol0 = domain->xprd * domain->yprd * domain->zprd; + } + + // set pressure compute ptr + + icompute = modify->find_compute(id_press); + if (icompute < 0) + error->all(FLERR,"Press ID for fix rigid npt/nph does not exist"); + pressure = modify->compute[icompute]; + + // detect if any rigid fixes exist so rigid bodies move on remap + // rfix[] = indices to each fix rigid + // this will include self + + if (rfix) delete [] rfix; + nrigidfix = 0; + rfix = NULL; + + for (int i = 0; i < modify->nfix; i++) + if (modify->fix[i]->rigid_flag) nrigidfix++; + if (nrigidfix) { + rfix = new int[nrigidfix]; + nrigidfix = 0; + for (int i = 0; i < modify->nfix; i++) + if (modify->fix[i]->rigid_flag) rfix[nrigidfix++] = i; + } + } +} + +/* ---------------------------------------------------------------------- */ + +void FixRigidNH::setup(int vflag) +{ + FixRigid::setup(vflag); + + double mbody[3]; + akin_t = akin_r = 0.0; + for (int ibody = 0; ibody < nbody; ibody++) { + MathExtra::transpose_matvec(ex_space[ibody],ey_space[ibody],ez_space[ibody], + angmom[ibody],mbody); + MathExtra::quatvec(quat[ibody],mbody,conjqm[ibody]); + conjqm[ibody][0] *= 2.0; + conjqm[ibody][1] *= 2.0; + conjqm[ibody][2] *= 2.0; + conjqm[ibody][3] *= 2.0; + + if (tstat_flag || pstat_flag) { + akin_t += masstotal[ibody]*(vcm[ibody][0]*vcm[ibody][0] + + vcm[ibody][1]*vcm[ibody][1] + vcm[ibody][2]*vcm[ibody][2]); + akin_r += angmom[ibody][0]*omega[ibody][0] + + angmom[ibody][1]*omega[ibody][1] + angmom[ibody][2]*omega[ibody][2]; + } + } + + // compute target temperature + + if (tstat_flag) compute_temp_target(); + else if (pstat_flag) { + t0 = temperature->compute_scalar(); + if (t0 == 0.0) { + if (strcmp(update->unit_style,"lj") == 0) t0 = 1.0; + else t0 = 300.0; + } + t_target = t0; + } + + // compute target pressure + // compute current pressure + // trigger virial computation on next timestep + + if (pstat_flag) { + compute_press_target(); + + temperature->compute_scalar(); + if (pstyle == ISO) pressure->compute_scalar(); + else pressure->compute_vector(); + couple(); + pressure->addstep(update->ntimestep+1); + } + + // initialize thermostat/barostat settings + + double kt, t_mass, tb_mass; + kt = boltz * t_target; + + if (tstat_flag) { + t_mass = kt / (t_freq*t_freq); + q_t[0] = nf_t * t_mass; + q_r[0] = nf_r * t_mass; + for (int i = 1; i < t_chain; i++) + q_t[i] = q_r[i] = t_mass; + + for (int i = 1; i < t_chain; i++) { + f_eta_t[i] = (q_t[i-1] * eta_dot_t[i-1] * eta_dot_t[i-1] - kt)/q_t[i]; + f_eta_r[i] = (q_r[i-1] * eta_dot_r[i-1] * eta_dot_r[i-1] - kt)/q_r[i]; + } + } + + // initial forces on barostat thermostat variables + + if (pstat_flag) { + for (int i = 0; i < 3; i++) + if (p_flag[i]) { + epsilon_mass[i] = (g_f + dimension) * kt / (p_freq[i]*p_freq[i]); + epsilon[i] = log(vol0)/dimension; + } + + tb_mass = kt / (p_freq_max * p_freq_max); + q_b[0] = dimension * dimension * tb_mass; + for (int i = 1; i < p_chain; i++) { + q_b[i] = tb_mass; + f_eta_b[i] = (q_b[i] * eta_dot_b[i-1] * eta_dot_b[i-1] - kt)/q_b[i]; + } + } + + // update order/timestep dependent coefficients + + if (tstat_flag || pstat_flag) { + for (int i = 0; i < t_order; i++) { + wdti1[i] = w[i] * dtv / t_iter; + wdti2[i] = wdti1[i] / 2.0; + wdti4[i] = wdti1[i] / 4.0; + } + } +} + +/* ---------------------------------------------------------------------- + perform preforce velocity Verlet integration + see Kamberaj paper for step references +------------------------------------------------------------------------- */ + +void FixRigidNH::initial_integrate(int vflag) +{ + double tmp,scale_r,scale_t[3],scale_v[3]; + double dtfm,mbody[3],tbody[3],fquat[4]; + double dtf2 = dtf * 2.0; + + // compute target temperature + // update thermostat chains coupled to particles + + if (tstat_flag) { + compute_temp_target(); + nhc_temp_integrate(); + } + + // compute target pressure + // update epsilon dot + // update thermostat coupled to barostat + + if (pstat_flag) { + nhc_press_integrate(); + + if (pstyle == ISO) { + temperature->compute_scalar(); + pressure->compute_scalar(); + } else { + temperature->compute_vector(); + pressure->compute_vector(); + } + couple(); + pressure->addstep(update->ntimestep+1); + + compute_press_target(); + nh_epsilon_dot(); + } + + // compute scale variables + + scale_t[0] = scale_t[1] = scale_t[2] = 1.0; + scale_v[0] = scale_v[1] = scale_v[2] = 1.0; + scale_r = 1.0; + + if (tstat_flag) { + akin_t = akin_r = 0.0; + tmp = exp(-1.0 * dtq * eta_dot_t[0]); + scale_t[0] = scale_t[1] = scale_t[2] = tmp; + tmp = exp(-1.0 * dtq * eta_dot_r[0]); + scale_r = tmp; + } + + if (pstat_flag) { + akin_t = akin_r = 0.0; + scale_t[0] *= exp(-dtq * (epsilon_dot[0] + mtk_term2)); + scale_t[1] *= exp(-dtq * (epsilon_dot[1] + mtk_term2)); + scale_t[2] *= exp(-dtq * (epsilon_dot[2] + mtk_term2)); + scale_r *= exp(-dtq * (pdim * mtk_term2)); + + tmp = dtq * epsilon_dot[0]; + scale_v[0] = dtv * exp(tmp) * maclaurin_series(tmp); + tmp = dtq * epsilon_dot[1]; + scale_v[1] = dtv * exp(tmp) * maclaurin_series(tmp); + tmp = dtq * epsilon_dot[2]; + scale_v[2] = dtv * exp(tmp) * maclaurin_series(tmp); + } + + // update xcm, vcm, quat, conjqm and angmom + + for (int ibody = 0; ibody < nbody; ibody++) { + + // step 1.1 - update vcm by 1/2 step + + dtfm = dtf / masstotal[ibody]; + vcm[ibody][0] += dtfm * fcm[ibody][0] * fflag[ibody][0]; + vcm[ibody][1] += dtfm * fcm[ibody][1] * fflag[ibody][1]; + vcm[ibody][2] += dtfm * fcm[ibody][2] * fflag[ibody][2]; + + if (tstat_flag || pstat_flag) { + vcm[ibody][0] *= scale_t[0]; + vcm[ibody][1] *= scale_t[1]; + vcm[ibody][2] *= scale_t[2]; + + tmp = vcm[ibody][0]*vcm[ibody][0] + vcm[ibody][1]*vcm[ibody][1] + + vcm[ibody][2]*vcm[ibody][2]; + akin_t += masstotal[ibody]*tmp; + } + + // step 1.2 - update xcm by full step + + if (!pstat_flag) { + xcm[ibody][0] += dtv * vcm[ibody][0]; + xcm[ibody][1] += dtv * vcm[ibody][1]; + xcm[ibody][2] += dtv * vcm[ibody][2]; + } else { + xcm[ibody][0] += scale_v[0] * vcm[ibody][0]; + xcm[ibody][1] += scale_v[1] * vcm[ibody][1]; + xcm[ibody][2] += scale_v[2] * vcm[ibody][2]; + } + + // step 1.3 - apply torque (body coords) to quaternion momentum + + torque[ibody][0] *= tflag[ibody][0]; + torque[ibody][1] *= tflag[ibody][1]; + torque[ibody][2] *= tflag[ibody][2]; + + MathExtra::transpose_matvec(ex_space[ibody],ey_space[ibody],ez_space[ibody], + torque[ibody],tbody); + MathExtra::quatvec(quat[ibody],tbody,fquat); + + conjqm[ibody][0] += dtf2 * fquat[0]; + conjqm[ibody][1] += dtf2 * fquat[1]; + conjqm[ibody][2] += dtf2 * fquat[2]; + conjqm[ibody][3] += dtf2 * fquat[3]; + + if (tstat_flag || pstat_flag) { + conjqm[ibody][0] *= scale_r; + conjqm[ibody][1] *= scale_r; + conjqm[ibody][2] *= scale_r; + conjqm[ibody][3] *= scale_r; + } + + // step 1.4 to 1.13 - use no_squish rotate to update p and q + + no_squish_rotate(3,conjqm[ibody],quat[ibody],inertia[ibody],dtq); + no_squish_rotate(2,conjqm[ibody],quat[ibody],inertia[ibody],dtq); + no_squish_rotate(1,conjqm[ibody],quat[ibody],inertia[ibody],dtv); + no_squish_rotate(2,conjqm[ibody],quat[ibody],inertia[ibody],dtq); + no_squish_rotate(3,conjqm[ibody],quat[ibody],inertia[ibody],dtq); + + // update exyz_space + // transform p back to angmom + // update angular velocity + + MathExtra::q_to_exyz(quat[ibody],ex_space[ibody],ey_space[ibody], + ez_space[ibody]); + MathExtra::invquatvec(quat[ibody],conjqm[ibody],mbody); + MathExtra::matvec(ex_space[ibody],ey_space[ibody],ez_space[ibody], + mbody,angmom[ibody]); + + angmom[ibody][0] *= 0.5; + angmom[ibody][1] *= 0.5; + angmom[ibody][2] *= 0.5; + + MathExtra::angmom_to_omega(angmom[ibody],ex_space[ibody],ey_space[ibody], + ez_space[ibody],inertia[ibody],omega[ibody]); + + if (tstat_flag || pstat_flag) { + akin_r += angmom[ibody][0]*omega[ibody][0] + + angmom[ibody][1]*omega[ibody][1] + angmom[ibody][2]*omega[ibody][2]; + } + } + + // virial setup before call to set_xv + + if (vflag) v_setup(vflag); + else evflag = 0; + + // remap simulation box by 1/2 step + + if (pstat_flag) remap(); + + // set coords/orient and velocity/rotation of atoms in rigid bodies + // from quarternion and omega + + set_xv(); + + // remap simulation box by full step + // redo KSpace coeffs since volume has changed + + if (pstat_flag) { + remap(); + if (kspace_flag) force->kspace->setup(); + } +} + +/* ---------------------------------------------------------------------- */ + +void FixRigidNH::final_integrate() +{ + int i,ibody; + double tmp,scale_t[3],scale_r; + double dtfm,xy,xz,yz; + double mbody[3],tbody[3],fquat[4]; + double dtf2 = dtf * 2.0; + + // compute scale variables + + scale_t[0] = scale_t[1] = scale_t[2] = 1.0; + scale_r = 1.0; + + if (tstat_flag) { + tmp = exp(-1.0 * dtq * eta_dot_t[0]); + scale_t[0] = scale_t[1] = scale_t[2] = tmp; + scale_r = exp(-1.0 * dtq * eta_dot_r[0]); + } + + if (pstat_flag) { + scale_t[0] *= exp(-dtq * (epsilon_dot[0] + mtk_term2)); + scale_t[1] *= exp(-dtq * (epsilon_dot[1] + mtk_term2)); + scale_t[2] *= exp(-dtq * (epsilon_dot[2] + mtk_term2)); + scale_r *= exp(-dtq * (pdim * mtk_term2)); + + akin_t = akin_r = 0.0; + } + + // sum over atoms to get force and torque on rigid body + + int *image = atom->image; + double **x = atom->x; + double **f = atom->f; + int nlocal = atom->nlocal; + + double xprd = domain->xprd; + double yprd = domain->yprd; + double zprd = domain->zprd; + if (triclinic) { + xy = domain->xy; + xz = domain->xz; + yz = domain->yz; + } + + int xbox,ybox,zbox; + double xunwrap,yunwrap,zunwrap,dx,dy,dz; + for (ibody = 0; ibody < nbody; ibody++) + for (i = 0; i < 6; i++) sum[ibody][i] = 0.0; + + for (i = 0; i < nlocal; i++) { + if (body[i] < 0) continue; + ibody = body[i]; + + sum[ibody][0] += f[i][0]; + sum[ibody][1] += f[i][1]; + sum[ibody][2] += f[i][2]; + + xbox = (image[i] & IMGMASK) - IMGMAX; + ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; + zbox = (image[i] >> IMG2BITS) - IMGMAX; + + if (triclinic == 0) { + xunwrap = x[i][0] + xbox*xprd; + yunwrap = x[i][1] + ybox*yprd; + zunwrap = x[i][2] + zbox*zprd; + } else { + xunwrap = x[i][0] + xbox*xprd + ybox*xy + zbox*xz; + yunwrap = x[i][1] + ybox*yprd + zbox*yz; + zunwrap = x[i][2] + zbox*zprd; + } + + dx = xunwrap - xcm[ibody][0]; + dy = yunwrap - xcm[ibody][1]; + dz = zunwrap - xcm[ibody][2]; + + sum[ibody][3] += dy*f[i][2] - dz*f[i][1]; + sum[ibody][4] += dz*f[i][0] - dx*f[i][2]; + sum[ibody][5] += dx*f[i][1] - dy*f[i][0]; + } + + // extended particles add their torque to torque of body + + if (extended) { + double **torque_one = atom->torque; + + for (i = 0; i < nlocal; i++) { + if (body[i] < 0) continue; + ibody = body[i]; + + if (eflags[i] & TORQUE) { + sum[ibody][3] += torque_one[i][0]; + sum[ibody][4] += torque_one[i][1]; + sum[ibody][5] += torque_one[i][2]; + } + } + } + + MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world); + + // update vcm and angmom + // include Langevin thermostat forces + // fflag,tflag = 0 for some dimensions in 2d + + for (ibody = 0; ibody < nbody; ibody++) { + fcm[ibody][0] = all[ibody][0] + langextra[ibody][0]; + fcm[ibody][1] = all[ibody][1] + langextra[ibody][1]; + fcm[ibody][2] = all[ibody][2] + langextra[ibody][2]; + torque[ibody][0] = all[ibody][3] + langextra[ibody][3]; + torque[ibody][1] = all[ibody][4] + langextra[ibody][4]; + torque[ibody][2] = all[ibody][5] + langextra[ibody][5]; + + // update vcm by 1/2 step + + dtfm = dtf / masstotal[ibody]; + if (tstat_flag || pstat_flag) { + vcm[ibody][0] *= scale_t[0]; + vcm[ibody][1] *= scale_t[1]; + vcm[ibody][2] *= scale_t[2]; + } + + vcm[ibody][0] += dtfm * fcm[ibody][0] * fflag[ibody][0]; + vcm[ibody][1] += dtfm * fcm[ibody][1] * fflag[ibody][1]; + vcm[ibody][2] += dtfm * fcm[ibody][2] * fflag[ibody][2]; + + if (pstat_flag) { + tmp = vcm[ibody][0]*vcm[ibody][0] + vcm[ibody][1]*vcm[ibody][1] + + vcm[ibody][2]*vcm[ibody][2]; + akin_t += masstotal[ibody]*tmp; + } + + // update conjqm, then transform to angmom, set velocity again + // virial is already setup from initial_integrate + + torque[ibody][0] *= tflag[ibody][0]; + torque[ibody][1] *= tflag[ibody][1]; + torque[ibody][2] *= tflag[ibody][2]; + + MathExtra::transpose_matvec(ex_space[ibody],ey_space[ibody], + ez_space[ibody],torque[ibody],tbody); + MathExtra::quatvec(quat[ibody],tbody,fquat); + + if (tstat_flag || pstat_flag) { + conjqm[ibody][0] = scale_r * conjqm[ibody][0] + dtf2 * fquat[0]; + conjqm[ibody][1] = scale_r * conjqm[ibody][1] + dtf2 * fquat[1]; + conjqm[ibody][2] = scale_r * conjqm[ibody][2] + dtf2 * fquat[2]; + conjqm[ibody][3] = scale_r * conjqm[ibody][3] + dtf2 * fquat[3]; + } else { + conjqm[ibody][0] += dtf2 * fquat[0]; + conjqm[ibody][1] += dtf2 * fquat[1]; + conjqm[ibody][2] += dtf2 * fquat[2]; + conjqm[ibody][3] += dtf2 * fquat[3]; + } + + MathExtra::invquatvec(quat[ibody],conjqm[ibody],mbody); + MathExtra::matvec(ex_space[ibody],ey_space[ibody],ez_space[ibody], + mbody,angmom[ibody]); + + angmom[ibody][0] *= 0.5; + angmom[ibody][1] *= 0.5; + angmom[ibody][2] *= 0.5; + + MathExtra::angmom_to_omega(angmom[ibody],ex_space[ibody],ey_space[ibody], + ez_space[ibody],inertia[ibody],omega[ibody]); + + if (pstat_flag) { + akin_r += angmom[ibody][0]*omega[ibody][0] + + angmom[ibody][1]*omega[ibody][1] + + angmom[ibody][2]*omega[ibody][2]; + } + } + + // set velocity/rotation of atoms in rigid bodies + // virial is already setup from initial_integrate + + set_v(); + + // compute temperature and pressure tensor + // couple to compute current pressure components + // trigger virial computation on next timestep + + if (tcomputeflag) t_current = temperature->compute_scalar(); + if (pstat_flag) { + if (pstyle == ISO) pressure->compute_scalar(); + else pressure->compute_vector(); + couple(); + pressure->addstep(update->ntimestep+1); + } + + if (pstat_flag) nh_epsilon_dot(); + + // update eta_dot_t and eta_dot_r + // update eta_dot_b + + if (tstat_flag) nhc_temp_integrate(); + if (pstat_flag) nhc_press_integrate(); +} + +/* ---------------------------------------------------------------------- */ + +void FixRigidNH::nhc_temp_integrate() +{ + int i,j,k; + double kt,gfkt_t,gfkt_r,tmp,ms,s,s2; + + kt = boltz * t_target; + gfkt_t = nf_t * kt; + gfkt_r = nf_r * kt; + + // update thermostat masses + + double t_mass = boltz * t_target / (t_freq * t_freq); + q_t[0] = nf_t * t_mass; + q_r[0] = nf_r * t_mass; + for (i = 1; i < t_chain; i++) + q_t[i] = q_r[i] = t_mass; + + // update force of thermostats coupled to particles + + f_eta_t[0] = (akin_t * mvv2e - gfkt_t) / q_t[0]; + f_eta_r[0] = (akin_r * mvv2e - gfkt_r) / q_r[0]; + + // multiple timestep iteration + + for (i = 0; i < t_iter; i++) { + for (j = 0; j < t_order; j++) { + + // update thermostat velocities half step + + eta_dot_t[t_chain-1] += wdti2[j] * f_eta_t[t_chain-1]; + eta_dot_r[t_chain-1] += wdti2[j] * f_eta_r[t_chain-1]; + + for (k = 1; k < t_chain; k++) { + tmp = wdti4[j] * eta_dot_t[t_chain-k]; + ms = maclaurin_series(tmp); + s = exp(-1.0 * tmp); + s2 = s * s; + eta_dot_t[t_chain-k-1] = eta_dot_t[t_chain-k-1] * s2 + + wdti2[j] * f_eta_t[t_chain-k-1] * s * ms; + + tmp = wdti4[j] * eta_dot_r[t_chain-k]; + ms = maclaurin_series(tmp); + s = exp(-1.0 * tmp); + s2 = s * s; + eta_dot_r[t_chain-k-1] = eta_dot_r[t_chain-k-1] * s2 + + wdti2[j] * f_eta_r[t_chain-k-1] * s * ms; + } + + // update thermostat positions a full step + + for (k = 0; k < t_chain; k++) { + eta_t[k] += wdti1[j] * eta_dot_t[k]; + eta_r[k] += wdti1[j] * eta_dot_r[k]; + } + + // update thermostat forces + + for (k = 1; k < t_chain; k++) { + f_eta_t[k] = q_t[k-1] * eta_dot_t[k-1] * eta_dot_t[k-1] - kt; + f_eta_t[k] /= q_t[k]; + f_eta_r[k] = q_r[k-1] * eta_dot_r[k-1] * eta_dot_r[k-1] - kt; + f_eta_r[k] /= q_r[k]; + } + + // update thermostat velocities a full step + + for (k = 0; k < t_chain-1; k++) { + tmp = wdti4[j] * eta_dot_t[k+1]; + ms = maclaurin_series(tmp); + s = exp(-1.0 * tmp); + s2 = s * s; + eta_dot_t[k] = eta_dot_t[k] * s2 + wdti2[j] * f_eta_t[k] * s * ms; + tmp = q_t[k] * eta_dot_t[k] * eta_dot_t[k] - kt; + f_eta_t[k+1] = tmp / q_t[k+1]; + + tmp = wdti4[j] * eta_dot_r[k+1]; + ms = maclaurin_series(tmp); + s = exp(-1.0 * tmp); + s2 = s * s; + eta_dot_r[k] = eta_dot_r[k] * s2 + wdti2[j] * f_eta_r[k] * s * ms; + tmp = q_r[k] * eta_dot_r[k] * eta_dot_r[k] - kt; + f_eta_r[k+1] = tmp / q_r[k+1]; + } + + eta_dot_t[t_chain-1] += wdti2[j] * f_eta_t[t_chain-1]; + eta_dot_r[t_chain-1] += wdti2[j] * f_eta_r[t_chain-1]; + } + } +} + +/* ---------------------------------------------------------------------- */ + +void FixRigidNH::nhc_press_integrate() +{ + int i,k; + double tmp,s,s2,ms,kecurrent; + double kt = boltz * t_target; + double lkt_press = kt; + + // update thermostat masses + + double tb_mass = kt / (p_freq_max * p_freq_max); + q_b[0] = tb_mass; + for (int i = 1; i < p_chain; i++) { + q_b[i] = tb_mass; + f_eta_b[i] = q_b[i-1] * eta_dot_b[i-1] * eta_dot_b[i-1] - kt; + f_eta_b[i] /= q_b[i]; + } + + // update forces acting on thermostat + + kecurrent = 0.0; + for (i = 0; i < 3; i++) + if (p_flag[i]) { + epsilon_mass[i] = (g_f + dimension) * kt / (p_freq[i] * p_freq[i]); + kecurrent += epsilon_mass[i] * epsilon_dot[i] * epsilon_dot[i]; + } + + f_eta_b[0] = (kecurrent - lkt_press) / q_b[0]; + + // update thermostat velocities a half step + + eta_dot_b[p_chain-1] += 0.5 * dtq * f_eta_b[p_chain-1]; + + for (k = 0; k < p_chain-1; k++) { + tmp = 0.5 * dtq * eta_dot_b[p_chain-k-1]; + ms = maclaurin_series(tmp); + s = exp(-0.5 * tmp); + s2 = s * s; + eta_dot_b[p_chain-k-2] = eta_dot_b[p_chain-k-2] * s2 + + dtq * f_eta_b[p_chain-k-2] * s * ms; + } + + // update thermostat positions + + for (k = 0; k < p_chain; k++) + eta_b[k] += dtv * eta_dot_b[k]; + + // update epsilon dot + + s = exp(-1.0 * dtq * eta_dot_b[0]); + for (i = 0; i < 3; i++) + if (p_flag[i]) epsilon_dot[i] *= s; + + kecurrent = 0.0; + for (i = 0; i < 3; i++) + if (p_flag[i]) + kecurrent += epsilon_mass[i] * epsilon_dot[i] * epsilon_dot[i]; + + f_eta_b[0] = (kecurrent - lkt_press) / q_b[0]; + + // update thermostat velocites a full step + + for (k = 0; k < p_chain-1; k++) { + tmp = 0.5 * dtq * eta_dot_b[k+1]; + ms = maclaurin_series(tmp); + s = exp(-0.5 * tmp); + s2 = s * s; + eta_dot_b[k] = eta_dot_b[k] * s2 + dtq * f_eta_b[k] * s * ms; + tmp = q_b[k] * eta_dot_b[k] * eta_dot_b[k] - kt; + f_eta_b[k+1] = tmp / q_b[k+1]; + } + + eta_dot_b[p_chain-1] += 0.5 * dtq * f_eta_b[p_chain-1]; + +} + +/* ---------------------------------------------------------------------- + compute kinetic energy in the extended Hamiltonian + conserved quantity = sum of returned energy and potential energy +-----------------------------------------------------------------------*/ + +double FixRigidNH::compute_scalar() +{ + int i,k,ibody; + double kt = boltz * t_target; + double energy,ke_t,ke_q,tmp,Pkq[4]; + + // compute the kinetic parts of H_NVE in Kameraj et al (JCP 2005, pp 224114) + + // translational kinetic energy + + ke_t = 0.0; + for (ibody = 0; ibody < nbody; ibody++) + ke_t += 0.5 * masstotal[ibody] * (vcm[ibody][0]*vcm[ibody][0] + + vcm[ibody][1]*vcm[ibody][1] + + vcm[ibody][2]*vcm[ibody][2]); + + // rotational kinetic energy + + ke_q = 0.0; + for (ibody = 0; ibody < nbody; ibody++) { + for (k = 1; k < 4; k++) { + if (k == 1) { + Pkq[0] = -quat[ibody][1]; + Pkq[1] = quat[ibody][0]; + Pkq[2] = quat[ibody][3]; + Pkq[3] = -quat[ibody][2]; + } else if (k == 2) { + Pkq[0] = -quat[ibody][2]; + Pkq[1] = -quat[ibody][3]; + Pkq[2] = quat[ibody][0]; + Pkq[3] = quat[ibody][1]; + } else if (k == 3) { + Pkq[0] = -quat[ibody][3]; + Pkq[1] = quat[ibody][2]; + Pkq[2] = -quat[ibody][1]; + Pkq[3] = quat[ibody][0]; + } + + tmp = conjqm[ibody][0]*Pkq[0] + conjqm[ibody][1]*Pkq[1] + + conjqm[ibody][2]*Pkq[2] + conjqm[ibody][3]*Pkq[3]; + tmp *= tmp; + + if (fabs(inertia[ibody][k-1]) < 1e-6) tmp = 0.0; + else tmp /= (8.0 * inertia[ibody][k-1]); + ke_q += tmp; + } + } + + energy = (ke_t + ke_q) * mvv2e; + + if (tstat_flag) { + + // thermostat chain energy: from equation 12 in Kameraj et al (JCP 2005) + + energy += kt * (nf_t * eta_t[0] + nf_r * eta_r[0]); + + for (i = 1; i < t_chain; i++) + energy += kt * (eta_t[i] + eta_r[i]); + + for (i = 0; i < t_chain; i++) { + energy += 0.5 * q_t[i] * (eta_dot_t[i] * eta_dot_t[i]); + energy += 0.5 * q_r[i] * (eta_dot_r[i] * eta_dot_r[i]); + } + } + + if (pstat_flag) { + + // using equation 22 in Kameraj et al for H_NPT + + for (i = 0; i < 3; i++) + energy += 0.5 * epsilon_mass[i] * epsilon_dot[i] * epsilon_dot[i]; + + double vol; + if (dimension == 2) vol = domain->xprd * domain->yprd; + else vol = domain->xprd * domain->yprd * domain->zprd; + + double p0 = (p_target[0] + p_target[1] + p_target[2]) / 3.0; + energy += p0 * vol / nktv2p; + + for (i = 0; i < p_chain; i++) { + energy += kt * eta_b[i]; + energy += 0.5 * q_b[i] * (eta_dot_b[i] * eta_dot_b[i]); + } + } + + return energy; +} + +/* ---------------------------------------------------------------------- */ + +void FixRigidNH::couple() +{ + double *tensor = pressure->vector; + + if (pstyle == ISO) { + p_current[0] = p_current[1] = p_current[2] = pressure->scalar; + } else if (pcouple == XYZ) { + double ave = 1.0/3.0 * (tensor[0] + tensor[1] + tensor[2]); + p_current[0] = p_current[1] = p_current[2] = ave; + } else if (pcouple == XY) { + double ave = 0.5 * (tensor[0] + tensor[1]); + p_current[0] = p_current[1] = ave; + p_current[2] = tensor[2]; + } else if (pcouple == YZ) { + double ave = 0.5 * (tensor[1] + tensor[2]); + p_current[1] = p_current[2] = ave; + p_current[0] = tensor[0]; + } else if (pcouple == XZ) { + double ave = 0.5 * (tensor[0] + tensor[2]); + p_current[0] = p_current[2] = ave; + p_current[1] = tensor[1]; + } else { + p_current[0] = tensor[0]; + p_current[1] = tensor[1]; + p_current[2] = tensor[2]; + } +} + +/* ---------------------------------------------------------------------- */ + +void FixRigidNH::remap() +{ + int i; + double oldlo,oldhi,ctr,expfac; + + double **x = atom->x; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + // epsilon is not used, except for book-keeping + + for (i = 0; i < 3; i++) epsilon[i] += dtq * epsilon_dot[i]; + + // convert pertinent atoms and rigid bodies to lamda coords + + if (allremap) domain->x2lamda(nlocal); + else { + for (i = 0; i < nlocal; i++) + if (mask[i] & dilate_group_bit) + domain->x2lamda(x[i],x[i]); + } + + if (nrigid) + for (i = 0; i < nrigidfix; i++) + modify->fix[rfix[i]]->deform(0); + + // reset global and local box to new size/shape + + for (i = 0; i < 3; i++) { + if (p_flag[i]) { + oldlo = domain->boxlo[i]; + oldhi = domain->boxhi[i]; + ctr = 0.5 * (oldlo + oldhi); + expfac = exp(dtq * epsilon_dot[i]); + domain->boxlo[i] = (oldlo-ctr)*expfac + ctr; + domain->boxhi[i] = (oldhi-ctr)*expfac + ctr; + } + } + + domain->set_global_box(); + domain->set_local_box(); + + // convert pertinent atoms and rigid bodies back to box coords + + if (allremap) domain->lamda2x(nlocal); + else { + for (i = 0; i < nlocal; i++) + if (mask[i] & dilate_group_bit) + domain->lamda2x(x[i],x[i]); + } + + if (nrigid) + for (i = 0; i< nrigidfix; i++) + modify->fix[rfix[i]]->deform(1); +} + +/* ---------------------------------------------------------------------- + compute target temperature and kinetic energy +-----------------------------------------------------------------------*/ + +void FixRigidNH::compute_temp_target() +{ + double delta = update->ntimestep - update->beginstep; + if (update->endstep > update->beginstep) + delta /= update->endstep - update->beginstep; + else delta = 0.0; + + t_target = t_start + delta * (t_stop-t_start); +} + +/* ---------------------------------------------------------------------- + compute hydrostatic target pressure +-----------------------------------------------------------------------*/ + +void FixRigidNH::compute_press_target() +{ + double delta = update->ntimestep - update->beginstep; + if (update->endstep > update->beginstep) + delta /= update->endstep - update->beginstep; + else delta = 0.0; + + p_hydro = 0.0; + for (int i = 0; i < 3; i++) + if (p_flag[i]) { + p_target[i] = p_start[i] + delta * (p_stop[i]-p_start[i]); + p_hydro += p_target[i]; + } + p_hydro /= pdim; +} + +/* ---------------------------------------------------------------------- + update epsilon_dot +-----------------------------------------------------------------------*/ + +void FixRigidNH::nh_epsilon_dot() +{ + int i; + double volume,scale,f_epsilon; + + if (dimension == 2) volume = domain->xprd*domain->yprd; + else volume = domain->xprd*domain->yprd*domain->zprd; + + // MTK terms + + mtk_term1 = (akin_t + akin_r) * mvv2e / g_f; + + scale = exp(-1.0 * dtq * eta_dot_b[0]); + + for (i = 0; i < 3; i++) + if (p_flag[i]) { + f_epsilon = (p_current[i]-p_hydro)*volume / nktv2p + mtk_term1; + f_epsilon /= epsilon_mass[i]; + epsilon_dot[i] += dtq * f_epsilon; + epsilon_dot[i] *= scale; + } + + mtk_term2 = 0.0; + for (i = 0; i < 3; i++) + if (p_flag[i]) mtk_term2 += epsilon_dot[i]; + mtk_term2 /= g_f; +} + +/* ---------------------------------------------------------------------- + pack entire state of Fix into one write +------------------------------------------------------------------------- */ + +void FixRigidNH::write_restart(FILE *fp) +{ + if (tstat_flag == 0 && pstat_flag == 0) return; + + int nsize = 2; // tstat_flag and pstat_flag + + if (tstat_flag) { + nsize += 1; // t_chain + nsize += 4*t_chain; // eta_t, eta_r, eta_dot_t, eta_dot_r + } + + if (pstat_flag) { + nsize += 7; // p_chain, epsilon(3) and epsilon_dot(3) + nsize += 2*p_chain; + } + + double *list; + memory->create(list,nsize,"rigid_nh:list"); + + int n = 0; + + list[n++] = tstat_flag; + if (tstat_flag) { + list[n++] = t_chain; + for (int i = 0; i < t_chain; i++) { + list[n++] = eta_t[i]; + list[n++] = eta_r[i]; + list[n++] = eta_dot_t[i]; + list[n++] = eta_dot_r[i]; + } + } + + list[n++] = pstat_flag; + if (pstat_flag) { + list[n++] = epsilon[0]; + list[n++] = epsilon[1]; + list[n++] = epsilon[2]; + list[n++] = epsilon_dot[0]; + list[n++] = epsilon_dot[1]; + list[n++] = epsilon_dot[2]; + + list[n++] = p_chain; + for (int i = 0; i < p_chain; i++) { + list[n++] = eta_b[i]; + list[n++] = eta_dot_b[i]; + } + } + + if (comm->me == 0) { + int size = (nsize)*sizeof(double); + fwrite(&size,sizeof(int),1,fp); + fwrite(list,sizeof(double),nsize,fp); + } + + memory->destroy(list); +} + +/* ---------------------------------------------------------------------- + use state info from restart file to restart the Fix +------------------------------------------------------------------------- */ + +void FixRigidNH::restart(char *buf) +{ + int n = 0; + double *list = (double *) buf; + int flag = static_cast (list[n++]); + + if (flag) { + int m = static_cast (list[n++]); + if (tstat_flag && m == t_chain) { + for (int i = 0; i < t_chain; i++) { + eta_t[i] = list[n++]; + eta_r[i] = list[n++]; + eta_dot_t[i] = list[n++]; + eta_dot_r[i] = list[n++]; + } + } else n += 4*m; + } + + flag = static_cast (list[n++]); + if (flag) { + epsilon[0] = list[n++]; + epsilon[1] = list[n++]; + epsilon[2] = list[n++]; + epsilon_dot[0] = list[n++]; + epsilon_dot[1] = list[n++]; + epsilon_dot[2] = list[n++]; + + int m = static_cast (list[n++]); + if (pstat_flag && m == p_chain) { + for (int i = 0; i < p_chain; i++) { + eta_b[i] = list[n++]; + eta_dot_b[i] = list[n++]; + } + } else n += 2*m; + } +} + +/* ---------------------------------------------------------------------- */ + +int FixRigidNH::modify_param(int narg, char **arg) +{ + if (strcmp(arg[0],"temp") == 0) { + if (narg < 2) error->all(FLERR,"Illegal fix_modify command"); + if (!pstat_flag) error->all(FLERR,"Illegal fix_modify command"); + if (tcomputeflag) { + modify->delete_compute(id_temp); + tcomputeflag = 0; + } + delete [] id_temp; + int n = strlen(arg[1]) + 1; + id_temp = new char[n]; + strcpy(id_temp,arg[1]); + + int icompute = modify->find_compute(arg[1]); + if (icompute < 0) + error->all(FLERR,"Could not find fix_modify temperature ID"); + temperature = modify->compute[icompute]; + + if (temperature->tempflag == 0) + error->all(FLERR, + "Fix_modify temperature ID does not compute temperature"); + if (temperature->igroup != 0 && comm->me == 0) + error->warning(FLERR,"Temperature for fix modify is not for group all"); + + // reset id_temp of pressure to new temperature ID + + if (pstat_flag) { + icompute = modify->find_compute(id_press); + if (icompute < 0) + error->all(FLERR,"Pressure ID for fix modify does not exist"); + modify->compute[icompute]->reset_extra_compute_fix(id_temp); + } + + return 2; + + } else if (strcmp(arg[0],"press") == 0) { + if (narg < 2) error->all(FLERR,"Illegal fix_modify command"); + if (!pstat_flag) error->all(FLERR,"Illegal fix_modify command"); + if (pcomputeflag) { + modify->delete_compute(id_press); + pcomputeflag = 0; + } + delete [] id_press; + int n = strlen(arg[1]) + 1; + id_press = new char[n]; + strcpy(id_press,arg[1]); + + int icompute = modify->find_compute(arg[1]); + if (icompute < 0) error->all(FLERR,"Could not find fix_modify pressure ID"); + pressure = modify->compute[icompute]; + + if (pressure->pressflag == 0) + error->all(FLERR,"Fix_modify pressure ID does not compute pressure"); + return 2; + } + + return 0; +} + +/* ---------------------------------------------------------------------- */ + +void FixRigidNH::allocate_chain() +{ + if (tstat_flag) { + q_t = new double[t_chain]; + q_r = new double[t_chain]; + eta_t = new double[t_chain]; + eta_r = new double[t_chain]; + eta_dot_t = new double[t_chain]; + eta_dot_r = new double[t_chain]; + f_eta_t = new double[t_chain]; + f_eta_r = new double[t_chain]; + } + + if (pstat_flag) { + q_b = new double[p_chain]; + eta_b = new double[p_chain]; + eta_dot_b = new double[p_chain]; + f_eta_b = new double[p_chain]; + } +} + +/* ---------------------------------------------------------------------- */ + +void FixRigidNH::reset_target(double t_new) +{ + t_start = t_stop = t_new; +} + +/* ---------------------------------------------------------------------- */ + +void FixRigidNH::allocate_order() +{ + w = new double[t_order]; + wdti1 = new double[t_order]; + wdti2 = new double[t_order]; + wdti4 = new double[t_order]; +} + +/* ---------------------------------------------------------------------- */ + +void FixRigidNH::deallocate_chain() +{ + if (tstat_flag) { + delete [] q_t; + delete [] q_r; + delete [] eta_t; + delete [] eta_r; + delete [] eta_dot_t; + delete [] eta_dot_r; + delete [] f_eta_t; + delete [] f_eta_r; + } + + if (pstat_flag) { + delete [] q_b; + delete [] eta_b; + delete [] eta_dot_b; + delete [] f_eta_b; + } +} + +/* ---------------------------------------------------------------------- */ + +void FixRigidNH::deallocate_order() +{ + delete [] w; + delete [] wdti1; + delete [] wdti2; + delete [] wdti4; +} diff --git a/src/RIGID/fix_rigid_nh.h b/src/RIGID/fix_rigid_nh.h new file mode 100644 index 0000000000..a348c92aae --- /dev/null +++ b/src/RIGID/fix_rigid_nh.h @@ -0,0 +1,178 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifndef LMP_FIX_RIGID_NH_H +#define LMP_FIX_RIGID_NH_H + +#include "fix_rigid.h" + +namespace LAMMPS_NS { + +class FixRigidNH : public FixRigid { + public: + FixRigidNH(class LAMMPS *, int, char **); + virtual ~FixRigidNH(); + virtual int setmask(); + virtual void init(); + virtual void setup(int); + virtual void initial_integrate(int); + virtual void final_integrate(); + virtual double compute_scalar(); + int modify_param(int, char **); + void write_restart(FILE *); + void restart(char *buf); + void reset_target(double); + + protected: + double **conjqm; // conjugate quaternion momentum + double boltz,nktv2p,mvv2e; // boltzman constant, conversion factors + + int nf_t,nf_r; // trans/rot degrees of freedom + double onednft,onednfr; // factors 1 + dimension/trans(rot) degrees of freedom + double *w,*wdti1,*wdti2,*wdti4; // Yoshida-Suzuki coefficients + double *q_t,*q_r; // trans/rot thermostat masses + double *eta_t,*eta_r; // trans/rot thermostat positions + double *eta_dot_t,*eta_dot_r; // trans/rot thermostat velocities + double *f_eta_t,*f_eta_r; // trans/rot thermostat forces + + double epsilon_mass[3], *q_b; // baro/thermo masses + double epsilon[3],*eta_b; // baro/thermo positions + double epsilon_dot[3],*eta_dot_b; // baro/thermo velocities + double *f_eta_b; // thermo forces + double akin_t,akin_r; // translational/rotational kinetic energies + + int kspace_flag; // 1 if KSpace invoked, 0 if not + int nrigidfix; // number of rigid fixes + int *rfix; // indicies of rigid fixes + + double vol0; // reference volume + double t0; // reference temperature + int pdim,g_f; // number of barostatted dims, total DoFs + double p_hydro; // hydrostatic target pressure + double p_freq_max; // maximum barostat frequency + + double mtk_term1,mtk_term2; // Martyna-Tobias-Klein corrections + + double t_current,t_target; + double p_current[3],p_target[3]; + + char *id_temp,*id_press; + class Compute *temperature,*pressure; + int tcomputeflag,pcomputeflag; + + void couple(); + void remap(); + void nhc_temp_integrate(); + void nhc_press_integrate(); + + virtual void compute_temp_target(); + void compute_press_target(); + void nh_epsilon_dot(); + + void allocate_chain(); + void allocate_order(); + void deallocate_chain(); + void deallocate_order(); + + inline double maclaurin_series(double); +}; + +inline double FixRigidNH::maclaurin_series(double x) +{ + double x2,x4; + x2 = x * x; + x4 = x2 * x2; + return (1.0 + (1.0/6.0) * x2 + (1.0/120.0) * x4 + (1.0/5040.0) * x2 * x4 + + (1.0/362880.0) * x4 * x4); +} + +} + +#endif + +/* ERROR/WARNING messages: + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +E: Target temperature for fix rigid nvt/npt cannot be 0.0 + +Self-explanatory. + +E: Invalid fix rigid npt/nph command for a 2d simulation + +Cannot control z dimension in a 2d model. + +E: Fix rigid npt/nph dilate group ID does not exist + +Self-explanatory. + +E: Invalid fix rigid npt/nph command pressure settings + +If multiple dimensions are coupled, those dimensions must be +specified. + +E: Cannot use fix rigid npt/nph on a non-periodic dimension + +When specifying a diagonal pressure component, the dimension must be +periodic. + +E: Invalid fix rigid npt/nph pressure settings + +Settings for coupled dimensions must be the same. + +E: Fix rigid nvt/npt/nph damping parameters must be > 0.0 + +Self-explanatory. + +E: Cannot use fix rigid npt/nph and fix deform on same component of stress tensor + +This would be changing the same box dimension twice. + +E: Temperature ID for fix rigid npt/nph does not exist + +Self-explanatory. + +E: Pressure ID for fix rigid npt/nph does not exist + +Self-explanatory. + +E: Could not find fix_modify temperature ID + +The compute ID for computing temperature does not exist. + +E: Fix_modify temperature ID does not compute temperature + +The compute ID assigned to the fix must compute temperature. + +W: Temperature for fix modify is not for group all + +The temperature compute is being used with a pressure calculation +which does operate on group all, so this may be inconsistent. + +E: Pressure ID for fix modify does not exist + +Self-explanatory. + +E: Could not find fix_modify pressure ID + +The compute ID for computing pressure does not exist. + +E: Fix_modify pressure ID does not compute pressure + +The compute ID assigned to the fix must compute pressure. + +*/ diff --git a/src/RIGID/fix_rigid_nph.cpp b/src/RIGID/fix_rigid_nph.cpp new file mode 100644 index 0000000000..49be257e05 --- /dev/null +++ b/src/RIGID/fix_rigid_nph.cpp @@ -0,0 +1,92 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Tony Sheh (U Michigan), Trung Dac Nguyen (U Michigan) + references: Kamberaj et al., J. Chem. Phys. 122, 224114 (2005) + Miller et al., J Chem Phys. 116, 8649-8659 (2002) +------------------------------------------------------------------------- */ + +#include "string.h" +#include "fix_rigid_nph.h" +#include "domain.h" +#include "modify.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +FixRigidNPH::FixRigidNPH(LAMMPS *lmp, int narg, char **arg) : + FixRigidNH(lmp, narg, arg) +{ + // other setting are made by parent + + scalar_flag = 1; + restart_global = 1; + box_change = 1; + extscalar = 1; + + // error checks + + if (pstat_flag == 0) + error->all(FLERR,"Pressure control must be used with fix nph"); + if (tstat_flag == 1) + error->all(FLERR,"Temperature control must not be used with fix nph"); + if (p_start[0] < 0.0 || p_start[1] < 0.0 || p_start[2] < 0.0 || + p_stop[0] < 0.0 || p_stop[1] < 0.0 || p_stop[2] < 0.0) + error->all(FLERR,"Target pressure for fix rigid/nph cannot be 0.0"); + + // convert input periods to frequency + p_freq[0] = p_freq[1] = p_freq[2] = 0.0; + + if (p_flag[0]) p_freq[0] = 1.0 / p_period[0]; + if (p_flag[1]) p_freq[1] = 1.0 / p_period[1]; + if (p_flag[2]) p_freq[2] = 1.0 / p_period[2]; + + // create a new compute temp style + // id = fix-ID + temp + // compute group = all since pressure is always global (group all) + // and thus its KE/temperature contribution should use group all + + int n = strlen(id) + 6; + id_temp = new char[n]; + strcpy(id_temp,id); + strcat(id_temp,"_temp"); + + char **newarg = new char*[3]; + newarg[0] = id_temp; + newarg[1] = (char *) "all"; + newarg[2] = (char *) "temp"; + modify->add_compute(3,newarg); + delete [] newarg; + tcomputeflag = 1; + + // create a new compute pressure style + // id = fix-ID + press, compute group = all + // pass id_temp as 4th arg to pressure constructor + + n = strlen(id) + 7; + id_press = new char[n]; + strcpy(id_press,id); + strcat(id_press,"_press"); + + newarg = new char*[4]; + newarg[0] = id_press; + newarg[1] = (char *) "all"; + newarg[2] = (char *) "pressure"; + newarg[3] = id_temp; + modify->add_compute(4,newarg); + delete [] newarg; + pcomputeflag = 1; +} diff --git a/src/RIGID/fix_rigid_nph.h b/src/RIGID/fix_rigid_nph.h new file mode 100644 index 0000000000..5184952a29 --- /dev/null +++ b/src/RIGID/fix_rigid_nph.h @@ -0,0 +1,37 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef FIX_CLASS + +FixStyle(rigid/nph,FixRigidNPH) + +#else + +#ifndef LMP_FIX_RIGID_NPH_H +#define LMP_FIX_RIGID_NPH_H + +#include "fix_rigid_nh.h" + +namespace LAMMPS_NS { + +class FixRigidNPH : public FixRigidNH { + public: + FixRigidNPH(class LAMMPS *, int, char **); + ~FixRigidNPH() {} +}; + + +} + +#endif +#endif diff --git a/src/RIGID/fix_rigid_npt.cpp b/src/RIGID/fix_rigid_npt.cpp new file mode 100644 index 0000000000..235faaba51 --- /dev/null +++ b/src/RIGID/fix_rigid_npt.cpp @@ -0,0 +1,104 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Tony Sheh (U Michigan), Trung Dac Nguyen (U Michigan) + references: Kamberaj et al., J. Chem. Phys. 122, 224114 (2005) + Miller et al., J Chem Phys. 116, 8649-8659 (2002) +------------------------------------------------------------------------- */ + +#include "string.h" +#include "fix_rigid_npt.h" +#include "domain.h" +#include "modify.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +FixRigidNPT::FixRigidNPT(LAMMPS *lmp, int narg, char **arg) : + FixRigidNH(lmp, narg, arg) +{ + // other setting are made by parent + + scalar_flag = 1; + restart_global = 1; + box_change = 1; + extscalar = 1; + + // error checks + + if (tstat_flag == 0 || pstat_flag == 0) + error->all(FLERR,"Did not set temp or press for fix rigid/npt"); + if (t_start <= 0.0 || t_stop <= 0.0) + error->all(FLERR,"Target temperature for fix rigid/npt cannot be 0.0"); + if (p_start[0] < 0.0 || p_start[1] < 0.0 || p_start[2] < 0.0 || + p_stop[0] < 0.0 || p_stop[1] < 0.0 || p_stop[2] < 0.0) + error->all(FLERR,"Target pressure for fix rigid/npt cannot be 0.0"); + + if (t_period <= 0.0) error->all(FLERR,"Fix rigid/npt period must be > 0.0"); + + // thermostat chain parameters + + if (t_chain < 1) error->all(FLERR,"Illegal fix_modify command"); + if (t_iter < 1) error->all(FLERR,"Illegal fix_modify command"); + if (t_order != 3 && t_order != 5) + error->all(FLERR,"Fix_modify order must be 3 or 5"); + + // convert input periods to frequency + + t_freq = 0.0; + p_freq[0] = p_freq[1] = p_freq[2] = 0.0; + + t_freq = 1.0 / t_period; + if (p_flag[0]) p_freq[0] = 1.0 / p_period[0]; + if (p_flag[1]) p_freq[1] = 1.0 / p_period[1]; + if (p_flag[2]) p_freq[2] = 1.0 / p_period[2]; + + // create a new compute temp style + // id = fix-ID + temp + // compute group = all since pressure is always global (group all) + // and thus its KE/temperature contribution should use group all + + int n = strlen(id) + 6; + id_temp = new char[n]; + strcpy(id_temp,id); + strcat(id_temp,"_temp"); + + char **newarg = new char*[3]; + newarg[0] = id_temp; + newarg[1] = (char *) "all"; + newarg[2] = (char *) "temp"; + modify->add_compute(3,newarg); + delete [] newarg; + tcomputeflag = 1; + + // create a new compute pressure style + // id = fix-ID + press, compute group = all + // pass id_temp as 4th arg to pressure constructor + + n = strlen(id) + 7; + id_press = new char[n]; + strcpy(id_press,id); + strcat(id_press,"_press"); + + newarg = new char*[4]; + newarg[0] = id_press; + newarg[1] = (char *) "all"; + newarg[2] = (char *) "pressure"; + newarg[3] = id_temp; + modify->add_compute(4,newarg); + delete [] newarg; + pcomputeflag = 1; +} diff --git a/src/RIGID/fix_rigid_npt.h b/src/RIGID/fix_rigid_npt.h new file mode 100644 index 0000000000..044a685437 --- /dev/null +++ b/src/RIGID/fix_rigid_npt.h @@ -0,0 +1,37 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef FIX_CLASS + +FixStyle(rigid/npt,FixRigidNPT) + +#else + +#ifndef LMP_FIX_RIGID_NPT_H +#define LMP_FIX_RIGID_NPT_H + +#include "fix_rigid_nh.h" + +namespace LAMMPS_NS { + +class FixRigidNPT : public FixRigidNH { + public: + FixRigidNPT(class LAMMPS *, int, char **); + ~FixRigidNPT() {} +}; + + +} + +#endif +#endif diff --git a/src/RIGID/fix_rigid_nve.cpp b/src/RIGID/fix_rigid_nve.cpp new file mode 100644 index 0000000000..6cde617e59 --- /dev/null +++ b/src/RIGID/fix_rigid_nve.cpp @@ -0,0 +1,28 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: 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 "fix_rigid_nve.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +FixRigidNVE::FixRigidNVE(LAMMPS *lmp, int narg, char **arg) : + FixRigidNH(lmp, narg, arg) {} + diff --git a/src/RIGID/fix_rigid_nve.h b/src/RIGID/fix_rigid_nve.h new file mode 100644 index 0000000000..2f5d854c1b --- /dev/null +++ b/src/RIGID/fix_rigid_nve.h @@ -0,0 +1,36 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef FIX_CLASS + +FixStyle(rigid/nve,FixRigidNVE) + +#else + +#ifndef LMP_FIX_RIGID_NVE_H +#define LMP_FIX_RIGID_NVE_H + +#include "fix_rigid_nh.h" + +namespace LAMMPS_NS { + +class FixRigidNVE : public FixRigidNH { + public: + FixRigidNVE(class LAMMPS *, int, char **); + ~FixRigidNVE() {} +}; + +} + +#endif +#endif diff --git a/src/RIGID/fix_rigid_nvt.cpp b/src/RIGID/fix_rigid_nvt.cpp new file mode 100644 index 0000000000..abf681d905 --- /dev/null +++ b/src/RIGID/fix_rigid_nvt.cpp @@ -0,0 +1,50 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: 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 "fix_rigid_nvt.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +FixRigidNVT::FixRigidNVT(LAMMPS *lmp, int narg, char **arg) : + FixRigidNH(lmp, narg, arg) +{ + // other settings are made by parent + + scalar_flag = 1; + restart_global = 1; + extscalar = 1; + + // error checking + // convert input period to frequency + + if (tstat_flag == 0) + error->all(FLERR,"Did not set temp for fix rigid/nvt"); + if (t_start < 0.0 || t_stop <= 0.0) + error->all(FLERR,"Target temperature for fix rigid/nvt cannot be 0.0"); + if (t_period <= 0.0) error->all(FLERR,"Fix rigid/nvt period must be > 0.0"); + t_freq = 1.0 / t_period; + + if (t_chain < 1) error->all(FLERR,"Illegal fix_modify command"); + if (t_iter < 1) error->all(FLERR,"Illegal fix_modify command"); + if (t_order != 3 && t_order != 5) + error->all(FLERR,"Fix_modify order must be 3 or 5"); +} diff --git a/src/RIGID/fix_rigid_nvt.h b/src/RIGID/fix_rigid_nvt.h new file mode 100644 index 0000000000..337de6d9ed --- /dev/null +++ b/src/RIGID/fix_rigid_nvt.h @@ -0,0 +1,36 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef FIX_CLASS + +FixStyle(rigid/nvt,FixRigidNVT) + +#else + +#ifndef LMP_FIX_RIGID_NVT_H +#define LMP_FIX_RIGID_NVT_H + +#include "fix_rigid_nh.h" + +namespace LAMMPS_NS { + +class FixRigidNVT : public FixRigidNH { + public: + FixRigidNVT(class LAMMPS *, int, char **); + ~FixRigidNVT() {} +}; + +} + +#endif +#endif diff --git a/src/fix_langevin.cpp b/src/fix_langevin.cpp index 98f827dfdd..4e49b8ede2 100644 --- a/src/fix_langevin.cpp +++ b/src/fix_langevin.cpp @@ -805,7 +805,7 @@ double FixLangevin::compute_scalar() void *FixLangevin::extract(const char *str, int &dim) { - dim=0; + dim = 0; if (strcmp(str,"t_target") == 0) { return &t_target; }