forked from lijiext/lammps
git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@8681 f3b2605a-c512-4ea7-a41b-209d697bcdaa
This commit is contained in:
sjplimp
parent
4cb989a2d3
commit
da4dc73784
|
|
@ -39,6 +39,8 @@ using namespace FixConst;
|
|||
using namespace MathConst;
|
||||
|
||||
enum{SINGLE,MOLECULE,GROUP};
|
||||
enum{NONE,XYZ,XY,YZ,XZ};
|
||||
enum{ISO,ANISO,TRICLINIC};
|
||||
|
||||
#define MAXLINE 256
|
||||
#define CHUNK 1024
|
||||
|
|
@ -241,14 +243,25 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
|
|||
|
||||
int seed;
|
||||
langflag = 0;
|
||||
tempflag = 0;
|
||||
pressflag = 0;
|
||||
tstat_flag = 0;
|
||||
pstat_flag = 0;
|
||||
allremap = 1;
|
||||
id_dilate = NULL;
|
||||
t_chain = 10;
|
||||
t_iter = 1;
|
||||
t_order = 3;
|
||||
p_chain = 10;
|
||||
infile = NULL;
|
||||
|
||||
pcouple = NONE;
|
||||
pstyle = ANISO;
|
||||
dimension = domain->dimension;
|
||||
|
||||
for (int i = 0; i < 3; i++) {
|
||||
p_start[i] = p_stop[i] = p_period[i] = 0.0;
|
||||
p_flag[i] = 0;
|
||||
}
|
||||
|
||||
while (iarg < narg) {
|
||||
if (strcmp(arg[iarg],"force") == 0) {
|
||||
if (iarg+5 > narg) error->all(FLERR,"Illegal fix rigid command");
|
||||
|
|
@ -330,34 +343,104 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
|
|||
if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid command");
|
||||
if (strcmp(style,"rigid/nvt") != 0 && strcmp(style,"rigid/npt") != 0)
|
||||
error->all(FLERR,"Illegal fix rigid command");
|
||||
tempflag = 1;
|
||||
tstat_flag = 1;
|
||||
t_start = atof(arg[iarg+1]);
|
||||
t_stop = atof(arg[iarg+2]);
|
||||
t_period = atof(arg[iarg+3]);
|
||||
iarg += 4;
|
||||
|
||||
} else if (strcmp(arg[iarg],"press") == 0) {
|
||||
} else if (strcmp(arg[iarg],"iso") == 0) {
|
||||
if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid command");
|
||||
if (strcmp(style,"rigid/npt") != 0)
|
||||
if (strcmp(style,"rigid/npt") != 0 && strcmp(style,"rigid/nph") != 0)
|
||||
error->all(FLERR,"Illegal fix rigid command");
|
||||
pressflag = 1;
|
||||
p_start = atof(arg[iarg+1]);
|
||||
p_stop = atof(arg[iarg+2]);
|
||||
p_period = atof(arg[iarg+3]);
|
||||
pcouple = XYZ;
|
||||
p_start[0] = p_start[1] = p_start[2] = atof(arg[iarg+1]);
|
||||
p_stop[0] = p_stop[1] = p_stop[2] = atof(arg[iarg+2]);
|
||||
p_period[0] = p_period[1] = p_period[2] = atof(arg[iarg+3]);
|
||||
p_flag[0] = p_flag[1] = p_flag[2] = 1;
|
||||
if (dimension == 2) {
|
||||
p_start[2] = p_stop[2] = p_period[2] = 0.0;
|
||||
p_flag[2] = 0;
|
||||
}
|
||||
iarg += 4;
|
||||
|
||||
} else if (strcmp(arg[iarg],"aniso") == 0) {
|
||||
if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid command");
|
||||
if (strcmp(style,"rigid/npt") != 0 && strcmp(style,"rigid/nph") != 0)
|
||||
error->all(FLERR,"Illegal fix rigid command");
|
||||
p_start[0] = p_start[1] = p_start[2] = atof(arg[iarg+1]);
|
||||
p_stop[0] = p_stop[1] = p_stop[2] = atof(arg[iarg+2]);
|
||||
p_period[0] = p_period[1] = p_period[2] = atof(arg[iarg+3]);
|
||||
p_flag[0] = p_flag[1] = p_flag[2] = 1;
|
||||
if (dimension == 2) {
|
||||
p_start[2] = p_stop[2] = p_period[2] = 0.0;
|
||||
p_flag[2] = 0;
|
||||
}
|
||||
iarg += 4;
|
||||
|
||||
} else if (strcmp(arg[iarg],"x") == 0) {
|
||||
if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid command");
|
||||
p_start[0] = atof(arg[iarg+1]);
|
||||
p_stop[0] = atof(arg[iarg+2]);
|
||||
p_period[0] = atof(arg[iarg+3]);
|
||||
p_flag[0] = 1;
|
||||
iarg += 4;
|
||||
|
||||
} else if (strcmp(arg[iarg],"y") == 0) {
|
||||
if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid command");
|
||||
p_start[1] = atof(arg[iarg+1]);
|
||||
p_stop[1] = atof(arg[iarg+2]);
|
||||
p_period[1] = atof(arg[iarg+3]);
|
||||
p_flag[1] = 1;
|
||||
iarg += 4;
|
||||
|
||||
} else if (strcmp(arg[iarg],"z") == 0) {
|
||||
if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid command");
|
||||
p_start[2] = atof(arg[iarg+1]);
|
||||
p_stop[2] = atof(arg[iarg+2]);
|
||||
p_period[2] = atof(arg[iarg+3]);
|
||||
p_flag[2] = 1;
|
||||
iarg += 4;
|
||||
|
||||
} else if (strcmp(arg[iarg],"couple") == 0) {
|
||||
if (iarg+2 > narg) error->all(FLERR,"Illegal fix rigid command");
|
||||
if (strcmp(arg[iarg+1],"xyz") == 0) pcouple = XYZ;
|
||||
else if (strcmp(arg[iarg+1],"xy") == 0) pcouple = XY;
|
||||
else if (strcmp(arg[iarg+1],"yz") == 0) pcouple = YZ;
|
||||
else if (strcmp(arg[iarg+1],"xz") == 0) pcouple = XZ;
|
||||
else if (strcmp(arg[iarg+1],"none") == 0) pcouple = NONE;
|
||||
else error->all(FLERR,"Illegal fix rigid command");
|
||||
iarg += 2;
|
||||
|
||||
} else if (strcmp(arg[iarg],"dilate") == 0) {
|
||||
if (iarg+2 > narg)
|
||||
error->all(FLERR,"Illegal fix rigid nvt/npt/nph command");
|
||||
if (strcmp(arg[iarg+1],"all") == 0) allremap = 1;
|
||||
else {
|
||||
allremap = 0;
|
||||
delete [] id_dilate;
|
||||
int n = strlen(arg[iarg+1]) + 1;
|
||||
id_dilate = new char[n];
|
||||
strcpy(id_dilate,arg[iarg+1]);
|
||||
int idilate = group->find(id_dilate);
|
||||
if (idilate == -1)
|
||||
error->all(FLERR,
|
||||
"Fix rigid nvt/npt/nph dilate group ID does not exist");
|
||||
}
|
||||
iarg += 2;
|
||||
|
||||
} else if (strcmp(arg[iarg],"tparam") == 0) {
|
||||
if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid command");
|
||||
if (strcmp(style,"rigid/nvt") != 0)
|
||||
if (strcmp(style,"rigid/nvt") != 0 && strcmp(style,"rigid/npt") != 0)
|
||||
error->all(FLERR,"Illegal fix rigid command");
|
||||
t_chain = atoi(arg[iarg+1]);
|
||||
t_iter = atoi(arg[iarg+2]);
|
||||
t_order = atoi(arg[iarg+3]);
|
||||
iarg += 4;
|
||||
|
||||
} else if (strcmp(arg[iarg],"pparam") == 0) {
|
||||
} else if (strcmp(arg[iarg],"pchain") == 0) {
|
||||
if (iarg+2 > narg) error->all(FLERR,"Illegal fix rigid command");
|
||||
if (strcmp(style,"rigid/npt") != 0)
|
||||
if (strcmp(style,"rigid/npt") != 0 && strcmp(style,"rigid/nph") != 0)
|
||||
error->all(FLERR,"Illegal fix rigid command");
|
||||
p_chain = atoi(arg[iarg+1]);
|
||||
iarg += 2;
|
||||
|
|
@ -372,7 +455,15 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
|
|||
|
||||
} else error->all(FLERR,"Illegal fix rigid command");
|
||||
}
|
||||
t_target = t_start;
|
||||
|
||||
// set pstat_flag
|
||||
|
||||
pstat_flag = 0;
|
||||
for (int i = 0; i < 3; i++)
|
||||
if (p_flag[i]) pstat_flag = 1;
|
||||
|
||||
if (pcouple == XYZ || (dimension == 2 && pcouple == XY)) pstyle = ISO;
|
||||
else pstyle = ANISO;
|
||||
|
||||
// initialize Marsaglia RNG with processor-unique seed
|
||||
|
||||
|
|
|
|||
|
|
@ -41,6 +41,7 @@ class FixRigid : public Fix {
|
|||
void initial_integrate_respa(int, int, int);
|
||||
void final_integrate_respa(int, int);
|
||||
virtual double compute_scalar();
|
||||
virtual int modify_param(int, char **) {return 0;}
|
||||
|
||||
double memory_usage();
|
||||
void grow_arrays(int);
|
||||
|
|
@ -67,6 +68,7 @@ class FixRigid : public Fix {
|
|||
int firstflag; // 1 for first-time setup of rigid bodies
|
||||
char *infile; // file to read rigid body attributes from
|
||||
|
||||
int dimension; // # of dimensions
|
||||
int nbody; // # of rigid bodies
|
||||
int *nrigid; // # of atoms in each rigid body
|
||||
int *mol2body; // convert mol-ID to rigid body index
|
||||
|
|
@ -103,16 +105,22 @@ class FixRigid : public Fix {
|
|||
double tfactor; // scale factor on temperature of rigid bodies
|
||||
int langflag; // 0/1 = no/yes Langevin thermostat
|
||||
|
||||
int tempflag; // NVT settings
|
||||
int tstat_flag; // NVT settings
|
||||
double t_start,t_stop,t_target;
|
||||
double t_period,t_freq;
|
||||
int t_chain,t_iter,t_order;
|
||||
|
||||
int pressflag; // NPT settings
|
||||
double p_start,p_stop;
|
||||
double p_period,p_freq;
|
||||
int pstat_flag; // NPT settings
|
||||
double p_start[3],p_stop[3];
|
||||
double p_period[3],p_freq[3];
|
||||
int p_flag[3];
|
||||
int pcouple,pstyle;
|
||||
int p_chain;
|
||||
|
||||
int allremap; // remap all atoms
|
||||
int dilate_group_bit; // mask for dilation group
|
||||
char *id_dilate; // group name to dilate
|
||||
|
||||
class RanMars *random;
|
||||
class AtomVecEllipsoid *avec_ellipsoid;
|
||||
class AtomVecLine *avec_line;
|
||||
|
|
|
|||
File diff suppressed because it is too large
Load Diff
|
|
@ -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.
|
||||
|
||||
*/
|
||||
|
|
@ -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;
|
||||
}
|
||||
|
|
@ -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
|
||||
|
|
@ -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;
|
||||
}
|
||||
|
|
@ -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
|
||||
|
|
@ -14,268 +14,15 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
Contributing author: Tony Sheh (U Michigan), Trung Dac Nguyen (U Michigan)
|
||||
references: Kamberaj et al., J. Chem. Phys. 122, 224114 (2005)
|
||||
Miller et al., J Chem Phys. 116, 8649-8659
|
||||
Miller et al., J Chem Phys. 116, 8649-8659 (2002)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include "math.h"
|
||||
#include "stdio.h"
|
||||
#include "string.h"
|
||||
#include "fix_rigid_nve.h"
|
||||
#include "math_extra.h"
|
||||
#include "atom.h"
|
||||
#include "domain.h"
|
||||
#include "update.h"
|
||||
#include "modify.h"
|
||||
#include "group.h"
|
||||
#include "comm.h"
|
||||
#include "force.h"
|
||||
#include "output.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace FixConst;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixRigidNVE::FixRigidNVE(LAMMPS *lmp, int narg, char **arg) :
|
||||
FixRigid(lmp, narg, arg)
|
||||
{
|
||||
memory->create(conjqm,nbody,4,"rigid/nve:conjqm");
|
||||
}
|
||||
FixRigidNH(lmp, narg, arg) {}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixRigidNVE::~FixRigidNVE()
|
||||
{
|
||||
memory->destroy(conjqm);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixRigidNVE::setup(int vflag)
|
||||
{
|
||||
FixRigid::setup(vflag);
|
||||
|
||||
double mbody[3];
|
||||
for (int ibody = 0; ibody < nbody; ibody++) {
|
||||
MathExtra::transpose_matvec(ex_space[ibody],ey_space[ibody],ez_space[ibody],
|
||||
angmom[ibody],mbody);
|
||||
MathExtra::quatvec(quat[ibody],mbody,conjqm[ibody]);
|
||||
conjqm[ibody][0] *= 2.0;
|
||||
conjqm[ibody][1] *= 2.0;
|
||||
conjqm[ibody][2] *= 2.0;
|
||||
conjqm[ibody][3] *= 2.0;
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
perform preforce velocity Verlet integration
|
||||
see Kamberaj paper for step references
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void FixRigidNVE::initial_integrate(int vflag)
|
||||
{
|
||||
double dtfm,mbody[3],tbody[3],fquat[4];
|
||||
double dtf2 = dtf * 2.0;
|
||||
|
||||
for (int ibody = 0; ibody < nbody; ibody++) {
|
||||
|
||||
// step 1.1 - update vcm by 1/2 step
|
||||
|
||||
dtfm = dtf / masstotal[ibody];
|
||||
vcm[ibody][0] += dtfm * fcm[ibody][0] * fflag[ibody][0];
|
||||
vcm[ibody][1] += dtfm * fcm[ibody][1] * fflag[ibody][1];
|
||||
vcm[ibody][2] += dtfm * fcm[ibody][2] * fflag[ibody][2];
|
||||
|
||||
// step 1.2 - update xcm by full step
|
||||
|
||||
xcm[ibody][0] += dtv * vcm[ibody][0];
|
||||
xcm[ibody][1] += dtv * vcm[ibody][1];
|
||||
xcm[ibody][2] += dtv * vcm[ibody][2];
|
||||
|
||||
// step 1.3 - apply torque (body coords) to quaternion momentum
|
||||
|
||||
torque[ibody][0] *= tflag[ibody][0];
|
||||
torque[ibody][1] *= tflag[ibody][1];
|
||||
torque[ibody][2] *= tflag[ibody][2];
|
||||
|
||||
MathExtra::transpose_matvec(ex_space[ibody],ey_space[ibody],ez_space[ibody],
|
||||
torque[ibody],tbody);
|
||||
MathExtra::quatvec(quat[ibody],tbody,fquat);
|
||||
|
||||
conjqm[ibody][0] += dtf2 * fquat[0];
|
||||
conjqm[ibody][1] += dtf2 * fquat[1];
|
||||
conjqm[ibody][2] += dtf2 * fquat[2];
|
||||
conjqm[ibody][3] += dtf2 * fquat[3];
|
||||
|
||||
// step 1.4 to 1.13 - use no_squish rotate to update p and q
|
||||
|
||||
no_squish_rotate(3,conjqm[ibody],quat[ibody],inertia[ibody],dtq);
|
||||
no_squish_rotate(2,conjqm[ibody],quat[ibody],inertia[ibody],dtq);
|
||||
no_squish_rotate(1,conjqm[ibody],quat[ibody],inertia[ibody],dtv);
|
||||
no_squish_rotate(2,conjqm[ibody],quat[ibody],inertia[ibody],dtq);
|
||||
no_squish_rotate(3,conjqm[ibody],quat[ibody],inertia[ibody],dtq);
|
||||
|
||||
// update exyz_space
|
||||
// transform p back to angmom
|
||||
// update angular velocity
|
||||
|
||||
MathExtra::q_to_exyz(quat[ibody],ex_space[ibody],ey_space[ibody],
|
||||
ez_space[ibody]);
|
||||
MathExtra::invquatvec(quat[ibody],conjqm[ibody],mbody);
|
||||
MathExtra::matvec(ex_space[ibody],ey_space[ibody],ez_space[ibody],
|
||||
mbody,angmom[ibody]);
|
||||
|
||||
angmom[ibody][0] *= 0.5;
|
||||
angmom[ibody][1] *= 0.5;
|
||||
angmom[ibody][2] *= 0.5;
|
||||
|
||||
MathExtra::angmom_to_omega(angmom[ibody],ex_space[ibody],ey_space[ibody],
|
||||
ez_space[ibody],inertia[ibody],omega[ibody]);
|
||||
}
|
||||
|
||||
// virial setup before call to set_xv
|
||||
|
||||
if (vflag) v_setup(vflag);
|
||||
else evflag = 0;
|
||||
|
||||
// set coords/orient and velocity/rotation of atoms in rigid bodies
|
||||
// from quarternion and omega
|
||||
|
||||
set_xv();
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixRigidNVE::final_integrate()
|
||||
{
|
||||
int i,ibody;
|
||||
double dtfm,xy,xz,yz;
|
||||
|
||||
// sum over atoms to get force and torque on rigid body
|
||||
|
||||
tagint *image = atom->image;
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
double xprd = domain->xprd;
|
||||
double yprd = domain->yprd;
|
||||
double zprd = domain->zprd;
|
||||
if (triclinic) {
|
||||
xy = domain->xy;
|
||||
xz = domain->xz;
|
||||
yz = domain->yz;
|
||||
}
|
||||
|
||||
int xbox,ybox,zbox;
|
||||
double xunwrap,yunwrap,zunwrap,dx,dy,dz;
|
||||
for (ibody = 0; ibody < nbody; ibody++)
|
||||
for (i = 0; i < 6; i++) sum[ibody][i] = 0.0;
|
||||
|
||||
for (i = 0; i < nlocal; i++) {
|
||||
if (body[i] < 0) continue;
|
||||
ibody = body[i];
|
||||
|
||||
sum[ibody][0] += f[i][0];
|
||||
sum[ibody][1] += f[i][1];
|
||||
sum[ibody][2] += f[i][2];
|
||||
|
||||
xbox = (image[i] & IMGMASK) - IMGMAX;
|
||||
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
|
||||
zbox = (image[i] >> IMG2BITS) - IMGMAX;
|
||||
|
||||
if (triclinic == 0) {
|
||||
xunwrap = x[i][0] + xbox*xprd;
|
||||
yunwrap = x[i][1] + ybox*yprd;
|
||||
zunwrap = x[i][2] + zbox*zprd;
|
||||
} else {
|
||||
xunwrap = x[i][0] + xbox*xprd + ybox*xy + zbox*xz;
|
||||
yunwrap = x[i][1] + ybox*yprd + zbox*yz;
|
||||
zunwrap = x[i][2] + zbox*zprd;
|
||||
}
|
||||
|
||||
dx = xunwrap - xcm[ibody][0];
|
||||
dy = yunwrap - xcm[ibody][1];
|
||||
dz = zunwrap - xcm[ibody][2];
|
||||
|
||||
sum[ibody][3] += dy*f[i][2] - dz*f[i][1];
|
||||
sum[ibody][4] += dz*f[i][0] - dx*f[i][2];
|
||||
sum[ibody][5] += dx*f[i][1] - dy*f[i][0];
|
||||
}
|
||||
|
||||
// extended particles add their torque to torque of body
|
||||
|
||||
if (extended) {
|
||||
double **torque_one = atom->torque;
|
||||
|
||||
for (i = 0; i < nlocal; i++) {
|
||||
if (body[i] < 0) continue;
|
||||
ibody = body[i];
|
||||
|
||||
if (eflags[i] & TORQUE) {
|
||||
sum[ibody][3] += torque_one[i][0];
|
||||
sum[ibody][4] += torque_one[i][1];
|
||||
sum[ibody][5] += torque_one[i][2];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world);
|
||||
|
||||
// update vcm and angmom
|
||||
// include Langevin thermostat forces
|
||||
// fflag,tflag = 0 for some dimensions in 2d
|
||||
|
||||
double mbody[3],tbody[3],fquat[4];
|
||||
double dtf2 = dtf * 2.0;
|
||||
|
||||
for (ibody = 0; ibody < nbody; ibody++) {
|
||||
fcm[ibody][0] = all[ibody][0] + langextra[ibody][0];
|
||||
fcm[ibody][1] = all[ibody][1] + langextra[ibody][1];
|
||||
fcm[ibody][2] = all[ibody][2] + langextra[ibody][2];
|
||||
torque[ibody][0] = all[ibody][3] + langextra[ibody][3];
|
||||
torque[ibody][1] = all[ibody][4] + langextra[ibody][4];
|
||||
torque[ibody][2] = all[ibody][5] + langextra[ibody][5];
|
||||
|
||||
// update vcm by 1/2 step
|
||||
|
||||
dtfm = dtf / masstotal[ibody];
|
||||
vcm[ibody][0] += dtfm * fcm[ibody][0] * fflag[ibody][0];
|
||||
vcm[ibody][1] += dtfm * fcm[ibody][1] * fflag[ibody][1];
|
||||
vcm[ibody][2] += dtfm * fcm[ibody][2] * fflag[ibody][2];
|
||||
|
||||
// update conjqm, then transform to angmom, set velocity again
|
||||
// virial is already setup from initial_integrate
|
||||
|
||||
torque[ibody][0] *= tflag[ibody][0];
|
||||
torque[ibody][1] *= tflag[ibody][1];
|
||||
torque[ibody][2] *= tflag[ibody][2];
|
||||
|
||||
MathExtra::transpose_matvec(ex_space[ibody],ey_space[ibody],ez_space[ibody],
|
||||
torque[ibody],tbody);
|
||||
MathExtra::quatvec(quat[ibody],tbody,fquat);
|
||||
|
||||
conjqm[ibody][0] += dtf2 * fquat[0];
|
||||
conjqm[ibody][1] += dtf2 * fquat[1];
|
||||
conjqm[ibody][2] += dtf2 * fquat[2];
|
||||
conjqm[ibody][3] += dtf2 * fquat[3];
|
||||
|
||||
MathExtra::invquatvec(quat[ibody],conjqm[ibody],mbody);
|
||||
MathExtra::matvec(ex_space[ibody],ey_space[ibody],ez_space[ibody],
|
||||
mbody,angmom[ibody]);
|
||||
|
||||
angmom[ibody][0] *= 0.5;
|
||||
angmom[ibody][1] *= 0.5;
|
||||
angmom[ibody][2] *= 0.5;
|
||||
|
||||
MathExtra::angmom_to_omega(angmom[ibody],ex_space[ibody],ey_space[ibody],
|
||||
ez_space[ibody],inertia[ibody],omega[ibody]);
|
||||
}
|
||||
|
||||
// set velocity/rotation of atoms in rigid bodies
|
||||
// virial is already setup from initial_integrate
|
||||
|
||||
set_v();
|
||||
}
|
||||
|
|
|
|||
|
|
@ -20,26 +20,17 @@ FixStyle(rigid/nve,FixRigidNVE)
|
|||
#ifndef LMP_FIX_RIGID_NVE_H
|
||||
#define LMP_FIX_RIGID_NVE_H
|
||||
|
||||
#include "fix_rigid.h"
|
||||
#include "fix_rigid_nh.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class FixRigidNVE : public FixRigid {
|
||||
class FixRigidNVE : public FixRigidNH {
|
||||
public:
|
||||
FixRigidNVE(class LAMMPS *, int, char **);
|
||||
~FixRigidNVE();
|
||||
void setup(int);
|
||||
void initial_integrate(int);
|
||||
void final_integrate();
|
||||
|
||||
private:
|
||||
double **conjqm; // conjugate quaternion momentum
|
||||
~FixRigidNVE() {}
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
||||
/* ERROR/WARNING messages:
|
||||
|
||||
*/
|
||||
|
|
|
|||
|
|
@ -14,35 +14,20 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
Contributing author: Tony Sheh (U Michigan), Trung Dac Nguyen (U Michigan)
|
||||
references: Kamberaj et al., J. Chem. Phys. 122, 224114 (2005)
|
||||
Miller et al., J Chem Phys. 116, 8649-8659
|
||||
Miller et al., J Chem Phys. 116, 8649-8659 (2002)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include "math.h"
|
||||
#include "stdio.h"
|
||||
#include "stdlib.h"
|
||||
#include "string.h"
|
||||
#include "fix_rigid_nvt.h"
|
||||
#include "math_extra.h"
|
||||
#include "atom.h"
|
||||
#include "domain.h"
|
||||
#include "update.h"
|
||||
#include "modify.h"
|
||||
#include "group.h"
|
||||
#include "comm.h"
|
||||
#include "force.h"
|
||||
#include "output.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace FixConst;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixRigidNVT::FixRigidNVT(LAMMPS *lmp, int narg, char **arg) :
|
||||
FixRigid(lmp, narg, arg)
|
||||
FixRigidNH(lmp, narg, arg)
|
||||
{
|
||||
// other settings are made by FixRigid parent
|
||||
// other settings are made by parent
|
||||
|
||||
scalar_flag = 1;
|
||||
restart_global = 1;
|
||||
|
|
@ -51,7 +36,7 @@ FixRigidNVT::FixRigidNVT(LAMMPS *lmp, int narg, char **arg) :
|
|||
// error checking
|
||||
// convert input period to frequency
|
||||
|
||||
if (tempflag == 0)
|
||||
if (tstat_flag == 0)
|
||||
error->all(FLERR,"Did not set temp for fix rigid/nvt");
|
||||
if (t_start < 0.0 || t_stop <= 0.0)
|
||||
error->all(FLERR,"Target temperature for fix rigid/nvt cannot be 0.0");
|
||||
|
|
@ -62,655 +47,4 @@ FixRigidNVT::FixRigidNVT(LAMMPS *lmp, int narg, char **arg) :
|
|||
if (t_iter < 1) error->all(FLERR,"Illegal fix_modify command");
|
||||
if (t_order != 3 && t_order != 5)
|
||||
error->all(FLERR,"Fix_modify order must be 3 or 5");
|
||||
|
||||
allocate_chain();
|
||||
allocate_order();
|
||||
memory->create(conjqm,nbody,4,"nve_rigid:conjqm");
|
||||
|
||||
// one-time initialize of thermostat variables
|
||||
|
||||
eta_t[0] = eta_r[0] = 0.0;
|
||||
eta_dot_t[0] = eta_dot_r[0] = 0.0;
|
||||
f_eta_t[0] = f_eta_r[0] = 0.0;
|
||||
|
||||
for (int i = 1; i < t_chain; i++) {
|
||||
eta_t[i] = eta_r[i] = 0.0;
|
||||
eta_dot_t[i] = eta_dot_r[i] = 0.0;
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixRigidNVT::~FixRigidNVT()
|
||||
{
|
||||
deallocate_chain();
|
||||
deallocate_order();
|
||||
memory->destroy(conjqm);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int FixRigidNVT::setmask()
|
||||
{
|
||||
int mask = 0;
|
||||
mask |= INITIAL_INTEGRATE;
|
||||
mask |= FINAL_INTEGRATE;
|
||||
mask |= PRE_NEIGHBOR;
|
||||
mask |= THERMO_ENERGY;
|
||||
mask |= INITIAL_INTEGRATE_RESPA;
|
||||
mask |= FINAL_INTEGRATE_RESPA;
|
||||
return mask;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixRigidNVT::init()
|
||||
{
|
||||
FixRigid::init();
|
||||
|
||||
// initialize thermostats
|
||||
// set timesteps, constants
|
||||
// store Yoshida-Suzuki integrator parameters
|
||||
|
||||
dtv = update->dt;
|
||||
dtf = 0.5 * update->dt * force->ftm2v;
|
||||
dtq = 0.5 * update->dt;
|
||||
|
||||
boltz = force->boltz;
|
||||
nf_t = nf_r = domain->dimension * nbody;
|
||||
|
||||
for (int ibody = 0; ibody < nbody; ibody++)
|
||||
for (int k = 0; k < domain->dimension; k++)
|
||||
if (fabs(inertia[ibody][k]) < 1e-6) nf_r--;
|
||||
|
||||
// see Table 1 in Kamberaj et al
|
||||
|
||||
if (t_order == 3) {
|
||||
w[0] = 1.0 / (2.0 - pow(2.0, 1.0/3.0));
|
||||
w[1] = 1.0 - 2.0*w[0];
|
||||
w[2] = w[0];
|
||||
} else if (t_order == 5) {
|
||||
w[0] = 1.0 / (4.0 - pow(4.0, 1.0/3.0));
|
||||
w[1] = w[0];
|
||||
w[2] = 1.0 - 4.0 * w[0];
|
||||
w[3] = w[0];
|
||||
w[4] = w[0];
|
||||
}
|
||||
|
||||
// initialize thermostat settings
|
||||
|
||||
t_target = t_start;
|
||||
double kt = boltz * t_target;
|
||||
double t_mass = kt / (t_freq*t_freq);
|
||||
q_t[0] = nf_t * t_mass;
|
||||
q_r[0] = nf_r * t_mass;
|
||||
for (int i = 1; i < t_chain; i++)
|
||||
q_t[i] = q_r[i] = t_mass;
|
||||
|
||||
// initialize thermostat chain positions, velocites, forces
|
||||
|
||||
for (int i = 1; i < t_chain; i++) {
|
||||
f_eta_t[i] = q_t[i-1] * eta_dot_t[i-1] * eta_dot_t[i-1] - kt;
|
||||
f_eta_r[i] = q_r[i-1] * eta_dot_r[i-1] * eta_dot_r[i-1] - kt;
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixRigidNVT::setup(int vflag)
|
||||
{
|
||||
FixRigid::setup(vflag);
|
||||
|
||||
t_target = t_start;
|
||||
|
||||
double mbody[3];
|
||||
for (int ibody = 0; ibody < nbody; ibody++) {
|
||||
MathExtra::transpose_matvec(ex_space[ibody],ey_space[ibody],ez_space[ibody],
|
||||
angmom[ibody],mbody);
|
||||
MathExtra::quatvec(quat[ibody],mbody,conjqm[ibody]);
|
||||
conjqm[ibody][0] *= 2.0;
|
||||
conjqm[ibody][1] *= 2.0;
|
||||
conjqm[ibody][2] *= 2.0;
|
||||
conjqm[ibody][3] *= 2.0;
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
perform preforce velocity Verlet integration
|
||||
see Kamberaj paper for step references
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void FixRigidNVT::initial_integrate(int vflag)
|
||||
{
|
||||
double tmp,akin_t,akin_r,scale_t,scale_r;
|
||||
double dtfm,mbody[3],tbody[3],fquat[4];
|
||||
double dtf2 = dtf * 2.0;
|
||||
|
||||
double delta = update->ntimestep - update->beginstep;
|
||||
delta /= update->endstep - update->beginstep;
|
||||
t_target = t_start + delta * (t_stop - t_start);
|
||||
|
||||
// intialize velocity scale for translation and rotation
|
||||
|
||||
akin_t = akin_r = 0.0;
|
||||
tmp = -1.0 * dtq * eta_dot_t[0];
|
||||
scale_t = exp(tmp);
|
||||
tmp = -1.0 * dtq * eta_dot_r[0];
|
||||
scale_r = exp(tmp);
|
||||
|
||||
for (int ibody = 0; ibody < nbody; ibody++) {
|
||||
|
||||
// step 1.1 - update vcm by 1/2 step
|
||||
|
||||
dtfm = dtf / masstotal[ibody];
|
||||
vcm[ibody][0] += dtfm * fcm[ibody][0] * fflag[ibody][0];
|
||||
vcm[ibody][1] += dtfm * fcm[ibody][1] * fflag[ibody][1];
|
||||
vcm[ibody][2] += dtfm * fcm[ibody][2] * fflag[ibody][2];
|
||||
vcm[ibody][0] *= scale_t;
|
||||
vcm[ibody][1] *= scale_t;
|
||||
vcm[ibody][2] *= scale_t;
|
||||
|
||||
tmp = vcm[ibody][0]*vcm[ibody][0] + vcm[ibody][1]*vcm[ibody][1] +
|
||||
vcm[ibody][2]*vcm[ibody][2];
|
||||
akin_t += masstotal[ibody]*tmp;
|
||||
|
||||
// step 1.2 - update xcm by full step
|
||||
|
||||
xcm[ibody][0] += dtv * vcm[ibody][0];
|
||||
xcm[ibody][1] += dtv * vcm[ibody][1];
|
||||
xcm[ibody][2] += dtv * vcm[ibody][2];
|
||||
|
||||
torque[ibody][0] *= tflag[ibody][0];
|
||||
torque[ibody][1] *= tflag[ibody][1];
|
||||
torque[ibody][2] *= tflag[ibody][2];
|
||||
|
||||
// step 1.3 - apply torque (body coords) to quaternion momentum
|
||||
|
||||
MathExtra::transpose_matvec(ex_space[ibody],ey_space[ibody],ez_space[ibody],
|
||||
torque[ibody],tbody);
|
||||
MathExtra::quatvec(quat[ibody],tbody,fquat);
|
||||
|
||||
conjqm[ibody][0] += dtf2 * fquat[0];
|
||||
conjqm[ibody][1] += dtf2 * fquat[1];
|
||||
conjqm[ibody][2] += dtf2 * fquat[2];
|
||||
conjqm[ibody][3] += dtf2 * fquat[3];
|
||||
conjqm[ibody][0] *= scale_r;
|
||||
conjqm[ibody][1] *= scale_r;
|
||||
conjqm[ibody][2] *= scale_r;
|
||||
conjqm[ibody][3] *= scale_r;
|
||||
|
||||
// step 1.4 to 1.8 - use no_squish rotate to update p (i.e. conjqm) and q
|
||||
|
||||
no_squish_rotate(3,conjqm[ibody],quat[ibody],inertia[ibody],dtq);
|
||||
no_squish_rotate(2,conjqm[ibody],quat[ibody],inertia[ibody],dtq);
|
||||
no_squish_rotate(1,conjqm[ibody],quat[ibody],inertia[ibody],dtv);
|
||||
no_squish_rotate(2,conjqm[ibody],quat[ibody],inertia[ibody],dtq);
|
||||
no_squish_rotate(3,conjqm[ibody],quat[ibody],inertia[ibody],dtq);
|
||||
|
||||
// update the exyz_space from new quaternion
|
||||
// transform p back to angmom
|
||||
// update angular velocity
|
||||
|
||||
MathExtra::q_to_exyz(quat[ibody],ex_space[ibody],ey_space[ibody],
|
||||
ez_space[ibody]);
|
||||
MathExtra::invquatvec(quat[ibody],conjqm[ibody],mbody);
|
||||
MathExtra::matvec(ex_space[ibody],ey_space[ibody],ez_space[ibody],
|
||||
mbody,angmom[ibody]);
|
||||
|
||||
angmom[ibody][0] *= 0.5;
|
||||
angmom[ibody][1] *= 0.5;
|
||||
angmom[ibody][2] *= 0.5;
|
||||
|
||||
MathExtra::angmom_to_omega(angmom[ibody],ex_space[ibody],ey_space[ibody],
|
||||
ez_space[ibody],inertia[ibody],omega[ibody]);
|
||||
|
||||
akin_r += angmom[ibody][0]*omega[ibody][0] +
|
||||
angmom[ibody][1]*omega[ibody][1] + angmom[ibody][2]*omega[ibody][2];
|
||||
}
|
||||
|
||||
// update thermostat chains
|
||||
|
||||
update_nhcp(akin_t,akin_r);
|
||||
|
||||
// virial setup before call to set_xv
|
||||
|
||||
if (vflag) v_setup(vflag);
|
||||
else evflag = 0;
|
||||
|
||||
// set coords/orient and velocity/rotation of atoms in rigid bodies
|
||||
// from quarternion and omega
|
||||
|
||||
set_xv();
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixRigidNVT::final_integrate()
|
||||
{
|
||||
int i,ibody;
|
||||
double tmp,scale_t,scale_r;
|
||||
double dtfm,xy,xz,yz;
|
||||
|
||||
// compute velocity scales for translation and rotation
|
||||
|
||||
tmp = -1.0 * dtq * eta_dot_t[0];
|
||||
scale_t = exp(tmp);
|
||||
tmp = -1.0 * dtq * eta_dot_r[0];
|
||||
scale_r = exp(tmp);
|
||||
|
||||
// sum over atoms to get force and torque on rigid body
|
||||
|
||||
tagint *image = atom->image;
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
double xprd = domain->xprd;
|
||||
double yprd = domain->yprd;
|
||||
double zprd = domain->zprd;
|
||||
if (triclinic) {
|
||||
xy = domain->xy;
|
||||
xz = domain->xz;
|
||||
yz = domain->yz;
|
||||
}
|
||||
|
||||
int xbox,ybox,zbox;
|
||||
double xunwrap,yunwrap,zunwrap,dx,dy,dz;
|
||||
for (ibody = 0; ibody < nbody; ibody++)
|
||||
for (i = 0; i < 6; i++) sum[ibody][i] = 0.0;
|
||||
|
||||
for (i = 0; i < nlocal; i++) {
|
||||
if (body[i] < 0) continue;
|
||||
ibody = body[i];
|
||||
|
||||
sum[ibody][0] += f[i][0];
|
||||
sum[ibody][1] += f[i][1];
|
||||
sum[ibody][2] += f[i][2];
|
||||
|
||||
xbox = (image[i] & IMGMASK) - IMGMAX;
|
||||
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
|
||||
zbox = (image[i] >> IMG2BITS) - IMGMAX;
|
||||
|
||||
if (triclinic == 0) {
|
||||
xunwrap = x[i][0] + xbox*xprd;
|
||||
yunwrap = x[i][1] + ybox*yprd;
|
||||
zunwrap = x[i][2] + zbox*zprd;
|
||||
} else {
|
||||
xunwrap = x[i][0] + xbox*xprd + ybox*xy + zbox*xz;
|
||||
yunwrap = x[i][1] + ybox*yprd + zbox*yz;
|
||||
zunwrap = x[i][2] + zbox*zprd;
|
||||
}
|
||||
|
||||
dx = xunwrap - xcm[ibody][0];
|
||||
dy = yunwrap - xcm[ibody][1];
|
||||
dz = zunwrap - xcm[ibody][2];
|
||||
|
||||
sum[ibody][3] += dy*f[i][2] - dz*f[i][1];
|
||||
sum[ibody][4] += dz*f[i][0] - dx*f[i][2];
|
||||
sum[ibody][5] += dx*f[i][1] - dy*f[i][0];
|
||||
}
|
||||
|
||||
// extended particles add their torque to torque of body
|
||||
|
||||
if (extended) {
|
||||
double **torque_one = atom->torque;
|
||||
|
||||
for (i = 0; i < nlocal; i++) {
|
||||
if (body[i] < 0) continue;
|
||||
ibody = body[i];
|
||||
|
||||
if (eflags[i] & TORQUE) {
|
||||
sum[ibody][3] += torque_one[i][0];
|
||||
sum[ibody][4] += torque_one[i][1];
|
||||
sum[ibody][5] += torque_one[i][2];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world);
|
||||
|
||||
double mbody[3],tbody[3],fquat[4];
|
||||
double dtf2 = dtf * 2.0;
|
||||
|
||||
for (ibody = 0; ibody < nbody; ibody++) {
|
||||
fcm[ibody][0] = all[ibody][0];
|
||||
fcm[ibody][1] = all[ibody][1];
|
||||
fcm[ibody][2] = all[ibody][2];
|
||||
torque[ibody][0] = all[ibody][3];
|
||||
torque[ibody][1] = all[ibody][4];
|
||||
torque[ibody][2] = all[ibody][5];
|
||||
|
||||
// 2.5-2.6 update vcm by 1/2 step
|
||||
|
||||
dtfm = dtf / masstotal[ibody];
|
||||
vcm[ibody][0] *= scale_t;
|
||||
vcm[ibody][1] *= scale_t;
|
||||
vcm[ibody][2] *= scale_t;
|
||||
vcm[ibody][0] += dtfm * fcm[ibody][0] * fflag[ibody][0];
|
||||
vcm[ibody][1] += dtfm * fcm[ibody][1] * fflag[ibody][1];
|
||||
vcm[ibody][2] += dtfm * fcm[ibody][2] * fflag[ibody][2];
|
||||
|
||||
// 2.1-2.4 update conjqm, angular momentum and angular velocity
|
||||
// apply body torque flags
|
||||
|
||||
torque[ibody][0] *= tflag[ibody][0];
|
||||
torque[ibody][1] *= tflag[ibody][1];
|
||||
torque[ibody][2] *= tflag[ibody][2];
|
||||
|
||||
// convert torque to the body frame
|
||||
|
||||
MathExtra::transpose_matvec(ex_space[ibody],ey_space[ibody],ez_space[ibody],
|
||||
torque[ibody],tbody);
|
||||
|
||||
// compute "force" for quaternion
|
||||
|
||||
MathExtra::quatvec(quat[ibody],tbody,fquat);
|
||||
|
||||
// update the conjugate quaternion momentum (conjqm)
|
||||
|
||||
conjqm[ibody][0] = scale_r * conjqm[ibody][0] + dtf2 * fquat[0];
|
||||
conjqm[ibody][1] = scale_r * conjqm[ibody][1] + dtf2 * fquat[1];
|
||||
conjqm[ibody][2] = scale_r * conjqm[ibody][2] + dtf2 * fquat[2];
|
||||
conjqm[ibody][3] = scale_r * conjqm[ibody][3] + dtf2 * fquat[3];
|
||||
|
||||
// compute angular momentum in the body frame
|
||||
// then convert to the space-fixed frame
|
||||
|
||||
MathExtra::invquatvec(quat[ibody],conjqm[ibody],mbody);
|
||||
MathExtra::matvec(ex_space[ibody],ey_space[ibody],ez_space[ibody],
|
||||
mbody,angmom[ibody]);
|
||||
|
||||
angmom[ibody][0] *= 0.5;
|
||||
angmom[ibody][1] *= 0.5;
|
||||
angmom[ibody][2] *= 0.5;
|
||||
|
||||
// compute new angular velocity
|
||||
|
||||
MathExtra::angmom_to_omega(angmom[ibody],ex_space[ibody],ey_space[ibody],
|
||||
ez_space[ibody],inertia[ibody],omega[ibody]);
|
||||
}
|
||||
|
||||
// set velocity/rotation of atoms in rigid bodies
|
||||
// virial is already setup from initial_integrate
|
||||
|
||||
set_v();
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixRigidNVT::update_nhcp(double akin_t, double akin_r)
|
||||
{
|
||||
int i,j,k;
|
||||
double kt,gfkt_t,gfkt_r,tmp,ms,s,s2;
|
||||
|
||||
kt = boltz * t_target;
|
||||
gfkt_t = nf_t * kt;
|
||||
gfkt_r = nf_r * kt;
|
||||
|
||||
akin_t *= force->mvv2e;
|
||||
akin_r *= force->mvv2e;
|
||||
|
||||
// update thermostat masses
|
||||
|
||||
double t_mass = boltz * t_target / (t_freq * t_freq);
|
||||
q_t[0] = nf_t * t_mass;
|
||||
q_r[0] = nf_r * t_mass;
|
||||
for (i = 1; i < t_chain; i++)
|
||||
q_t[i] = q_r[i] = t_mass;
|
||||
|
||||
// update order/timestep dependent coefficients
|
||||
|
||||
for (i = 0; i < t_order; i++) {
|
||||
wdti1[i] = w[i] * dtv / t_iter;
|
||||
wdti2[i] = wdti1[i] / 2.0;
|
||||
wdti4[i] = wdti1[i] / 4.0;
|
||||
}
|
||||
|
||||
// update force of thermostats coupled to particles
|
||||
|
||||
f_eta_t[0] = (akin_t - gfkt_t) / q_t[0];
|
||||
f_eta_r[0] = (akin_r - gfkt_r) / q_r[0];
|
||||
|
||||
// multiple timestep iteration
|
||||
|
||||
for (i = 0; i < t_iter; i++) {
|
||||
for (j = 0; j < t_order; j++) {
|
||||
|
||||
// update thermostat velocities half step
|
||||
|
||||
eta_dot_t[t_chain-1] += wdti2[j] * f_eta_t[t_chain-1];
|
||||
eta_dot_r[t_chain-1] += wdti2[j] * f_eta_r[t_chain-1];
|
||||
|
||||
for (k = 1; k < t_chain; k++) {
|
||||
tmp = wdti4[j] * eta_dot_t[t_chain-k];
|
||||
ms = maclaurin_series(tmp);
|
||||
s = exp(-1.0 * tmp);
|
||||
s2 = s * s;
|
||||
eta_dot_t[t_chain-k-1] = eta_dot_t[t_chain-k-1] * s2 +
|
||||
wdti2[j] * f_eta_t[t_chain-k-1] * s * ms;
|
||||
|
||||
tmp = wdti4[j] * eta_dot_r[t_chain-k];
|
||||
ms = maclaurin_series(tmp);
|
||||
s = exp(-1.0 * tmp);
|
||||
s2 = s * s;
|
||||
eta_dot_r[t_chain-k-1] = eta_dot_r[t_chain-k-1] * s2 +
|
||||
wdti2[j] * f_eta_r[t_chain-k-1] * s * ms;
|
||||
}
|
||||
|
||||
// update thermostat positions a full step
|
||||
|
||||
for (k = 0; k < t_chain; k++) {
|
||||
eta_t[k] += wdti1[j] * eta_dot_t[k];
|
||||
eta_r[k] += wdti1[j] * eta_dot_r[k];
|
||||
}
|
||||
|
||||
// update thermostat forces
|
||||
|
||||
for (k = 1; k < t_chain; k++) {
|
||||
f_eta_t[k] = q_t[k-1] * eta_dot_t[k-1] * eta_dot_t[k-1] - kt;
|
||||
f_eta_t[k] /= q_t[k];
|
||||
f_eta_r[k] = q_r[k-1] * eta_dot_r[k-1] * eta_dot_r[k-1] - kt;
|
||||
f_eta_r[k] /= q_r[k];
|
||||
}
|
||||
|
||||
// update thermostat velocities a full step
|
||||
|
||||
for (k = 0; k < t_chain-1; k++) {
|
||||
tmp = wdti4[j] * eta_dot_t[k+1];
|
||||
ms = maclaurin_series(tmp);
|
||||
s = exp(-1.0 * tmp);
|
||||
s2 = s * s;
|
||||
eta_dot_t[k] = eta_dot_t[k] * s2 + wdti2[j] * f_eta_t[k] * s * ms;
|
||||
tmp = q_t[k] * eta_dot_t[k] * eta_dot_t[k] - kt;
|
||||
f_eta_t[k+1] = tmp / q_t[k+1];
|
||||
|
||||
tmp = wdti4[j] * eta_dot_r[k+1];
|
||||
ms = maclaurin_series(tmp);
|
||||
s = exp(-1.0 * tmp);
|
||||
s2 = s * s;
|
||||
eta_dot_r[k] = eta_dot_r[k] * s2 + wdti2[j] * f_eta_r[k] * s * ms;
|
||||
tmp = q_r[k] * eta_dot_r[k] * eta_dot_r[k] - kt;
|
||||
f_eta_r[k+1] = tmp / q_r[k+1];
|
||||
}
|
||||
|
||||
eta_dot_t[t_chain-1] += wdti2[j] * f_eta_t[t_chain-1];
|
||||
eta_dot_r[t_chain-1] += wdti2[j] * f_eta_r[t_chain-1];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
pack entire state of Fix into one write
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void FixRigidNVT::write_restart(FILE *fp)
|
||||
{
|
||||
int n = 0;
|
||||
double *list = new double[1+6*t_chain];
|
||||
list[n++] = t_chain;
|
||||
|
||||
for (int i = 0; i < t_chain; i++) {
|
||||
list[n++] = eta_t[i];
|
||||
list[n++] = eta_r[i];
|
||||
list[n++] = eta_dot_t[i];
|
||||
list[n++] = eta_dot_r[i];
|
||||
list[n++] = f_eta_t[i];
|
||||
list[n++] = f_eta_r[i];
|
||||
}
|
||||
|
||||
if (comm->me == 0) {
|
||||
int size = (1 + 6*t_chain)*sizeof(double);
|
||||
fwrite(&size,sizeof(int),1,fp);
|
||||
fwrite(list,sizeof(double),1+6*t_chain,fp);
|
||||
}
|
||||
|
||||
delete list;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
compute kinetic energy in the extended Hamiltonian
|
||||
conserved quantity = sum of returned energy and potential energy
|
||||
-----------------------------------------------------------------------*/
|
||||
|
||||
double FixRigidNVT::compute_scalar()
|
||||
{
|
||||
int i,k,ibody;
|
||||
double kt = boltz * t_target;
|
||||
double energy,ke_t,ke_q,tmp,Pkq[4];
|
||||
|
||||
// compute the kinetic parts of H_NVE in Kameraj et al (JCP 2005, pp 224114)
|
||||
|
||||
// translational kinetic energy
|
||||
|
||||
ke_t = 0.0;
|
||||
for (ibody = 0; ibody < nbody; ibody++)
|
||||
ke_t += 0.5 * masstotal[ibody] * (vcm[ibody][0]*vcm[ibody][0] +
|
||||
vcm[ibody][1]*vcm[ibody][1] +
|
||||
vcm[ibody][2]*vcm[ibody][2]);
|
||||
|
||||
// rotational kinetic energy
|
||||
|
||||
ke_q = 0.0;
|
||||
for (ibody = 0; ibody < nbody; ibody++) {
|
||||
for (k = 1; k < 4; k++) {
|
||||
if (k == 1) {
|
||||
Pkq[0] = -quat[ibody][1];
|
||||
Pkq[1] = quat[ibody][0];
|
||||
Pkq[2] = quat[ibody][3];
|
||||
Pkq[3] = -quat[ibody][2];
|
||||
} else if (k == 2) {
|
||||
Pkq[0] = -quat[ibody][2];
|
||||
Pkq[1] = -quat[ibody][3];
|
||||
Pkq[2] = quat[ibody][0];
|
||||
Pkq[3] = quat[ibody][1];
|
||||
} else if (k == 3) {
|
||||
Pkq[0] = -quat[ibody][3];
|
||||
Pkq[1] = quat[ibody][2];
|
||||
Pkq[2] = -quat[ibody][1];
|
||||
Pkq[3] = quat[ibody][0];
|
||||
}
|
||||
|
||||
tmp = conjqm[ibody][0]*Pkq[0] + conjqm[ibody][1]*Pkq[1] +
|
||||
conjqm[ibody][2]*Pkq[2] + conjqm[ibody][3]*Pkq[3];
|
||||
tmp *= tmp;
|
||||
|
||||
if (fabs(inertia[ibody][k-1]) < 1e-6) tmp = 0.0;
|
||||
else tmp /= (8.0 * inertia[ibody][k-1]);
|
||||
ke_q += tmp;
|
||||
}
|
||||
}
|
||||
|
||||
energy = ke_t + ke_q;
|
||||
|
||||
// thermostat chain energy: from equation 12 in Kameraj et al (JCP 2005)
|
||||
|
||||
energy += kt * (nf_t * eta_t[0] + nf_r * eta_r[0]);
|
||||
|
||||
for (i = 1; i < t_chain; i++)
|
||||
energy += kt * (eta_t[i] + eta_r[i]);
|
||||
|
||||
for (i = 0; i < t_chain; i++) {
|
||||
energy += 0.5 * q_t[i] * (eta_dot_t[i] * eta_dot_t[i]);
|
||||
energy += 0.5 * q_r[i] * (eta_dot_r[i] * eta_dot_r[i]);
|
||||
}
|
||||
|
||||
return energy;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
use state info from restart file to restart the Fix
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void FixRigidNVT::restart(char *buf)
|
||||
{
|
||||
int n = 0;
|
||||
double *list = (double *) buf;
|
||||
|
||||
int t_chain_prev = static_cast<int> (list[n++]);
|
||||
if (t_chain_prev != t_chain)
|
||||
error->all(FLERR,"Cannot restart fix rigid/nvt with different # of chains");
|
||||
|
||||
for (int i = 0; i < t_chain; i++) {
|
||||
eta_t[i] = list[n++];
|
||||
eta_r[i] = list[n++];
|
||||
eta_dot_t[i] = list[n++];
|
||||
eta_dot_r[i] = list[n++];
|
||||
f_eta_t[i] = list[n++];
|
||||
f_eta_r[i] = list[n++];
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixRigidNVT::reset_target(double t_new)
|
||||
{
|
||||
t_start = t_stop = t_new;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixRigidNVT::allocate_chain()
|
||||
{
|
||||
q_t = new double[t_chain];
|
||||
q_r = new double[t_chain];
|
||||
eta_t = new double[t_chain];
|
||||
eta_r = new double[t_chain];
|
||||
eta_dot_t = new double[t_chain];
|
||||
eta_dot_r = new double[t_chain];
|
||||
f_eta_t = new double[t_chain];
|
||||
f_eta_r = new double[t_chain];
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixRigidNVT::allocate_order()
|
||||
{
|
||||
w = new double[t_order];
|
||||
wdti1 = new double[t_order];
|
||||
wdti2 = new double[t_order];
|
||||
wdti4 = new double[t_order];
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixRigidNVT::deallocate_chain()
|
||||
{
|
||||
delete [] q_t;
|
||||
delete [] q_r;
|
||||
delete [] eta_t;
|
||||
delete [] eta_r;
|
||||
delete [] eta_dot_t;
|
||||
delete [] eta_dot_r;
|
||||
delete [] f_eta_t;
|
||||
delete [] f_eta_r;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixRigidNVT::deallocate_order()
|
||||
{
|
||||
delete [] w;
|
||||
delete [] wdti1;
|
||||
delete [] wdti2;
|
||||
delete [] wdti4;
|
||||
}
|
||||
|
|
|
|||
|
|
@ -20,85 +20,17 @@ FixStyle(rigid/nvt,FixRigidNVT)
|
|||
#ifndef LMP_FIX_RIGID_NVT_H
|
||||
#define LMP_FIX_RIGID_NVT_H
|
||||
|
||||
#include "fix_rigid.h"
|
||||
#include "fix_rigid_nh.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class FixRigidNVT : public FixRigid {
|
||||
class FixRigidNVT : public FixRigidNH {
|
||||
public:
|
||||
FixRigidNVT(class LAMMPS *, int, char **);
|
||||
~FixRigidNVT();
|
||||
int setmask();
|
||||
void init();
|
||||
void setup(int);
|
||||
void initial_integrate(int);
|
||||
void final_integrate();
|
||||
double compute_scalar();
|
||||
void write_restart(FILE *);
|
||||
void restart(char *);
|
||||
void reset_target(double);
|
||||
|
||||
private:
|
||||
double **conjqm; // conjugate quaternion momentum
|
||||
double boltz; // boltzman constant
|
||||
double t_target;
|
||||
|
||||
int nf_t,nf_r; // trans/rot degrees of freedom
|
||||
double *w,*wdti1,*wdti2,*wdti4; // Yoshida-Suzuki coefficients
|
||||
double *q_t,*q_r; // trans/rot thermostat masses
|
||||
double *eta_t,*eta_r; // trans/rot thermostat positions
|
||||
double *eta_dot_t,*eta_dot_r; // trans/rot thermostat velocities
|
||||
double *f_eta_t,*f_eta_r; // trans/rot thermostat forces
|
||||
|
||||
void update_nhcp(double, double);
|
||||
inline double maclaurin_series(double);
|
||||
|
||||
void allocate_chain();
|
||||
void allocate_order();
|
||||
void deallocate_chain();
|
||||
void deallocate_order();
|
||||
~FixRigidNVT() {}
|
||||
};
|
||||
|
||||
inline double FixRigidNVT::maclaurin_series(double x)
|
||||
{
|
||||
double x2,x4;
|
||||
x2 = x * x;
|
||||
x4 = x2 * x2;
|
||||
return (1.0 + (1.0/6.0) * x2 + (1.0/120.0) * x4 + (1.0/5040.0) * x2 * x4 +
|
||||
(1.0/362880.0) * x4 * x4);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
||||
/* ERROR/WARNING messages:
|
||||
|
||||
E: Did not set temp for fix rigid/nvt
|
||||
|
||||
The temp keyword must be used.
|
||||
|
||||
E: Target temperature for fix rigid/nvt cannot be 0.0
|
||||
|
||||
Self-explanatory.
|
||||
|
||||
E: Fix rigid/nvt period must be > 0.0
|
||||
|
||||
Self-explanatory
|
||||
|
||||
E: Illegal ... command
|
||||
|
||||
Self-explanatory. Check the input script syntax and compare to the
|
||||
documentation for the command. You can use -echo screen as a
|
||||
command-line option when running LAMMPS to see the offending line.
|
||||
|
||||
E: Fix_modify order must be 3 or 5
|
||||
|
||||
Self-explanatory.
|
||||
|
||||
E: Cannot restart fix rigid/nvt with different # of chains
|
||||
|
||||
This is because the restart file contains per-chain info.
|
||||
|
||||
*/
|
||||
|
|
|
|||
Loading…
Reference in New Issue