git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@11837 f3b2605a-c512-4ea7-a41b-209d697bcdaa

src/RIGID/fix_rigid_nh_small.h

@@ -0,0 +1,195 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifndef LMP_FIX_RIGID_NH_SMALL_H
+#define LMP_FIX_RIGID_NH_SMALL_H
+
+#include "fix_rigid_small.h"
+
+namespace LAMMPS_NS {
+
+class FixRigidNHSmall : public FixRigidSmall {
+ public:
+  FixRigidNHSmall(class LAMMPS *, int, char **);
+  virtual ~FixRigidNHSmall();
+  virtual int setmask();
+  virtual void init();
+  virtual void setup(int);
+  virtual void initial_integrate(int);
+  virtual void final_integrate();
+  virtual double compute_scalar();
+  int modify_param(int, char **);
+  void write_restart(FILE *);
+  void restart(char *buf);
+  void reset_target(double);
+  
+ protected:
+  double boltz,nktv2p,mvv2e;          // boltzman constant, conversion factors
+
+  int dimension;                      // # of dimensions  
+  int nf_t,nf_r;                      // trans/rot degrees of freedom
+  double onednft,onednfr;             // factors 1 + dimension/trans(rot) degrees of freedom
+  double *w,*wdti1,*wdti2,*wdti4;     // Yoshida-Suzuki coefficients
+  double *q_t,*q_r;                   // trans/rot thermostat masses
+  double *eta_t,*eta_r;               // trans/rot thermostat positions
+  double *eta_dot_t,*eta_dot_r;       // trans/rot thermostat velocities
+  double *f_eta_t,*f_eta_r;           // trans/rot thermostat forces
+  
+  double epsilon_mass[3], *q_b;       // baro/thermo masses
+  double epsilon[3],*eta_b;           // baro/thermo positions
+  double epsilon_dot[3],*eta_dot_b;   // baro/thermo velocities
+  double *f_eta_b;                    // thermo forces
+  double akin_t,akin_r;               // translational/rotational kinetic energies
+  
+  int kspace_flag;                    // 1 if KSpace invoked, 0 if not
+  int nrigidfix;                      // number of rigid fixes
+  int *rfix;                          // indicies of rigid fixes
+
+  double vol0;                        // reference volume
+  double t0;                          // reference temperature
+  int pdim,g_f;                       // number of barostatted dims, total DoFs
+  double p_hydro;                     // hydrostatic target pressure
+  double p_freq_max;                  // maximum barostat frequency
+  
+  double mtk_term1,mtk_term2;         // Martyna-Tobias-Klein corrections
+
+  double t_target,t_current;
+  double t_freq;
+  
+  char *id_temp,*id_press;
+  class Compute *temperature,*pressure;
+  int tcomputeflag,pcomputeflag;
+
+  void couple();
+  void remap();  
+  void nhc_temp_integrate();
+  void nhc_press_integrate();
+
+  virtual void compute_temp_target();
+  void compute_press_target();
+  void nh_epsilon_dot();
+
+  void allocate_chain();
+  void allocate_order();
+  void deallocate_chain();
+  void deallocate_order();
+
+  void no_squish_rotate(int, double *, double *, double *, double);
+  inline double maclaurin_series(double);
+};
+
+inline double FixRigidNHSmall::maclaurin_series(double x)
+{
+  double x2,x4;
+  x2 = x * x;
+  x4 = x2 * x2;
+  return (1.0 + (1.0/6.0) * x2 + (1.0/120.0) * x4 + (1.0/5040.0) * x2 * x4 +
+         (1.0/362880.0) * x4 * x4);
+}
+
+}
+
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Fix rigid npt/nph period must be > 0.0
+
+Self-explanatory.
+
+E: Invalid fix rigid npt/nph command for a 2d simulation
+
+Cannot control z dimension in a 2d model.
+
+E: Invalid fix rigid npt/nph command pressure settings
+
+If multiple dimensions are coupled, those dimensions must be
+specified.
+
+E: Cannot use fix rigid npt/nph on a non-periodic dimension
+
+When specifying a diagonal pressure component, the dimension must be
+periodic.
+
+E: Invalid fix rigid npt/nph pressure settings
+
+Settings for coupled dimensions must be the same.
+
+E: Fix rigid nvt/npt/nph damping parameters must be > 0.0
+
+Self-explanatory.
+
+E: Fix rigid npt/nph dilate group ID does not exist
+
+Self-explanatory.
+
+E: Temp ID for fix rigid npt/nph does not exist
+
+Specified compute temperature must be valid.
+
+E: fix rigid npt/nph does not yet allow triclinic box
+
+Self-explanatory.
+
+E: Cannot use fix rigid npt/nph and fix deform on same component of stress tensor
+
+This would be changing the same box dimension twice.
+
+E: Press ID for fix rigid npt/nph does not exist
+
+Specified compute pressure must be valid.
+
+E: Illegal ... command
+
+Self-explanatory.  Check the input script syntax and compare to the
+documentation for the command.  You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+E: Could not find fix_modify temperature ID
+
+The compute ID for computing temperature does not exist.
+
+E: Fix_modify temperature ID does not compute temperature
+
+The compute ID assigned to the fix must compute temperature.
+
+W: Temperature for fix modify is not for group all
+
+The temperature compute is being used with a pressure calculation
+which does operate on group all, so this may be inconsistent.
+
+E: Pressure ID for fix modify does not exist
+
+Self-explanatory.
+
+E: Could not find fix_modify pressure ID
+
+The compute ID for computing pressure does not exist.
+
+E: Fix_modify pressure ID does not compute pressure
+
+The compute ID assigned to the fix must compute pressure.
+
+U: Target temperature for fix rigid nvt/npt cannot be 0.0
+
+Self-explanatory.
+
+U: Temperature ID for fix rigid npt/nph does not exist
+
+Self-explanatory.
+
+U: Pressure ID for fix rigid npt/nph does not exist
+
+Self-explanatory.
+
+*/

