1 /* -*- c++ -*- ---------------------------------------------------------- 2 LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator 3 https://www.lammps.org/, Sandia National Laboratories 4 Steve Plimpton, sjplimp@sandia.gov 5 6 Copyright (2003) Sandia Corporation. Under the terms of Contract 7 DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains 8 certain rights in this software. This software is distributed under 9 the GNU General Public License. 10 11 See the README file in the top-level LAMMPS directory. 12 ------------------------------------------------------------------------- */ 13 14 #ifdef FIX_CLASS 15 // clang-format off 16 FixStyle(npt/cauchy,FixNPTCauchy); 17 // clang-format on 18 #else 19 20 #ifndef LMP_FIX_NPT_CAUCHY_H 21 #define LMP_FIX_NPT_CAUCHY_H 22 23 #include "fix.h" 24 25 namespace LAMMPS_NS { 26 27 class FixNPTCauchy : public Fix { 28 public: 29 FixNPTCauchy(class LAMMPS *, int, char **); 30 virtual void post_constructor(); 31 virtual ~FixNPTCauchy(); 32 int setmask(); 33 virtual void init(); 34 virtual void setup(int); 35 virtual void initial_integrate(int); 36 virtual void final_integrate(); 37 void initial_integrate_respa(int, int, int); 38 void pre_force_respa(int, int, int); 39 void final_integrate_respa(int, int); 40 virtual void pre_exchange(); 41 double compute_scalar(); 42 virtual double compute_vector(int); 43 void write_restart(FILE *); 44 virtual int pack_restart_data(double *); // pack restart data 45 virtual void restart(char *); 46 int modify_param(int, char **); 47 void reset_target(double); 48 void reset_dt(); 49 virtual void *extract(const char *, int &); 50 double memory_usage(); 51 52 protected: 53 int dimension, which; 54 double dtv, dtf, dthalf, dt4, dt8, dto; 55 double boltz, nktv2p, tdof; 56 double vol0; // reference volume 57 double t0; // reference temperature 58 // used for barostat mass 59 double t_start, t_stop; 60 double t_current, t_target, ke_target; 61 double t_freq; 62 63 int tstat_flag; // 1 if control T 64 int pstat_flag; // 1 if control P 65 66 int pstyle, pcouple, allremap; 67 int p_flag[6]; // 1 if control P on this dim, 0 if not 68 double p_start[6], p_stop[6]; 69 double p_freq[6], p_target[6]; 70 double omega[6], omega_dot[6]; 71 double omega_mass[6]; 72 double p_current[6]; 73 double drag, tdrag_factor; // drag factor on particle thermostat 74 double pdrag_factor; // drag factor on barostat 75 int kspace_flag; // 1 if KSpace invoked, 0 if not 76 int nrigid; // number of rigid fixes 77 int dilate_group_bit; // mask for dilation group 78 int *rfix; // indices of rigid fixes 79 char *id_dilate; // group name to dilate 80 class Irregular *irregular; // for migrating atoms after box flips 81 82 int nlevels_respa; 83 double *step_respa; 84 85 char *id_temp, *id_press; 86 class Compute *temperature, *pressure; 87 int tcomputeflag, pcomputeflag; // 1 = compute was created by fix 88 // 0 = created externally 89 90 double *eta, *eta_dot; // chain thermostat for particles 91 double *eta_dotdot; 92 double *eta_mass; 93 int mtchain; // length of chain 94 int mtchain_default_flag; // 1 = mtchain is default 95 96 double *etap; // chain thermostat for barostat 97 double *etap_dot; 98 double *etap_dotdot; 99 double *etap_mass; 100 int mpchain; // length of chain 101 102 int mtk_flag; // 0 if using Hoover barostat 103 int pdim; // number of barostatted dims 104 double p_freq_max; // maximum barostat frequency 105 106 double p_hydro; // hydrostatic target pressure 107 108 int nc_tchain, nc_pchain; 109 double factor_eta; 110 double sigma[6]; // scaled target stress 111 double fdev[6]; // deviatoric force on barostat 112 int deviatoric_flag; // 0 if target stress tensor is hydrostatic 113 double h0_inv[6]; // h_inv of reference (zero strain) box 114 int nreset_h0; // interval for resetting h0 115 116 double mtk_term1, mtk_term2; // Martyna-Tobias-Klein corrections 117 118 int eta_mass_flag; // 1 if eta_mass updated, 0 if not. 119 int omega_mass_flag; // 1 if omega_mass updated, 0 if not. 120 int etap_mass_flag; // 1 if etap_mass updated, 0 if not. 121 int dipole_flag; // 1 if dipole is updated, 0 if not. 122 int dlm_flag; // 1 if using the DLM rotational integrator, 0 if not 123 124 int scaleyz; // 1 if yz scaled with lz 125 int scalexz; // 1 if xz scaled with lz 126 int scalexy; // 1 if xy scaled with ly 127 int flipflag; // 1 if box flips are invoked as needed 128 129 int pre_exchange_flag; // set if pre_exchange needed for box flips 130 131 double fixedpoint[3]; // location of dilation fixed-point 132 133 void couple(); 134 virtual void remap(); 135 void nhc_temp_integrate(); 136 void nhc_press_integrate(); 137 138 virtual void nve_x(); // may be overwritten by child classes 139 virtual void nve_v(); 140 virtual void nh_v_press(); 141 virtual void nh_v_temp(); 142 virtual void compute_temp_target(); 143 virtual int size_restart_global(); 144 145 void compute_sigma(); 146 void compute_deviatoric(); 147 double compute_strain_energy(); 148 void compute_press_target(); 149 void nh_omega_dot(); 150 151 // Implementation of CauchyStat 152 char *id_store; // fix id of the STORE fix for retaining data 153 class FixStore *init_store; // fix pointer to STORE fix 154 double H0[3][3]; // shape matrix for the undeformed cell 155 double h_old[6]; // previous time step shape matrix for 156 // the undeformed cell 157 double invH0[3][3]; // inverse of H0; 158 double CSvol0; // volume of undeformed cell 159 double setPK[3][3]; // current set values of the PK stress 160 // (this is modified until the cauchy 161 // stress converges) 162 double alpha; // integration parameter for the cauchystat 163 int initPK; // 1 if setPK needs to be initialized either 164 // from cauchy or restart, else 0 165 int restartPK; // Read PK stress from the previous run 166 int restart_stored; // values of PK stress from the previous step stored 167 int initRUN; // 0 if run not initialized 168 // (pressure->vector not computed yet), 169 // else 1 (pressure->vector available) 170 171 void CauchyStat_init(); 172 void CauchyStat_cleanup(); 173 void CauchyStat(); 174 void CauchyStat_Step(double (&Fi)[3][3], double (&Fdot)[3][3], double (&cauchy)[3][3], 175 double (&setcauchy)[3][3], double (&setPK)[3][3], double volume, 176 double volume0, double deltat, double alpha); 177 }; 178 179 } // namespace LAMMPS_NS 180 181 #endif 182 #endif 183 184 /* ERROR/WARNING messages: 185 186 E: Illegal ... command 187 188 Self-explanatory. Check the input script syntax and compare to the 189 documentation for the command. You can use -echo screen as a 190 command-line option when running LAMMPS to see the offending line. 191 192 E: Target temperature for fix npt/cauchy cannot be 0.0 193 194 Self-explanatory. 195 196 E: Invalid fix npt/cauchy command for a 2d simulation 197 198 Cannot control z dimension in a 2d model. 199 200 E: Fix npt/cauchy dilate group ID does not exist 201 202 Self-explanatory. 203 204 E: Invalid fix npt/cauchy command pressure settings 205 206 If multiple dimensions are coupled, those dimensions must be 207 specified. 208 209 E: Cannot use fix npt/cauchy on a non-periodic dimension 210 211 When specifying a diagonal pressure component, the dimension must be 212 periodic. 213 214 E: Cannot use fix npt/cauchy on a 2nd non-periodic dimension 215 216 When specifying an off-diagonal pressure component, the 2nd of the two 217 dimensions must be periodic. E.g. if the xy component is specified, 218 then the y dimension must be periodic. 219 220 E: Cannot use fix npt/cauchy with yz scaling when z is non-periodic dimension 221 222 The 2nd dimension in the barostatted tilt factor must be periodic. 223 224 E: Cannot use fix npt/cauchy with xz scaling when z is non-periodic dimension 225 226 The 2nd dimension in the barostatted tilt factor must be periodic. 227 228 E: Cannot use fix npt/cauchy with xy scaling when y is non-periodic dimension 229 230 The 2nd dimension in the barostatted tilt factor must be periodic. 231 232 E: Cannot use fix npt/cauchy with both yz dynamics and yz scaling 233 234 Self-explanatory. 235 236 E: Cannot use fix npt/cauchy with both xz dynamics and xz scaling 237 238 Self-explanatory. 239 240 E: Cannot use fix npt/cauchy with both xy dynamics and xy scaling 241 242 Self-explanatory. 243 244 E: Can not specify Pxy/Pxz/Pyz in fix npt/cauchy with non-triclinic box 245 246 Only triclinic boxes can be used with off-diagonal pressure components. 247 See the region prism command for details. 248 249 E: Invalid fix npt/cauchy pressure settings 250 251 Settings for coupled dimensions must be the same. 252 253 E: Using update dipole flag requires atom style sphere 254 255 Self-explanatory. 256 257 E: Using update dipole flag requires atom attribute mu 258 259 Self-explanatory. 260 261 E: Fix npt/cauchy damping parameters must be > 0.0 262 263 Self-explanatory. 264 265 E: Cannot use fix npt/cauchy and fix deform on same component of stress tensor 266 267 This would be changing the same box dimension twice. 268 269 E: Temperature ID for fix npt/cauchy does not exist 270 271 Self-explanatory. 272 273 E: Pressure ID for fix npt/cauchy does not exist 274 275 Self-explanatory. 276 277 E: Non-numeric pressure - simulation unstable 278 279 UNDOCUMENTED 280 281 E: Fix npt/cauchy has tilted box too far in one step - periodic cell is too far from equilibrium state 282 283 Self-explanatory. The change in the box tilt is too extreme 284 on a short timescale. 285 286 E: Could not find fix_modify temperature ID 287 288 The compute ID for computing temperature does not exist. 289 290 E: Fix_modify temperature ID does not compute temperature 291 292 The compute ID assigned to the fix must compute temperature. 293 294 W: Temperature for fix modify is not for group all 295 296 The temperature compute is being used with a pressure calculation 297 which does operate on group all, so this may be inconsistent. 298 299 E: Pressure ID for fix modify does not exist 300 301 Self-explanatory. 302 303 E: Could not find fix_modify pressure ID 304 305 The compute ID for computing pressure does not exist. 306 307 E: Fix_modify pressure ID does not compute pressure 308 309 The compute ID assigned to the fix must compute pressure. 310 311 U: The dlm flag must be used with update dipole 312 313 Self-explanatory. 314 315 */ 316