src/RIGID/fix_rigid_nph_small.cpp

@@ -0,0 +1,92 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+   Contributing author: Trung Dac Nguyen (ORNL)
+   references: Kamberaj et al., J. Chem. Phys. 122, 224114 (2005)
+               Miller et al., J Chem Phys. 116, 8649-8659 (2002)
+------------------------------------------------------------------------- */
+
+#include "string.h"
+#include "fix_rigid_nph_small.h"
+#include "domain.h"
+#include "modify.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+FixRigidNPHSmall::FixRigidNPHSmall(LAMMPS *lmp, int narg, char **arg) :
+  FixRigidNHSmall(lmp, narg, arg)
+{ 
+  // other setting are made by parent
+
+  scalar_flag = 1;
+  restart_global = 1;
+  box_change_size = 1;
+  extscalar = 1;
+    
+  // error checks
+
+  if (pstat_flag == 0)
+    error->all(FLERR,"Pressure control must be used with fix nph");
+  if (tstat_flag == 1)
+    error->all(FLERR,"Temperature control must not be used with fix nph");
+  if (p_start[0] < 0.0 || p_start[1] < 0.0 || p_start[2] < 0.0 || 
+      p_stop[0] < 0.0 || p_stop[1] < 0.0 || p_stop[2] < 0.0)
+    error->all(FLERR,"Target pressure for fix rigid/nph cannot be 0.0");
+  
+  // convert input periods to frequency
+  p_freq[0] = p_freq[1] = p_freq[2] = 0.0;
+
+  if (p_flag[0]) p_freq[0] = 1.0 / p_period[0];
+  if (p_flag[1]) p_freq[1] = 1.0 / p_period[1];
+  if (p_flag[2]) p_freq[2] = 1.0 / p_period[2];
+
+  // create a new compute temp style
+  // id = fix-ID + temp
+  // compute group = all since pressure is always global (group all)
+  //   and thus its KE/temperature contribution should use group all
+
+  int n = strlen(id) + 6;
+  id_temp = new char[n];
+  strcpy(id_temp,id);
+  strcat(id_temp,"_temp");
+
+  char **newarg = new char*[3];
+  newarg[0] = id_temp;
+  newarg[1] = (char *) "all";
+  newarg[2] = (char *) "temp";  
+  modify->add_compute(3,newarg);
+  delete [] newarg;
+  tcomputeflag = 1;
+  
+  // create a new compute pressure style
+  // id = fix-ID + press, compute group = all
+  // pass id_temp as 4th arg to pressure constructor
+  
+  n = strlen(id) + 7;
+  id_press = new char[n];
+  strcpy(id_press,id);
+  strcat(id_press,"_press");
+
+  newarg = new char*[4];
+  newarg[0] = id_press;
+  newarg[1] = (char *) "all";
+  newarg[2] = (char *) "pressure";
+  newarg[3] = id_temp;
+  modify->add_compute(4,newarg);
+  delete [] newarg;
+  pcomputeflag = 1;
+}

src/RIGID/fix_rigid_nph_small.h

@@ -0,0 +1,53 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef FIX_CLASS
+
+FixStyle(rigid/nph/small,FixRigidNPHSmall)
+
+#else
+
+#ifndef LMP_FIX_RIGID_NPH_SMALL_H
+#define LMP_FIX_RIGID_NPH_SMALL_H
+
+#include "fix_rigid_nh_small.h"
+
+namespace LAMMPS_NS {
+
+class FixRigidNPHSmall : public FixRigidNHSmall {
+ public:
+  FixRigidNPHSmall(class LAMMPS *, int, char **);
+  ~FixRigidNPHSmall() {}
+};
+
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Pressure control must be used with fix rigid nph/small
+
+UNDOCUMENTED
+
+E: Temperature control must not be used with fix rigid/nph/small
+
+UNDOCUMENTED
+
+E: Target pressure for fix rigid/nph/small cannot be 0.0
+
+UNDOCUMENTED
+
+*/

src/RIGID/fix_rigid_npt_small.cpp

@@ -0,0 +1,104 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+   Contributing author: Trung Dac Nguyen (ORNL)
+   references: Kamberaj et al., J. Chem. Phys. 122, 224114 (2005)
+               Miller et al., J Chem Phys. 116, 8649-8659 (2002)
+------------------------------------------------------------------------- */
+
+#include "string.h"
+#include "fix_rigid_npt_small.h"
+#include "domain.h"
+#include "modify.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+FixRigidNPTSmall::FixRigidNPTSmall(LAMMPS *lmp, int narg, char **arg) :
+  FixRigidNHSmall(lmp, narg, arg)
+{ 
+  // other setting are made by parent
+
+  scalar_flag = 1;
+  restart_global = 1;
+  box_change_size = 1;
+  extscalar = 1;
+  
+  // error checks
+
+  if (tstat_flag == 0 || pstat_flag == 0)
+    error->all(FLERR,"Did not set temp or press for fix rigid/npt/small");
+  if (t_start <= 0.0 || t_stop <= 0.0)
+    error->all(FLERR,"Target temperature for fix rigid/npt/small cannot be 0.0");
+  if (p_start[0] < 0.0 || p_start[1] < 0.0 || p_start[2] < 0.0 || 
+      p_stop[0] < 0.0 || p_stop[1] < 0.0 || p_stop[2] < 0.0)
+    error->all(FLERR,"Target pressure for fix rigid/npt/small cannot be 0.0");
+
+  if (t_period <= 0.0) error->all(FLERR,"Fix rigid/npt/small period must be > 0.0");
+
+  // thermostat chain parameters
+
+  if (t_chain < 1) error->all(FLERR,"Illegal fix_modify command");
+  if (t_iter < 1) error->all(FLERR,"Illegal fix_modify command");
+  if (t_order != 3 && t_order != 5) 
+    error->all(FLERR,"Fix_modify order must be 3 or 5");
+
+  // convert input periods to frequency
+
+  t_freq = 0.0;
+  p_freq[0] = p_freq[1] = p_freq[2] = 0.0;
+
+  t_freq = 1.0 / t_period;
+  if (p_flag[0]) p_freq[0] = 1.0 / p_period[0];
+  if (p_flag[1]) p_freq[1] = 1.0 / p_period[1];
+  if (p_flag[2]) p_freq[2] = 1.0 / p_period[2];
+
+  // create a new compute temp style
+  // id = fix-ID + temp
+  // compute group = all since pressure is always global (group all)
+  //   and thus its KE/temperature contribution should use group all
+
+  int n = strlen(id) + 6;
+  id_temp = new char[n];
+  strcpy(id_temp,id);
+  strcat(id_temp,"_temp");
+
+  char **newarg = new char*[3];
+  newarg[0] = id_temp;
+  newarg[1] = (char *) "all";
+  newarg[2] = (char *) "temp";  
+  modify->add_compute(3,newarg);
+  delete [] newarg;
+  tcomputeflag = 1;
+  
+  // create a new compute pressure style
+  // id = fix-ID + press, compute group = all
+  // pass id_temp as 4th arg to pressure constructor
+  
+  n = strlen(id) + 7;
+  id_press = new char[n];
+  strcpy(id_press,id);
+  strcat(id_press,"_press");
+
+  newarg = new char*[4];
+  newarg[0] = id_press;
+  newarg[1] = (char *) "all";
+  newarg[2] = (char *) "pressure";
+  newarg[3] = id_temp;
+  modify->add_compute(4,newarg);
+  delete [] newarg;
+  pcomputeflag = 1;
+}

src/RIGID/fix_rigid_npt_small.h

@@ -0,0 +1,65 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef FIX_CLASS
+
+FixStyle(rigid/npt/small,FixRigidNPTSmall)
+
+#else
+
+#ifndef LMP_FIX_RIGID_NPT_SMALL_H
+#define LMP_FIX_RIGID_NPT_SMALL_H
+
+#include "fix_rigid_nh_small.h"
+
+namespace LAMMPS_NS {
+
+class FixRigidNPTSmall : public FixRigidNHSmall {
+ public:
+  FixRigidNPTSmall(class LAMMPS *, int, char **);
+  ~FixRigidNPTSmall() {}
+};
+
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Did not set temp or press for fix rigid/npt/small
+
+UNDOCUMENTED
+
+E: Target temperature for fix rigid/npt/small cannot be 0.0
+
+UNDOCUMENTED
+
+E: Target pressure for fix rigid/npt/small cannot be 0.0
+
+UNDOCUMENTED
+
+E: Fix rigid/npt/small period must be > 0.0
+
+UNDOCUMENTED
+
+E: Illegal ... command
+
+UNDOCUMENTED
+
+E: Fix_modify order must be 3 or 5
+
+UNDOCUMENTED
+
+*/

src/RIGID/fix_rigid_nve_small.cpp

@@ -0,0 +1,28 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+   Contributing author: Trung Dac Nguyen (ORNL)
+   references: Kamberaj et al., J. Chem. Phys. 122, 224114 (2005)
+               Miller et al., J Chem Phys. 116, 8649-8659 (2002)
+------------------------------------------------------------------------- */
+
+#include "fix_rigid_nve_small.h"
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+FixRigidNVESmall::FixRigidNVESmall(LAMMPS *lmp, int narg, char **arg) :
+  FixRigidNHSmall(lmp, narg, arg) {}
+

src/RIGID/fix_rigid_nve_small.h

@@ -0,0 +1,36 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef FIX_CLASS
+
+FixStyle(rigid/nve/small,FixRigidNVESmall)
+
+#else
+
+#ifndef LMP_FIX_RIGID_NVE_SMALL_H
+#define LMP_FIX_RIGID_NVE_SMALL_H
+
+#include "fix_rigid_nh_small.h"
+
+namespace LAMMPS_NS {
+
+class FixRigidNVESmall : public FixRigidNHSmall {
+ public:
+  FixRigidNVESmall(class LAMMPS *, int, char **);
+  ~FixRigidNVESmall() {}
+};
+
+}
+
+#endif
+#endif

src/RIGID/fix_rigid_nvt_small.cpp

@@ -0,0 +1,50 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+   Contributing author: Trung Dac Nguyen (ORNL)
+   references: Kamberaj et al., J. Chem. Phys. 122, 224114 (2005)
+               Miller et al., J Chem Phys. 116, 8649-8659 (2002)
+------------------------------------------------------------------------- */
+
+#include "fix_rigid_nvt_small.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+FixRigidNVTSmall::FixRigidNVTSmall(LAMMPS *lmp, int narg, char **arg) :
+  FixRigidNHSmall(lmp, narg, arg)
+{ 
+  // other settings are made by parent
+
+  scalar_flag = 1;
+  restart_global = 1;
+  extscalar = 1;
+  
+  // error checking
+  // convert input period to frequency
+
+  if (tstat_flag == 0)
+    error->all(FLERR,"Did not set temp for fix rigid/nvt/small");
+  if (t_start < 0.0 || t_stop <= 0.0)
+    error->all(FLERR,"Target temperature for fix rigid/nvt/small cannot be 0.0");
+  if (t_period <= 0.0) error->all(FLERR,"Fix rigid/nvt/small period must be > 0.0");
+  t_freq = 1.0 / t_period;
+
+  if (t_chain < 1) error->all(FLERR,"Fix rigid nvt/small t_chain should not be less than 1");
+  if (t_iter < 1) error->all(FLERR,"Fix rigid nvt/small t_iter should not be less than 1");
+  if (t_order != 3 && t_order != 5) 
+    error->all(FLERR,"Fix rigid nvt/small t_order must be 3 or 5"); 
+}

src/RIGID/fix_rigid_nvt_small.h

@@ -0,0 +1,60 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef FIX_CLASS
+
+FixStyle(rigid/nvt/small,FixRigidNVTSmall)
+
+#else
+
+#ifndef LMP_FIX_RIGID_NVT_SMALL_H
+#define LMP_FIX_RIGID_NVT_SMALL_H
+
+#include "fix_rigid_nh_small.h"
+
+namespace LAMMPS_NS {
+
+class FixRigidNVTSmall : public FixRigidNHSmall {
+ public:
+  FixRigidNVTSmall(class LAMMPS *, int, char **);
+  ~FixRigidNVTSmall() {}
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Did not set temp for fix rigid/nvt/small
+
+UNDOCUMENTED
+
+E: Target temperature for fix rigid/nvt/small cannot be 0.0
+
+UNDOCUMENTED
+
+E: Fix rigid/nvt/small period must be > 0.0
+
+UNDOCUMENTED
+
+E: Illegal ... command
+
+UNDOCUMENTED
+
+E: Fix_modify order must be 3 or 5
+
+UNDOCUMENTED
+
+*/

src/RIGID/fix_rigid_small.cpp

@@ -59,6 +59,9 @@ FixRigidSmall *FixRigidSmall::frsptr;
 
 #define DELTA_BODY 10000
 
+enum{NONE,XYZ,XY,YZ,XZ};        // same as in FixRigid
+enum{ISO,ANISO,TRICLINIC};      // same as in FixRigid
+
 enum{FULL_BODY,INITIAL,FINAL,FORCE_TORQUE,VCM_ANGMOM,XCM_MASS,ITENSOR,DOF};
 
 /* ---------------------------------------------------------------------- */
@@ -125,6 +128,23 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) :
   infile = NULL;
   onemol = NULL;
 
+  tstat_flag = 0;
+  pstat_flag = 0;
+  allremap = 1;
+  id_dilate = NULL;
+  t_chain = 10;
+  t_iter = 1;
+  t_order = 3;
+  p_chain = 10;
+
+  pcouple = NONE;
+  pstyle = ANISO;
+
+  for (int i = 0; i < 3; i++) {
+    p_start[i] = p_stop[i] = p_period[i] = 0.0;
+    p_flag[i] = 0;
+  }
+
   int iarg = 4;
   while (iarg < narg) {
     if (strcmp(arg[iarg],"langevin") == 0) {
@@ -159,6 +179,121 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) :
                        "fix rigid/small has multiple molecules");
       onemol = atom->molecules[imol];
       iarg += 2;
+
+    } else if (strcmp(arg[iarg],"temp") == 0) {
+      if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid/small command");
+      if (strcmp(style,"rigid/nvt/small") != 0 && 
+          strcmp(style,"rigid/npt/small") != 0)
+        error->all(FLERR,"Illegal fix rigid command");
+      tstat_flag = 1;
+      t_start = force->numeric(FLERR,arg[iarg+1]);
+      t_stop = force->numeric(FLERR,arg[iarg+2]);
+      t_period = force->numeric(FLERR,arg[iarg+3]);
+      iarg += 4;
+
+    } else if (strcmp(arg[iarg],"iso") == 0) {
+      if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid/small command");
+      if (strcmp(style,"rigid/npt/small") != 0 && 
+          strcmp(style,"rigid/nph/small") != 0)
+	      error->all(FLERR,"Illegal fix rigid/small command");
+      pcouple = XYZ;
+      p_start[0] = p_start[1] = p_start[2] = force->numeric(FLERR,arg[iarg+1]);
+      p_stop[0] = p_stop[1] = p_stop[2] = force->numeric(FLERR,arg[iarg+2]);
+      p_period[0] = p_period[1] = p_period[2] = 
+        force->numeric(FLERR,arg[iarg+3]);
+      p_flag[0] = p_flag[1] = p_flag[2] = 1;
+      if (domain->dimension == 2) {
+	      p_start[2] = p_stop[2] = p_period[2] = 0.0;
+      	p_flag[2] = 0;
+      }
+      iarg += 4;
+
+    } else if (strcmp(arg[iarg],"aniso") == 0) {
+      if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid/small command");
+      if (strcmp(style,"rigid/npt/small") != 0 && 
+          strcmp(style,"rigid/nph/small") != 0)
+	      error->all(FLERR,"Illegal fix rigid/small command");
+      p_start[0] = p_start[1] = p_start[2] = force->numeric(FLERR,arg[iarg+1]);
+      p_stop[0] = p_stop[1] = p_stop[2] = force->numeric(FLERR,arg[iarg+2]);
+      p_period[0] = p_period[1] = p_period[2] = 
+        force->numeric(FLERR,arg[iarg+3]);
+      p_flag[0] = p_flag[1] = p_flag[2] = 1;
+      if (domain->dimension == 2) {
+      	p_start[2] = p_stop[2] = p_period[2] = 0.0;
+	      p_flag[2] = 0;
+      }
+      iarg += 4;
+
+    } else if (strcmp(arg[iarg],"x") == 0) {
+      if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid/small command");
+      p_start[0] = force->numeric(FLERR,arg[iarg+1]);
+      p_stop[0] = force->numeric(FLERR,arg[iarg+2]);
+      p_period[0] = force->numeric(FLERR,arg[iarg+3]);
+      p_flag[0] = 1;
+      iarg += 4;
+
+    } else if (strcmp(arg[iarg],"y") == 0) {
+      if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid/small command");
+      p_start[1] = force->numeric(FLERR,arg[iarg+1]);
+      p_stop[1] = force->numeric(FLERR,arg[iarg+2]);
+      p_period[1] = force->numeric(FLERR,arg[iarg+3]);
+      p_flag[1] = 1;
+      iarg += 4;
+
+    } else if (strcmp(arg[iarg],"z") == 0) {
+      if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid/small command");
+      p_start[2] = force->numeric(FLERR,arg[iarg+1]);
+      p_stop[2] = force->numeric(FLERR,arg[iarg+2]);
+      p_period[2] = force->numeric(FLERR,arg[iarg+3]);
+      p_flag[2] = 1;
+      iarg += 4;
+
+    } else if (strcmp(arg[iarg],"couple") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal fix rigid/small command");
+      if (strcmp(arg[iarg+1],"xyz") == 0) pcouple = XYZ;
+      else if (strcmp(arg[iarg+1],"xy") == 0) pcouple = XY;
+      else if (strcmp(arg[iarg+1],"yz") == 0) pcouple = YZ;
+      else if (strcmp(arg[iarg+1],"xz") == 0) pcouple = XZ;
+      else if (strcmp(arg[iarg+1],"none") == 0) pcouple = NONE;
+      else error->all(FLERR,"Illegal fix rigid/small command");
+      iarg += 2;
+
+    } else if (strcmp(arg[iarg],"dilate") == 0) {
+      if (iarg+2 > narg) 
+        error->all(FLERR,"Illegal fix rigid/small nvt/npt/nph command");
+      if (strcmp(arg[iarg+1],"all") == 0) allremap = 1;
+      else {
+        allremap = 0;
+        delete [] id_dilate;
+        int n = strlen(arg[iarg+1]) + 1;
+        id_dilate = new char[n];
+        strcpy(id_dilate,arg[iarg+1]);
+        int idilate = group->find(id_dilate);
+        if (idilate == -1)
+          error->all(FLERR,"Fix rigid/small nvt/npt/nph dilate group ID "
+                     "does not exist");
+      }
+      iarg += 2;
+
+    } else if (strcmp(arg[iarg],"tparam") == 0) {
+      if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid/small command");
+      if (strcmp(style,"rigid/nvt/small") != 0 && 
+          strcmp(style,"rigid/npt/small") != 0)
+        error->all(FLERR,"Illegal fix rigid/small command");
+      t_chain = force->numeric(FLERR,arg[iarg+1]);
+      t_iter = force->numeric(FLERR,arg[iarg+2]);
+      t_order = force->numeric(FLERR,arg[iarg+3]);
+      iarg += 4;
+
+    } else if (strcmp(arg[iarg],"pchain") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal fix rigid/small command");
+      if (strcmp(style,"rigid/npt/small") != 0 && 
+          strcmp(style,"rigid/nph/small") != 0)
+        error->all(FLERR,"Illegal fix rigid/small command");
+      p_chain = force->numeric(FLERR,arg[iarg+1]);
+      iarg += 2;
+
+
     } else error->all(FLERR,"Illegal fix rigid/small command");
   }
 

src/RIGID/fix_rigid_small.h

@@ -97,7 +97,8 @@ class FixRigidSmall : public Fix {
     double ez_space[3];
     double angmom[3];         // space-frame angular momentum of body
     double omega[3];          // space-frame omega of body
-    imageint image;             // image flags of xcm
+    double conjqm[4];         // conjugate quaternion momentum
+    imageint image;           // image flags of xcm
     int remapflag[4];         // PBC remap flags
     int ilocal;               // index of owning atom
   };
@@ -151,6 +152,22 @@ class FixRigidSmall : public Fix {
   int maxlang;                      // max size of langextra
   class RanMars *random;            // RNG
 
+  int tstat_flag,pstat_flag;        // 0/1 = no/yes thermostat/barostat
+  
+  int t_chain,t_iter,t_order;
+
+  double p_start[3],p_stop[3];
+  double p_period[3],p_freq[3];
+  int p_flag[3];  
+  int pcouple,pstyle;
+  int p_chain;
+
+  int allremap;              // remap all atoms
+  int dilate_group_bit;      // mask for dilation group
+  char *id_dilate;           // group name to dilate
+
+  double p_current[3],p_target[3];
+
   // molecules added on-the-fly as rigid bodies
 
   class Molecule *onemol;