NAMD
SimParameters.C
Go to the documentation of this file.
1 
7 /*****************************************************************************
8  * $Source: /home/cvs/namd/cvsroot/namd2/src/SimParameters.C,v $
9  * $Author: jim $
10  * $Date: 2017/03/30 20:06:17 $
11  * $Revision: 1.1478 $
12  *****************************************************************************/
13 
20 #include "InfoStream.h"
21 #include "ComputeNonbondedUtil.h"
22 #include "LJTable.h"
23 #include "ComputePme.h"
24 #include "ComputeCUDAMgr.h"
25 #include "ConfigList.h"
26 #include "SimParameters.h"
27 #include "ParseOptions.h"
28 #include "structures.h"
29 #include "Communicate.h"
30 #include "MStream.h"
31 #include "Output.h"
32 #include <stdio.h>
33 #include <time.h>
34 #ifdef NAMD_FFTW
35 #ifdef NAMD_FFTW_3
36 #include <fftw3.h>
37 #else
38 // fftw2 doesn't have these defined
39 #define fftwf_malloc fftw_malloc
40 #define fftwf_free fftw_free
41 #ifdef NAMD_FFTW_NO_TYPE_PREFIX
42 #include <fftw.h>
43 #include <rfftw.h>
44 #else
45 #include <sfftw.h>
46 #include <srfftw.h>
47 #endif
48 #endif
49 #endif
50 #if defined(WIN32) && !defined(__CYGWIN__)
51 #include <direct.h>
52 #define CHDIR _chdir
53 #define MKDIR(X) mkdir(X)
54 #define PATHSEP '\\'
55 #define PATHSEPSTR "\\"
56 #else
57 #include <unistd.h>
58 #define CHDIR chdir
59 #define MKDIR(X) mkdir(X,0777)
60 #define PATHSEP '/'
61 #define PATHSEPSTR "/"
62 #endif
63 #include <fcntl.h>
64 #include <sys/stat.h>
65 #ifdef WIN32
66 #include <io.h>
67 #define access(PATH,MODE) _access(PATH,00)
68 #endif
69 #include <fstream>
70 using namespace std;
71 
72 #ifdef WIN32
73 extern "C" {
74  double erfc(double);
75 }
76 #endif
77 
78 #include "strlib.h" // For strcasecmp and strncasecmp
79 
81 
82 //#ifdef NAMD_CUDA
83 //bool one_cuda_device_per_node();
84 //#endif
85 #include "DeviceCUDA.h"
86 #if defined(NAMD_CUDA) || defined(NAMD_HIP)
87 #ifdef WIN32
88 #define __thread __declspec(thread)
89 #endif
90 extern __thread DeviceCUDA *deviceCUDA;
91 #endif
92 
93 #ifdef NAMD_AVXTILES
94 // For "+notiles" commandline option to disable "Tiles" algorithm (BackEnd.C)
95 extern int avxTilesCommandLineDisable;
96 #endif
97 
98 //#define DEBUGM
99 #include "Debug.h"
100 
101 #define XXXBIGREAL 1.0e32
102 
103 
104 char* SimParameters::getfromparseopts(const char* name, char *outbuf) {
105  if ( parseopts ) return parseopts->getfromptr(name,outbuf);
106  else return 0;
107 }
108 
109 int SimParameters::istrueinparseopts(const char* name) {
110  if ( parseopts ) return parseopts->istruefromptr(name);
111  else return -1;
112 }
113 
114 int SimParameters::issetinparseopts(const char* name) {
115  if ( parseopts ) return parseopts->issetfromptr(name);
116  else return -1;
117 }
118 
119 
120 /************************************************************************/
121 /* */
122 /* FUNCTION initialize_config_data */
123 /* */
124 /* This function is used by the master process to populate the */
125 /* simulation parameters from a ConfigList object that is passed to */
126 /* it. Each parameter is checked to make sure that it has a value */
127 /* that makes sense, and that it doesn't conflict with any other */
128 /* values that have been given. */
129 /* */
130 /************************************************************************/
131 
133 
134 {
135 
136  parseopts = new ParseOptions; // Object to check consistency of config file
137  ParseOptions &opts = *parseopts;
138 
139  config_parser(opts);
140 
142  if (!opts.check_consistency())
143  {
144  NAMD_die("Internal error in configuration file parser");
145  }
146 
147  // Now, feed the object with the actual configuration options through the
148  // ParseOptions file and make sure everything is OK
149  if (!opts.set(*config))
150  {
151  NAMD_die("ERROR(S) IN THE CONFIGURATION FILE");
152  }
153 
155 
156  check_config(opts,config,cwd);
157 
158  print_config(opts,config,cwd);
159 
160 }
161 
162 /************************************************************************/
163 /* */
164 /* FUNCTION scriptSet */
165 /* */
166 /************************************************************************/
167 
168 int atobool(const char *s) {
169  return ( (! strncasecmp(s,"yes",8)) ||
170  (! strncasecmp(s,"on",8)) ||
171  (! strncasecmp(s,"true",8)) );
172 };
173 
174 void SimParameters::scriptSet(const char *param, const char *value) {
175 
176  if ( CkMyRank() ) return;
177 
178 #define MAX_SCRIPT_PARAM_SIZE 128
179 #define SCRIPT_PARSE_BOOL(NAME,VAR) { if ( ! strncasecmp(param,(NAME),MAX_SCRIPT_PARAM_SIZE) ) { (VAR) = atobool(value); return; } }
180 #define SCRIPT_PARSE_INT(NAME,VAR) { if ( ! strncasecmp(param,(NAME),MAX_SCRIPT_PARAM_SIZE) ) { (VAR) = atoi(value); return; } }
181 #define SCRIPT_PARSE_FLOAT(NAME,VAR) { if ( ! strncasecmp(param,(NAME),MAX_SCRIPT_PARAM_SIZE) ) { (VAR) = atof(value); return; } }
182 #define SCRIPT_PARSE_MOD_FLOAT(NAME,VAR,MOD) { if ( ! strncasecmp(param,(NAME),MAX_SCRIPT_PARAM_SIZE) ) { (VAR) = atof(value) MOD; return; } }
183 #define SCRIPT_PARSE_VECTOR(NAME,VAR) { if ( ! strncasecmp(param,(NAME),MAX_SCRIPT_PARAM_SIZE) ) { (VAR).set(value); return; } }
184 #define SCRIPT_PARSE_STRING(NAME,VAR) { if ( ! strncasecmp(param,(NAME),MAX_SCRIPT_PARAM_SIZE) ) { strcpy(VAR,value); return; } }
185 
186  SCRIPT_PARSE_FLOAT("scriptArg1",scriptArg1)
187  SCRIPT_PARSE_FLOAT("scriptArg2",scriptArg2)
188  SCRIPT_PARSE_FLOAT("scriptArg3",scriptArg3)
189  SCRIPT_PARSE_FLOAT("scriptArg4",scriptArg4)
190  SCRIPT_PARSE_FLOAT("scriptArg5",scriptArg5)
191  SCRIPT_PARSE_INT("scriptIntArg1",scriptIntArg1)
192  SCRIPT_PARSE_INT("scriptIntArg2",scriptIntArg2)
193  SCRIPT_PARSE_STRING("scriptStringArg1",scriptStringArg1)
194  SCRIPT_PARSE_STRING("scriptStringArg2",scriptStringArg2)
195  SCRIPT_PARSE_INT("numsteps",N)
196  if ( ! strncasecmp(param,"firsttimestep",MAX_SCRIPT_PARAM_SIZE) ) {
197  N = firstTimestep = atoi(value); return;
198  }
199  SCRIPT_PARSE_FLOAT("reassignTemp",reassignTemp)
200  SCRIPT_PARSE_FLOAT("rescaleTemp",rescaleTemp)
201  SCRIPT_PARSE_BOOL("velocityQuenching",minimizeOn)
202  SCRIPT_PARSE_BOOL("maximumMove",maximumMove)
203  // SCRIPT_PARSE_BOOL("Langevin",langevinOn)
204  if ( ! strncasecmp(param,"Langevin",MAX_SCRIPT_PARAM_SIZE) ) {
205  langevinOn = atobool(value);
206  if ( langevinOn && ! langevinOnAtStartup ) {
207  NAMD_die("Langevin must be enabled at startup to disable and re-enable in script.");
208  }
209  return;
210  }
211  SCRIPT_PARSE_FLOAT("langevinTemp",langevinTemp)
212  SCRIPT_PARSE_BOOL("langevinBAOAB",langevin_useBAOAB) // [!!] Use the BAOAB integrator or not
213  SCRIPT_PARSE_FLOAT("loweAndersenTemp",loweAndersenTemp) // BEGIN LA, END LA
214  SCRIPT_PARSE_BOOL("stochRescale",stochRescaleOn)
215  SCRIPT_PARSE_FLOAT("stochRescaleTemp",stochRescaleTemp)
216  SCRIPT_PARSE_FLOAT("initialTemp",initialTemp)
217  SCRIPT_PARSE_BOOL("useGroupPressure",useGroupPressure)
218  SCRIPT_PARSE_BOOL("useFlexibleCell",useFlexibleCell)
219  SCRIPT_PARSE_BOOL("useConstantArea",useConstantArea)
220  SCRIPT_PARSE_BOOL("fixCellDims",fixCellDims)
221  SCRIPT_PARSE_BOOL("fixCellDimX",fixCellDimX)
222  SCRIPT_PARSE_BOOL("fixCellDimY",fixCellDimY)
223  SCRIPT_PARSE_BOOL("fixCellDimZ",fixCellDimZ)
224  SCRIPT_PARSE_BOOL("useConstantRatio",useConstantRatio)
225  SCRIPT_PARSE_BOOL("LangevinPiston",langevinPistonOn)
226  SCRIPT_PARSE_MOD_FLOAT("LangevinPistonTarget",
227  langevinPistonTarget,/PRESSUREFACTOR)
228  SCRIPT_PARSE_FLOAT("LangevinPistonPeriod",langevinPistonPeriod)
229  SCRIPT_PARSE_FLOAT("LangevinPistonDecay",langevinPistonDecay)
230  SCRIPT_PARSE_FLOAT("LangevinPistonTemp",langevinPistonTemp)
231  SCRIPT_PARSE_MOD_FLOAT("SurfaceTensionTarget",
232  surfaceTensionTarget,*(100.0/PRESSUREFACTOR))
233  SCRIPT_PARSE_BOOL("BerendsenPressure",berendsenPressureOn)
234  SCRIPT_PARSE_MOD_FLOAT("BerendsenPressureTarget",
235  berendsenPressureTarget,/PRESSUREFACTOR)
236  SCRIPT_PARSE_MOD_FLOAT("BerendsenPressureCompressibility",
237  berendsenPressureCompressibility,*PRESSUREFACTOR)
238  SCRIPT_PARSE_FLOAT("BerendsenPressureRelaxationTime",
239  berendsenPressureRelaxationTime)
240  SCRIPT_PARSE_FLOAT("constraintScaling",constraintScaling)
241  SCRIPT_PARSE_FLOAT("consForceScaling",consForceScaling)
242  SCRIPT_PARSE_BOOL("drudeHardWall",drudeHardWallOn)
243  SCRIPT_PARSE_FLOAT("drudeBondConst",drudeBondConst)
244  SCRIPT_PARSE_FLOAT("drudeBondLen",drudeBondLen)
245  SCRIPT_PARSE_STRING("outputname",outputFilename)
246  SCRIPT_PARSE_INT("outputEnergies",outputEnergies)
247  SCRIPT_PARSE_STRING("restartname",restartFilename)
248  SCRIPT_PARSE_INT("DCDfreq",dcdFrequency)
249  if ( ! strncasecmp(param,"DCDfile",MAX_SCRIPT_PARAM_SIZE) ) {
250  close_dcdfile(); // *** implemented in Output.C ***
251  strcpy(dcdFilename,value);
252  return;
253  }
254  if ( ! strncasecmp(param,"velDCDfile",MAX_SCRIPT_PARAM_SIZE) ) {
255  close_veldcdfile(); // *** implemented in Output.C ***
256  strcpy(velDcdFilename,value);
257  return;
258  }
259  SCRIPT_PARSE_STRING("tclBCArgs",tclBCArgs)
260  SCRIPT_PARSE_VECTOR("eField",eField)
261  SCRIPT_PARSE_FLOAT("eFieldFreq",eFieldFreq)
262  SCRIPT_PARSE_FLOAT("eFieldPhase",eFieldPhase)
263  SCRIPT_PARSE_FLOAT("accelMDE",accelMDE)
264  SCRIPT_PARSE_FLOAT("accelMDalpha",accelMDalpha)
265  SCRIPT_PARSE_FLOAT("accelMDTE",accelMDTE)
266  SCRIPT_PARSE_FLOAT("accelMDTalpha",accelMDTalpha)
267  SCRIPT_PARSE_FLOAT("accelMDGSigma0P",accelMDGSigma0P)
268  SCRIPT_PARSE_FLOAT("accelMDGSigma0D",accelMDGSigma0D)
269  SCRIPT_PARSE_STRING("accelMDGRestartFile",accelMDGRestartFile)
270  SCRIPT_PARSE_VECTOR("stirAxis",stirAxis)
271  SCRIPT_PARSE_VECTOR("stirPivot",stirPivot)
272  if ( ! strncasecmp(param,"mgridforcescale",MAX_SCRIPT_PARAM_SIZE) ) {
273  NAMD_die("Can't yet modify mgridforcescale in a script");
274  return;
275  }
276  if ( ! strncasecmp(param,"mgridforcevoff",MAX_SCRIPT_PARAM_SIZE) ) {
277  NAMD_die("Can't yet modify mgridforcevoff in a script");
278  return;
279  }
280 
281  if ( ! strncasecmp(param,"fixedatoms",MAX_SCRIPT_PARAM_SIZE) ) {
282  if ( ! fixedAtomsOn )
283  NAMD_die("FixedAtoms may not be enabled in a script.");
284  if ( ! fixedAtomsForces )
285  NAMD_die("To use fixedAtoms in script first use fixedAtomsForces yes.");
286  fixedAtomsOn = atobool(value);
287  return;
288  }
289 
290 //fepb
291  if ( ! strncasecmp(param,"alch",MAX_SCRIPT_PARAM_SIZE) ) {
292  alchOn = atobool(value);
293  if ( alchOn && ! alchOnAtStartup ) {
294  NAMD_die("Alchemy must be enabled at startup to disable and re-enable in script.");
295  }
296  alchFepOn = alchOn && alchFepOnAtStartup;
297  alchThermIntOn = alchOn && alchThermIntOnAtStartup;
299  if ( PMEOn ) ComputePmeUtil::select();
300  return;
301  }
302  SCRIPT_PARSE_INT("alchEquilSteps",alchEquilSteps)
303 
304  if ( ! strncasecmp(param,"alchLambda",MAX_SCRIPT_PARAM_SIZE) ) {
305  alchLambda = atof(value);
306  if ( alchLambda < 0.0 || 1.0 < alchLambda ) {
307  NAMD_die("Alchemical lambda values should be in the range [0.0, 1.0]\n");
308  }
310  return;
311  }
312 
313  if ( ! strncasecmp(param,"alchLambda2",MAX_SCRIPT_PARAM_SIZE) ) {
314  alchLambda2 = atof(value);
315  if ( alchLambda2 < 0.0 || 1.0 < alchLambda2 ) {
316  NAMD_die("Alchemical lambda values should be in the range [0.0, 1.0]\n");
317  }
319  return;
320  }
321 
322  if ( ! strncasecmp(param,"alchLambdaIDWS",MAX_SCRIPT_PARAM_SIZE) ) {
323  alchLambdaIDWS = atof(value);
324  setupIDWS();
326  return;
327  }
328 
329  if ( ! strncasecmp(param,"alchLambdaFreq",MAX_SCRIPT_PARAM_SIZE) ) {
330  alchLambdaFreq = atoi(value);
331  if ( alchLambdaIDWS >= 0 ) {
332  NAMD_die("alchLambdaIDWS and alchLambdaFreq are not compatible.\n");
333  }
335  return;
336  }
337 //fepe
338 
339  if ( ! strncasecmp(param,"nonbondedScaling",MAX_SCRIPT_PARAM_SIZE) ) {
340  nonbondedScaling = atof(value);
342  return;
343  }
344 
345  if ( ! strncasecmp(param,"commOnly",MAX_SCRIPT_PARAM_SIZE) ) {
346  commOnly = atobool(value);
348  return;
349  }
350 
351  // REST2 - We don't have to make any changes to ComputeNonbondedUtil
352  // other than recalculating the LJTable. Skip doing the other stuff.
353  if ( ! strncasecmp(param,"soluteScalingFactor",MAX_SCRIPT_PARAM_SIZE)) {
354  soluteScalingFactor = atof(value);
355  if (soluteScalingFactor < 0.0) {
356  NAMD_die("Solute scaling factor should be non-negative\n");
357  }
358  soluteScalingFactorCharge = soluteScalingFactor;
359  soluteScalingFactorVdw = soluteScalingFactor;
360  // update LJTable for CPU
362 #if defined(NAMD_CUDA) || defined(NAMD_HIP)
363  // update LJTable for GPU, needs CPU update first
365 #endif
366  return;
367  }
368  if ( ! strncasecmp(param,"soluteScalingFactorVdw",MAX_SCRIPT_PARAM_SIZE)) {
369  soluteScalingFactorVdw = atof(value);
370  if (soluteScalingFactorVdw < 0.0) {
371  NAMD_die("Solute scaling factor for van der Waals "
372  "should be non-negative\n");
373  }
374  // update LJTable for CPU
376 #if defined(NAMD_CUDA) || defined(NAMD_HIP)
377  // update LJTable for GPU, needs CPU update first
379 #endif
380  return;
381  }
382  if ( ! strncasecmp(param,"soluteScalingFactorCharge",MAX_SCRIPT_PARAM_SIZE)) {
383  soluteScalingFactorCharge = atof(value);
384  if (soluteScalingFactorCharge < 0.0) {
385  NAMD_die("Solute scaling factor for electrostatics "
386  "should be non-negative\n");
387  }
388  return;
389  }
390 
391  char *error = new char[2 * MAX_SCRIPT_PARAM_SIZE + 100];
392  sprintf(error,"Setting parameter %s from script failed!\n",param);
393  NAMD_die(error);
394 
395 }
396 
399 }
400 
403 }
404 
405 /************************************************************************/
406 /* */
407 /* FUNCTION config_parser */
408 /* */
409 /************************************************************************/
410 
411 void SimParameters::config_parser(ParseOptions &opts) {
412 
413  // Set all variable to fallback default values. This is not really
414  // necessary, as we give default values when we set up the ParseOptions
415  // object, but it helps the debuggers figure out we've initialized the
416  // variables.
417  HydrogenBonds = FALSE;
418  useAntecedent = TRUE;
419  aaAngleExp = 2;
420  haAngleExp = 4;
421  distAttExp = 4;
422  distRepExp = 6;
423  dhaCutoffAngle = 100.0;
424  dhaOnAngle = 60.0;
425  dhaOffAngle = 80.0;
426  daCutoffDist = 7.5;
427  daOnDist = 5.5;
428  daOffDist = 6.5;
429 
430  config_parser_basic(opts);
431  config_parser_fileio(opts);
432  config_parser_fullelect(opts);
433  config_parser_methods(opts);
434  config_parser_constraints(opts);
435  #ifdef OPENATOM_VERSION
436  config_parser_openatom(opts);
437  #endif // OPENATOM_VERSION
438 
439  config_parser_gridforce(opts);
440  config_parser_mgridforce(opts);
441  config_parser_movdrag(opts);
442  config_parser_rotdrag(opts);
443  config_parser_constorque(opts);
444  config_parser_boundary(opts);
445  config_parser_misc(opts);
446 
447 }
448 
449 void SimParameters::config_parser_basic(ParseOptions &opts) {
450 
451  // So first we set up the ParseOptions objects so that it will check
452  // all of the logical rules that the configuration file must follow.
453 
454  opts.optional("main", "obsolete", "used to flag obsolete options",
455  PARSE_STRING);
456 
458  opts.require("main", "timestep", "size of the timestep, in fs",
459  &dt, 1.0);
460  opts.range("timestep", NOT_NEGATIVE);
461  opts.units("timestep", N_FSEC);
462 
463  opts.optional("main", "numsteps", "number of timesteps to perform",
464  &N,0);
465  opts.range("numsteps", NOT_NEGATIVE);
466 
467  opts.optional("main", "stepspercycle",
468  "Number of steps between atom migrations",
469  &stepsPerCycle, 20);
470  opts.range("stepspercycle", POSITIVE);
471 
472  opts.require("main", "cutoff", "local electrostatic and Vdw distance",
473  &cutoff);
474  opts.range("cutoff", POSITIVE);
475  opts.units("cutoff", N_ANGSTROM);
476 
477  opts.optional("main", "nonbondedScaling", "nonbonded scaling factor",
478  &nonbondedScaling, 1.0);
479  opts.range("nonbondedScaling", NOT_NEGATIVE);
480 
481  opts.optional("main", "limitDist", "limit nonbonded below this distance",
482  &limitDist, 0.0);
483  opts.range("limitDist", NOT_NEGATIVE);
484 
485  opts.require("main", "exclude", "Electrostatic and VDW exclusion policy",
486  PARSE_STRING);
487 
488  opts.optional("exclude", "1-4scaling", "1-4 electrostatic scaling factor",
489  &scale14, 1.0);
490  opts.range("1-4scaling", POSITIVE);
491 
492  opts.optionalB("main", "switching",
493  "Should a smoothing function be used?", &switchingActive, TRUE);
494 
495  opts.optionalB("switching", "vdwForceSwitching",
496  "Use force switching for vdw?", &vdwForceSwitching, FALSE);
497 
498  opts.optional("switching", "switchdist",
499  "Distance for switching function activation",
500  &switchingDist);
501  opts.range("switchdist", POSITIVE);
502  opts.units("switchdist", N_ANGSTROM);
503 
504  opts.optionalB("main", "martiniSwitching",
505  "Use Martini residue-based coarse-grain switching?", &martiniSwitching, FALSE);
506  opts.optionalB("main", "martiniDielAllow",
507  "Allow use of dielectric != 15.0 when using Martini", &martiniDielAllow, FALSE);
508 
509  opts.optional("main", "pairlistdist", "Pairlist inclusion distance",
510  &pairlistDist);
511  opts.range("pairlistdist", POSITIVE);
512  opts.units("pairlistdist", N_ANGSTROM);
513 
514  opts.optional("main", "pairlistMinProcs", "Min procs for pairlists",
515  &pairlistMinProcs,1);
516  opts.range("pairlistMinProcs", POSITIVE);
517 
518  opts.optional("main", "pairlistsPerCycle", "regenerate x times per cycle",
519  &pairlistsPerCycle,2);
520  opts.range("pairlistsPerCycle", POSITIVE);
521 
522  opts.optional("main", "outputPairlists", "how often to print warnings",
523  &outputPairlists, 0);
524  opts.range("outputPairlists", NOT_NEGATIVE);
525 
526  opts.optional("main", "pairlistShrink", "tol *= (1 - x) on regeneration",
527  &pairlistShrink,0.01);
528  opts.range("pairlistShrink", NOT_NEGATIVE);
529 
530  opts.optional("main", "pairlistGrow", "tol *= (1 + x) on trigger",
531  &pairlistGrow, 0.01);
532  opts.range("pairlistGrow", NOT_NEGATIVE);
533 
534  opts.optional("main", "pairlistTrigger", "trigger is atom > (1 - x) * tol",
535  &pairlistTrigger, 0.3);
536  opts.range("pairlistTrigger", NOT_NEGATIVE);
537 
538  opts.optional("main", "temperature", "initial temperature",
539  &initialTemp);
540  opts.range("temperature", NOT_NEGATIVE);
541  opts.units("temperature", N_KELVIN);
542 
543  opts.optionalB("main", "COMmotion", "allow initial center of mass movement",
544  &comMove, FALSE);
545 
546  opts.optionalB("main", "zeroMomentum", "constrain center of mass",
547  &zeroMomentum, FALSE);
548  opts.optionalB("zeroMomentum", "zeroMomentumAlt", "constrain center of mass",
549  &zeroMomentumAlt, FALSE);
550 
551  opts.optionalB("main", "wrapWater", "wrap waters around periodic boundaries on output",
552  &wrapWater, FALSE);
553  opts.optionalB("main", "wrapAll", "wrap all clusters around periodic boundaries on output",
554  &wrapAll, FALSE);
555  opts.optionalB("main", "wrapNearest", "wrap to nearest image to cell origin",
556  &wrapNearest, FALSE);
557 
558  opts.optional("main", "dielectric", "dielectric constant",
559  &dielectric, 1.0);
560  opts.range("dielectric", POSITIVE); // Hmmm, dielectric < 1 ...
561 
562  opts.optional("main", "margin", "Patch width margin", &margin, XXXBIGREAL);
563  opts.range("margin", NOT_NEGATIVE);
564  opts.units("margin", N_ANGSTROM);
565 
566  opts.optional("main", "seed", "Initial random number seed", &randomSeed);
567  opts.range("seed", POSITIVE);
568 
569  opts.optional("main", "outputEnergies", "How often to print energies in timesteps",
570  &outputEnergies, 1);
571  opts.range("outputEnergies", POSITIVE);
572 
573  opts.optional("main", "outputMomenta", "How often to print linear and angular momenta in timesteps",
574  &outputMomenta, 0);
575  opts.range("outputMomenta", NOT_NEGATIVE);
576 
577  opts.optional("main", "outputTiming", "How often to print timing data in timesteps",
578  &outputTiming);
579  opts.range("outputTiming", NOT_NEGATIVE);
580 
581  opts.optional("main", "outputCudaTiming", "How often to print CUDA timing data in timesteps",
582  &outputCudaTiming, 0);
583  opts.range("outputCudaTiming", NOT_NEGATIVE);
584 
585  opts.optional("main", "outputPressure", "How often to print pressure data in timesteps",
586  &outputPressure, 0);
587  opts.range("outputPressure", NOT_NEGATIVE);
588 
589  opts.optionalB("main", "mergeCrossterms", "merge crossterm energy with dihedral when printing?",
590  &mergeCrossterms, TRUE);
591 
592  opts.optional("main", "MTSAlgorithm", "Multiple timestep algorithm",
593  PARSE_STRING);
594 
595  opts.optional("main", "longSplitting", "Long range force splitting option",
596  PARSE_STRING);
597 
598  opts.optionalB("main", "ignoreMass", "Do not use masses to find hydrogen atoms",
599  &ignoreMass, FALSE);
600 
601  opts.optional("main", "splitPatch", "Atom into patch splitting option",
602  PARSE_STRING);
603  opts.optional("main", "hgroupCutoff", "Hydrogen margin", &hgroupCutoff, 2.5);
604 
605  opts.optional("main", "extendedSystem",
606  "Initial configuration of extended system variables and periodic cell",
607  PARSE_STRING);
608 
609  opts.optional("main", "cellBasisVector1", "Basis vector for periodic cell",
610  &cellBasisVector1);
611  opts.optional("main", "cellBasisVector2", "Basis vector for periodic cell",
612  &cellBasisVector2);
613  opts.optional("main", "cellBasisVector3", "Basis vector for periodic cell",
614  &cellBasisVector3);
615  opts.optional("main", "cellOrigin", "Fixed center of periodic cell",
616  &cellOrigin);
617 
618  opts.optionalB("main", "molly", "Rigid bonds to hydrogen",&mollyOn,FALSE);
619  opts.optional("main", "mollyTolerance", "Error tolerance for MOLLY",
620  &mollyTol, 0.00001);
621  opts.optional("main", "mollyIterations",
622  "Max number of iterations for MOLLY", &mollyIter, 100);
623 
624  opts.optional("main", "rigidBonds", "Rigid bonds to hydrogen",PARSE_STRING);
625  opts.optional("main", "rigidTolerance",
626  "Error tolerance for rigid bonds to hydrogen",
627  &rigidTol, 1.0e-8);
628  opts.optional("main", "rigidIterations",
629  "Max number of SHAKE iterations for rigid bonds to hydrogen",
630  &rigidIter, 100);
631  opts.optionalB("main", "rigidDieOnError",
632  "Die if rigidTolerance is not achieved after rigidIterations",
633  &rigidDie, TRUE);
634  opts.optionalB("main", "useSettle",
635  "Use the SETTLE algorithm for rigid waters",
636  &useSettle, TRUE);
637 
638  opts.optional("main", "nonbondedFreq", "Nonbonded evaluation frequency",
639  &nonbondedFrequency, 1);
640  opts.range("nonbondedFreq", POSITIVE);
641 
642  opts.optionalB("main", "outputPatchDetails", "print number of atoms in each patch",
643  &outputPatchDetails, FALSE);
644  opts.optionalB("main", "staticAtomAssignment", "never migrate atoms",
645  &staticAtomAssignment, FALSE);
646  opts.optionalB("main", "replicaUniformPatchGrids", "same patch grid size on all replicas",
647  &replicaUniformPatchGrids, FALSE);
648 #ifndef MEM_OPT_VERSION
649  // in standard (non-mem-opt) version, enable lone pairs by default
650  // for compatibility with recent force fields
651  opts.optionalB("main", "lonePairs", "Enable lone pairs", &lonepairs, TRUE);
652 #else
653  // in mem-opt version, disable lone pairs by default
654  // because they are not supported
655  opts.optionalB("main", "lonePairs", "Enable lone pairs", &lonepairs, FALSE);
656 #endif
657  opts.optional("main", "waterModel", "Water model to use", PARSE_STRING);
658  opts.optionalB("main", "LJcorrection", "Apply analytical tail corrections for energy and virial", &LJcorrection, FALSE);
659 #ifdef TIMER_COLLECTION
660  opts.optional("main", "TimerBinWidth",
661  "Bin width of timer histogram collection in microseconds",
662  &timerBinWidth, 1.0);
663 #endif
664 #if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED) || defined(NAMD_ROCTX_ENABLED)
665  // default NVTX or Projections profiling is up to the first 1000 patches
666  opts.optional("main", "beginEventPatchID","Beginning patch ID for profiling",
667  &beginEventPatchID, 0);
668  opts.optional("main", "endEventPatchID", "Ending patch ID for profiling",
669  &endEventPatchID, 5000);
670  // default NVTX or Projections profiling is up to the first 1000 time steps
671  opts.optional("main", "beginEventStep", "Beginning time step for profiling",
672  &beginEventStep, 0);
673  opts.optional("main", "endEventStep", "Ending time step for profiling",
674  &endEventStep, 1000);
675 #endif
676 }
677 
678 void SimParameters::config_parser_fileio(ParseOptions &opts) {
679 
681 
682  opts.optional("main", "cwd", "current working directory", PARSE_STRING);
683 
684 // In order to include AMBER options, "coordinates", "structure"
685 // and "parameters" are now optional, not required. The presence
686 // of them will be checked later in check_config()
687 
688 // opts.require("main", "coordinates", "initial PDB coordinate file",
689 // PARSE_STRING);
690  opts.optional("main", "coordinates", "initial PDB coordinate file",
691  PARSE_STRING);
692 
693  opts.optional("main", "velocities",
694  "initial velocities, given as a PDB file", PARSE_STRING);
695  opts.optional("main", "binvelocities",
696  "initial velocities, given as a binary restart", PARSE_STRING);
697  opts.optional("main", "bincoordinates",
698  "initial coordinates in a binary restart file", PARSE_STRING);
699 #ifdef MEM_OPT_VERSION
700  opts.optional("main", "binrefcoords",
701  "reference coordinates in a binary restart file", PARSE_STRING);
702 #endif
703 
704 // opts.require("main", "structure", "initial PSF structure file",
705 // PARSE_STRING);
706  opts.optional("main", "structure", "initial PSF structure file",
707  PARSE_STRING);
708 
709 // opts.require("main", "parameters",
710 //"CHARMm 19 or CHARMm 22 compatable force field file (multiple "
711 //"inputs allowed)", PARSE_MULTIPLES);
712  opts.optional("main", "parameters",
713 "CHARMm 19 or CHARMm 22 compatable force field file (multiple "
714 "inputs allowed)", PARSE_MULTIPLES);
715 
716 
717  //****** BEGIN CHARMM/XPLOR type changes
719  opts.optionalB("parameters", "paraTypeXplor", "Parameter file in Xplor format?", &paraTypeXplorOn, FALSE);
720  opts.optionalB("parameters", "paraTypeCharmm", "Parameter file in Charmm format?", &paraTypeCharmmOn, FALSE);
721  //****** END CHARMM/XPLOR type changes
722 
723  // Ported by JLai -- JE - Go parameters
724  opts.optionalB("main", "GromacsPair", "Separately calculate pair interactions", &goGroPair, FALSE);
725  opts.optionalB("main", "GoForcesOn", "Go forces will be calculated", &goForcesOn, FALSE);
726  opts.require("GoForcesOn", "GoParameters", "Go parameter file", goParameters);
727  opts.require("GoForcesOn", "GoCoordinates", "target coordinates for Go forces", goCoordinates);
728  // Added by JLai -- Go-method parameter -- switches between using the matrix and sparse matrix representations [6.3.11]
729  // opts.optional("GoForcesOn", "GoMethod", "Which type of matrix should be used to store Go contacts?", &goMethod);
730  // Added by JLai -- Go-method parameter [8.2.11]
731  //opts.optional("GoForcesOn", "GoMethod", "Which type of matrix should be used to store Go contacts?", PARSE_STRING);
732  opts.require("GoForcesOn", "GoMethod", "Which type of matrix should be used to store Go contacts?", PARSE_STRING);
733  // End of Port -- JL
734 
735  opts.require("main", "outputname",
736  "prefix for the final PDB position and velocity filenames",
737  outputFilename);
738 
739  opts.optional("main", "auxFile", "Filename for data stream output",
740  auxFilename);
741 
742  opts.optional("main", "numinputprocs", "Number of pes to use for parallel input",
743  &numinputprocs, 0);
744  opts.range("numinputprocs", NOT_NEGATIVE);
745 
746  opts.optional("main", "numoutputprocs", "Number of pes to use for parallel output",
747  &numoutputprocs, 0);
748  opts.range("numoutputprocs", NOT_NEGATIVE);
749  opts.optional("main", "numoutputwriters", "Number of output processors that simultaneously write to an output file",
750  &numoutputwrts, 1);
751  opts.range("numoutputwriters", NOT_NEGATIVE);
752 
753  opts.optional("main", "DCDfreq", "Frequency of DCD trajectory output, in "
754  "timesteps", &dcdFrequency, 0);
755  opts.range("DCDfreq", NOT_NEGATIVE);
756  opts.optional("DCDfreq", "DCDfile", "DCD trajectory output file name",
757  dcdFilename);
758  opts.optionalB("DCDfreq", "DCDunitcell", "Store unit cell in dcd timesteps?",
759  &dcdUnitCell);
760 
761  opts.optional("main", "velDCDfreq", "Frequency of velocity "
762  "DCD output, in timesteps", &velDcdFrequency, 0);
763  opts.range("velDCDfreq", NOT_NEGATIVE);
764  opts.optional("velDCDfreq", "velDCDfile", "velocity DCD output file name",
765  velDcdFilename);
766 
767  opts.optional("main", "forceDCDfreq", "Frequency of force"
768  "DCD output, in timesteps", &forceDcdFrequency, 0);
769  opts.range("forceDCDfreq", NOT_NEGATIVE);
770  opts.optional("forceDCDfreq", "forceDCDfile", "force DCD output file name",
771  forceDcdFilename);
772 
773  opts.optional("main", "XSTfreq", "Frequency of XST trajectory output, in "
774  "timesteps", &xstFrequency, 0);
775  opts.range("XSTfreq", NOT_NEGATIVE);
776  opts.optional("XSTfreq", "XSTfile", "Extended sytem trajectory output "
777  "file name", xstFilename);
778 
779  opts.optional("main", "restartfreq", "Frequency of restart file "
780  "generation", &restartFrequency, 0);
781  opts.range("restartfreq", NOT_NEGATIVE);
782  opts.optional("restartfreq", "restartname", "Prefix for the position and "
783  "velocity PDB files used for restarting", restartFilename);
784  opts.optionalB("restartfreq", "restartsave", "Save restart files with "
785  "unique filenames rather than overwriting", &restartSave, FALSE);
786  opts.optionalB("restartfreq", "restartsavedcd", "Save DCD files with "
787  "unique filenames at each restart", &restartSaveDcd, FALSE);
788 
789  opts.optionalB("restartfreq", "binaryrestart", "Specify use of binary restart files ",
790  &binaryRestart, TRUE);
791 
792  opts.optionalB("outputname", "binaryoutput", "Specify use of binary output files ",
793  &binaryOutput, TRUE);
794 
795  opts.optionalB("main", "amber", "Is it AMBER force field?",
796  &amberOn, FALSE);
797  opts.optionalB("amber", "readexclusions", "Read exclusions from parm file?",
798  &readExclusions, TRUE);
799  opts.require("amber", "scnb", "1-4 VDW interactions are divided by scnb",
800  &vdwscale14, 2.0);
801  opts.require("amber", "parmfile", "AMBER parm file", PARSE_STRING);
802  opts.optional("amber", "ambercoor", "AMBER coordinate file", PARSE_STRING);
803 
804  /* GROMACS options */
805  opts.optionalB("main", "gromacs", "Use GROMACS-like force field?",
806  &gromacsOn, FALSE);
807  opts.require("gromacs", "grotopfile", "GROMACS topology file",
808  PARSE_STRING);
809  opts.optional("gromacs", "grocoorfile","GROMACS coordinate file",
810  PARSE_STRING);
811 
812  // OPLS options
813  opts.optionalB("main", "vdwGeometricSigma",
814  "Use geometric mean to combine L-J sigmas, as for OPLS",
815  &vdwGeometricSigma, FALSE);
816 
817  // load/store computeMap
818  opts.optional("main", "computeMapFile", "Filename for computeMap",
819  computeMapFilename);
820  opts.optionalB("main", "storeComputeMap", "store computeMap?",
821  &storeComputeMap, FALSE);
822  opts.optionalB("main", "loadComputeMap", "load computeMap?",
823  &loadComputeMap, FALSE);
824 }
825 
826 
827 void SimParameters::config_parser_fullelect(ParseOptions &opts) {
828 
830 #ifdef DPMTA
831  DebugM(1,"DPMTA setup start\n");
832  // PMTA is included, so really get these values
833  opts.optionalB("main", "FMA", "Should FMA be used?", &FMAOn, FALSE);
834  opts.optional("FMA", "FMALevels", "Tree levels to use in FMA", &FMALevels,
835  5);
836  opts.range("FMALevels", POSITIVE);
837  opts.optional("FMA", "FMAMp", "Number of FMA multipoles", &FMAMp, 8);
838  opts.range("FMAMp", POSITIVE);
839  opts.optionalB("FMA", "FMAFFT", "Use FFT enhancement in FMA?", &FMAFFTOn, TRUE);
840  opts.optional("FMAFFT", "FMAFFTBlock", "FFT blocking factor",
841  &FMAFFTBlock, 4);
842  opts.range("FMAFFTBlock", POSITIVE);
843  DebugM(1,"DPMTA setup end\n");
844 #else
845  // PMTA is NOT included. So just set all the values to 0.
846  FMAOn = FALSE;
847  FMALevels = 0;
848  FMAMp = 0;
849  FMAFFTOn = FALSE;
850  FMAFFTBlock = 0;
851 #endif
852 
853  opts.optional("main", "fullElectFrequency",
854  "Number of steps between full electrostatic executions",
855  &fullElectFrequency);
856  opts.range("fullElectFrequency", POSITIVE);
857 
858  // USE OF THIS PARAMETER DISCOURAGED
859  opts.optional("main", "fmaFrequency",
860  "Number of steps between full electrostatic executions",
861  &fmaFrequency);
862  opts.range("fmaFrequency", POSITIVE);
863 
864  opts.optional("main", "fmaTheta",
865  "FMA theta parameter value",
866  &fmaTheta,0.715);
867  opts.range("fmaTheta", POSITIVE);
868 
869  opts.optionalB("main", "FullDirect", "Should direct calculations of full electrostatics be performed?",
870  &fullDirectOn, FALSE);
871 
872 
874 
875  opts.optionalB("main", "MSM",
876  "Use multilevel summation method for electrostatics?",
877  &MSMOn, FALSE);
878  opts.optional("MSM", "MSMQuality", "MSM quality",
879  &MSMQuality, 0);
880  opts.optional("MSM", "MSMApprox", "MSM approximation",
881  &MSMApprox, 0);
882  opts.optional("MSM", "MSMSplit", "MSM splitting",
883  &MSMSplit, 0);
884  opts.optional("MSM", "MSMLevels", "MSM maximum number of levels",
885  &MSMLevels, 0); // set to 0 adapts to as many as needed
886  opts.optional("MSM", "MSMGridSpacing", "MSM grid spacing (Angstroms)",
887  &MSMGridSpacing, 2.5);
888  opts.optional("MSM", "MSMPadding", "MSM padding (Angstroms)",
889  &MSMPadding, 2.5);
890  opts.optional("MSM", "MSMxmin", "MSM x minimum (Angstroms)", &MSMxmin, 0);
891  opts.optional("MSM", "MSMxmax", "MSM x maximum (Angstroms)", &MSMxmax, 0);
892  opts.optional("MSM", "MSMymin", "MSM y minimum (Angstroms)", &MSMymin, 0);
893  opts.optional("MSM", "MSMymax", "MSM y maximum (Angstroms)", &MSMymax, 0);
894  opts.optional("MSM", "MSMzmin", "MSM z minimum (Angstroms)", &MSMzmin, 0);
895  opts.optional("MSM", "MSMzmax", "MSM z maximum (Angstroms)", &MSMzmax, 0);
896  opts.optional("MSM", "MSMBlockSizeX",
897  "MSM grid block size along X direction (for decomposing parallel work)",
898  &MSMBlockSizeX, 8);
899  opts.optional("MSM", "MSMBlockSizeY",
900  "MSM grid block size along Y direction (for decomposing parallel work)",
901  &MSMBlockSizeY, 8);
902  opts.optional("MSM", "MSMBlockSizeZ",
903  "MSM grid block size along Z direction (for decomposing parallel work)",
904  &MSMBlockSizeZ, 8);
905 
906  opts.optionalB("MSM", "MsmSerial",
907  "Use MSM serial version for long-range calculation?",
908  &MsmSerialOn, FALSE);
909 
910 
912 
913  opts.optionalB("main", "FMM",
914  "Use fast multipole method for electrostatics?",
915  &FMMOn, FALSE);
916  opts.optional("FMM", "FMMLevels", "FMM number of levels",
917  &FMMLevels, 0);
918  opts.optional("FMM", "FMMPadding", "FMM padding margin (Angstroms)",
919  &FMMPadding, 0);
920 
921  opts.optionalB("main", "useCUDA2", "Use new CUDA code", &useCUDA2, TRUE);
922 
924 
925  opts.optionalB("main", "PME", "Use particle mesh Ewald for electrostatics?",
926  &PMEOn, FALSE);
927  opts.optional("PME", "PMETolerance", "PME direct space tolerance",
928  &PMETolerance, 1.e-6);
929  opts.optional("PME", "PMEInterpOrder", "PME interpolation order",
930  &PMEInterpOrder, 4); // cubic interpolation is default
931  opts.optional("PME", "PMEGridSizeX", "PME grid in x dimension",
932  &PMEGridSizeX, 0);
933  opts.optional("PME", "PMEGridSizeY", "PME grid in y dimension",
934  &PMEGridSizeY, 0);
935  opts.optional("PME", "PMEGridSizeZ", "PME grid in z dimension",
936  &PMEGridSizeZ, 0);
937  opts.optional("PME", "PMEGridSpacing", "Maximum PME grid spacing (Angstroms)",
938  &PMEGridSpacing, 0.);
939  opts.range("PMEGridSpacing", NOT_NEGATIVE);
940  opts.optional("PME", "PMEProcessors",
941  "PME FFT and reciprocal sum processor count", &PMEProcessors, 0);
942  opts.optional("PME", "PMEMinSlices",
943  "minimum thickness of PME reciprocal sum slab", &PMEMinSlices, 2);
944  opts.range("PMEMinSlices", NOT_NEGATIVE);
945  opts.optional("PME", "PMEPencils",
946  "PME FFT and reciprocal sum pencil grid size", &PMEPencils, -1);
947  opts.optional("PME", "PMEPencilsX",
948  "PME FFT and reciprocal sum pencil grid size X", &PMEPencilsX, 0);
949  opts.optional("PME", "PMEPencilsY",
950  "PME FFT and reciprocal sum pencil grid size Y", &PMEPencilsY, 0);
951  opts.optional("PME", "PMEPencilsZ",
952  "PME FFT and reciprocal sum pencil grid size Z", &PMEPencilsZ, 0);
953  opts.range("PMEPencilsX", NOT_NEGATIVE);
954  opts.range("PMEPencilsY", NOT_NEGATIVE);
955  opts.range("PMEPencilsZ", NOT_NEGATIVE);
956  opts.optional("PME", "PMEPencilsYLayout",
957  "PME FFT and reciprocal sum Y pencil layout strategy", &PMEPencilsYLayout, 0);
958  opts.optional("PME", "PMEPencilsXLayout",
959  "PME FFT and reciprocal sum X pencil layout strategy", &PMEPencilsXLayout, 1);
960  opts.range("PMEPencilsYLayout", NOT_NEGATIVE);
961  opts.range("PMEPencilsXLayout", NOT_NEGATIVE);
962  opts.optional("PME", "PMESendOrder",
963  "PME message ordering control", &PMESendOrder, 0);
964  opts.range("PMESendOrder", NOT_NEGATIVE);
965  opts.optional("PME", "PMEMinPoints",
966  "minimum points per PME reciprocal sum pencil", &PMEMinPoints, 10000);
967  opts.range("PMEMinPoints", NOT_NEGATIVE);
968  opts.optionalB("main", "PMEBarrier", "Use barrier in PME?",
969  &PMEBarrier, FALSE);
970  opts.optionalB("main", "PMEOffload", "Offload PME to accelerator?",
971  &PMEOffload);
972 
973  opts.optionalB("PME", "usePMECUDA", "Use the PME CUDA version", &usePMECUDA, CmiNumPhysicalNodes() < 5);
974 
975 #ifdef DPME
976  opts.optionalB("PME", "useDPME", "Use old DPME code?", &useDPME, FALSE);
977 #else
978  useDPME = 0;
979 #endif
980  opts.optionalB("main", "FFTWPatient", "Use intensive plan creation to optimize FFTW?",
981 #ifdef WIN32
982  &FFTWPatient, TRUE);
983 #else
984  &FFTWPatient, FALSE);
985 #endif
986 
987  opts.optionalB("main", "FFTWEstimate", "Use estimates to optimize FFTW?",
988 #ifdef WIN32
989  &FFTWEstimate, TRUE);
990 #else
991  &FFTWEstimate, FALSE);
992 #endif
993  opts.optionalB("main", "FFTWUseWisdom", "Read/save wisdom file for FFTW?",
994 #ifdef WIN32
995  &FFTWUseWisdom, FALSE);
996 #else
997  &FFTWUseWisdom, TRUE);
998 #endif
999  opts.optional("FFTWUseWisdom", "FFTWWisdomFile", "File for FFTW wisdom",
1000  FFTWWisdomFile);
1001 
1002  opts.optionalB("main", "useAVXTiles",
1003  "Use \"Tiles\" mode for AVX-512 optimized calculations",
1004  &useAVXTiles, TRUE);
1005 }
1006 
1007 void SimParameters::config_parser_methods(ParseOptions &opts) {
1008 
1010  opts.optionalB("main", "minimization", "Should minimization be performed?",
1011  &minimizeCGOn, FALSE);
1012  opts.optionalB("main", "minVerbose", "Print extra minimization diagnostics?",
1013  &minVerbose, FALSE);
1014  opts.optional("main", "minTinyStep", "very first minimization steps",
1015  &minTinyStep, 1.0e-6);
1016  opts.range("minTinyStep", POSITIVE);
1017  opts.optional("main", "minBabyStep", "initial minimization steps",
1018  &minBabyStep, 1.0e-2);
1019  opts.range("minBabyStep", POSITIVE);
1020  opts.optional("main", "minLineGoal", "line minimization gradient reduction",
1021  &minLineGoal, 1.0e-3);
1022  opts.range("minLineGoal", POSITIVE);
1023 
1024  opts.optionalB("main", "velocityQuenching",
1025  "Should old-style minimization be performed?", &minimizeOn, FALSE);
1026 
1027  opts.optional("main", "maximumMove", "Maximum atom movement per step", &maximumMove, 0.0);
1028  opts.range("maximumMove", NOT_NEGATIVE);
1029  opts.units("maximumMove", N_ANGSTROM);
1030 
1031  opts.optionalB("main", "Langevin", "Should Langevin dynamics be performed?",
1032  &langevinOn, FALSE);
1033  opts.require("Langevin", "langevinTemp", "Temperature for heat bath in Langevin "
1034  "dynamics", &langevinTemp);
1035  opts.range("langevinTemp", NOT_NEGATIVE);
1036  opts.units("langevinTemp", N_KELVIN);
1037  opts.optional("Langevin", "langevinDamping", "Damping coefficient (1/ps)",
1038  &langevinDamping);
1039  opts.range("langevinDamping", POSITIVE);
1040  opts.optionalB("Langevin", "langevinHydrogen", "Should Langevin dynamics be applied to hydrogen atoms?",
1041  &langevinHydrogen);
1042  opts.optional("Langevin", "langevinFile", "PDB file with temperature "
1043  "coupling terms (B(i)) (default is the PDB input file)",
1044  PARSE_STRING);
1045  opts.optional("Langevin", "langevinCol", "Column in the langevinFile "
1046  "containing the temperature coupling term B(i);\n"
1047  "default is 'O'", PARSE_STRING);
1048 
1049  // use BAOAB integration instead of BBK
1050  opts.optionalB("Langevin", "langevinBAOAB",
1051  "Should Langevin dynamics be performed using BAOAB integration?",
1052  &langevin_useBAOAB, FALSE);
1053 
1054 // BEGIN LA
1055  opts.optionalB("main", "LoweAndersen", "Should Lowe-Andersen dynamics be performed?",
1056  &loweAndersenOn, FALSE);
1057  opts.require("LoweAndersen", "loweAndersenTemp", "Temperature for heat bath in Lowe-Andersen "
1058  "dynamics", &loweAndersenTemp);
1059  opts.range("loweAndersenTemp", NOT_NEGATIVE);
1060  opts.units("loweAndersenTemp", N_KELVIN);
1061  opts.optional("LoweAndersen", "loweAndersenRate", "Collision rate (1/ps)",
1062  &loweAndersenRate, 50);
1063  opts.range("loweAndersenRate", POSITIVE);
1064  opts.optional("LoweAndersen", "loweAndersenCutoff", "Cutoff radius",
1065  &loweAndersenCutoff, 2.7);
1066  opts.range("loweAndersenCutoff", POSITIVE);
1067  opts.units("loweAndersenCutoff", N_ANGSTROM);
1068 // END LA
1069 
1070 //fepb
1071  opts.optionalB("main", "alch", "Is achemical simulation being performed?",
1072  &alchOn, FALSE);
1073  opts.require("alch", "alchLambda", "Alchemical coupling parameter value",
1074  &alchLambda);
1075  opts.range("alchLambda", NOT_NEGATIVE);
1076 
1077  opts.optionalB("alch", "singleTopology",
1078  "Is single topology used for relative free energy?", &singleTopology, FALSE);
1079 
1080  opts.optionalB("alch", "sdBondScaling",
1081  "Is S-D bonded terms scaling for relative free energy?", &sdScaling, FALSE);
1082 
1083  opts.optional("alch", "alchFile", "PDB file with perturbation flags "
1084  "default is the input PDB file", PARSE_STRING);
1085  opts.optional("alch", "alchCol", "Column in the alchFile with the "
1086  "perturbation flag", PARSE_STRING);
1087 
1088  opts.optional("alch", "unperturbedBondFile", "mini psf file with unperturbed bond info"
1089  " ", PARSE_STRING);
1090  opts.optional("alch", "alchOutFreq", "Frequency of alchemical energy"
1091  "output in timesteps", &alchOutFreq, 5);
1092  opts.range("alchoutfreq", NOT_NEGATIVE);
1093  opts.optional("alch", "alchOutFile", "Alchemical energy output filename",
1094  alchOutFile);
1095 
1096  // soft-core parameters
1097  opts.optional("alch", "alchVdwShiftCoeff", "Coeff used for generating"
1098  "the altered alchemical vDW interactions", &alchVdwShiftCoeff, 5.);
1099  opts.range("alchVdwShiftCoeff", NOT_NEGATIVE);
1100 
1101  opts.optionalB("alch", "alchWCA", "Is WCA decomposition being performed?",
1102  &alchWCAOn, FALSE);
1103 
1104  // scheduling options for different interaction types
1105  opts.optional("alch", "alchElecLambdaStart", "Lambda at which electrostatic"
1106  "scaling of exnihilated particles begins", &alchElecLambdaStart, 0.5);
1107  opts.range("alchElecLambdaStart", NOT_NEGATIVE);
1108 
1109  opts.optional("alch", "alchVdwLambdaEnd", "Lambda at which vdW"
1110  "scaling of exnihilated particles ends", &alchVdwLambdaEnd, 1.0);
1111  opts.range("alchVdwLambdaEnd", NOT_NEGATIVE);
1112 
1113  opts.optional("alch", "alchRepLambdaEnd", "Lambda at which repulsive vdW"
1114  "scaling of exnihilated particles ends and attractive vdW scaling"
1115  "begins", &alchRepLambdaEnd, 0.5);
1116  opts.range("alchRepLambdaEnd", NOT_NEGATIVE);
1117 
1118  opts.optional("alch", "alchBondLambdaEnd", "Lambda at which bonded"
1119  "scaling of exnihilated particles begins", &alchBondLambdaEnd, 0.0);
1120  opts.range("alchBondLambdaEnd", NOT_NEGATIVE);
1121 
1122  opts.optionalB("alch", "alchDecouple", "Enable alchemical decoupling?",
1123  &alchDecouple, FALSE);
1124  opts.optionalB("alch", "alchBondDecouple", "Enable decoupling of purely "
1125  "alchemical bonds?", &alchBondDecouple, FALSE);
1126 
1127  // parameters for alchemical analysis options
1128  opts.optional("alch", "alchType", "Which alchemical method to use?",
1129  PARSE_STRING);
1130  opts.optional("alch", "alchLambda2", "Alchemical coupling comparison value",
1131  &alchLambda2, -1);
1132  opts.optional("alch", "alchLambdaIDWS", "Alchemical coupling comparison value for interleaved double-wide sampling",
1133  &alchLambdaIDWS, -1);
1134  opts.optional("alch", "alchLambdaFreq",
1135  "Frequency of increasing coupling parameter value", &alchLambdaFreq, 0);
1136  opts.range("alchLambdaFreq", NOT_NEGATIVE);
1137  opts.optional("alch", "alchSwitchType", "Switching type flag",
1138  PARSE_STRING);
1139  opts.optional("alch", "alchEquilSteps", "Equilibration steps, before "
1140  "data collection in the alchemical window", &alchEquilSteps, 0);
1141  opts.range("alchEquilSteps", NOT_NEGATIVE);
1142 
1143  opts.optionalB("alch", "alchEnsembleAvg", "Ensemble Average in use?",
1144  &alchEnsembleAvg, TRUE);
1145 //fepe
1146 
1147  opts.optionalB("main", "les", "Is locally enhanced sampling enabled?",
1148  &lesOn, FALSE);
1149  opts.require("les", "lesFactor", "Local enhancement factor", &lesFactor);
1150  opts.optional("les", "lesFile", "PDB file with enhancement flags "
1151  "default is the input PDB file", PARSE_STRING);
1152  opts.optional("les", "lesCol", "Column in the lesFile with the "
1153  "enhancement flag", PARSE_STRING);
1154  opts.optionalB("les", "lesReduceTemp", "Reduce enhanced atom temperature?",
1155  &lesReduceTemp, FALSE);
1156  opts.optionalB("les", "lesReduceMass", "Reduce enhanced atom mass?",
1157  &lesReduceMass, FALSE);
1158 
1159  // REST2 (replica exchange solute tempering) parameters
1160  opts.optionalB("main", "soluteScaling",
1161  "Is replica exchange solute tempering enabled?",
1162  &soluteScalingOn, FALSE);
1163  opts.require("soluteScaling", "soluteScalingFactor",
1164  "Solute scaling factor",
1165  &soluteScalingFactor);
1166  opts.range("soluteScalingFactor", NOT_NEGATIVE);
1167  opts.optional("soluteScaling", "soluteScalingFactorCharge",
1168  "Solute scaling factor for electrostatic interactions",
1169  &soluteScalingFactorCharge);
1170  opts.range("soluteScalingFactorCharge", NOT_NEGATIVE);
1171  opts.optional("soluteScaling", "soluteScalingFactorVdw",
1172  "Solute scaling factor for van der Waals interactions",
1173  &soluteScalingFactorVdw);
1174  opts.range("soluteScalingFactorVdw", NOT_NEGATIVE);
1175  opts.optional("soluteScaling", "soluteScalingFile",
1176  "PDB file with scaling flags; if undefined, defaults to main PDB file",
1177  PARSE_STRING);
1178  opts.optional("soluteScaling", "soluteScalingCol",
1179  "Column in the soluteScalingFile providing the scaling flag",
1180  PARSE_STRING);
1181  opts.optionalB("main", "soluteScalingAll",
1182  "Apply scaling also to bond and angle interactions?",
1183  &soluteScalingAll, FALSE);
1184 
1185  // Drude oscillators
1186  opts.optionalB("main", "drude", "Perform integration of Drude oscillators?",
1187  &drudeOn, FALSE);
1188  opts.require("drude", "drudeTemp", "Temperature for freezing "
1189  "Drude oscillators", &drudeTemp);
1190  opts.range("drudeTemp", NOT_NEGATIVE);
1191  opts.units("drudeTemp", N_KELVIN);
1192  opts.optional("drude", "drudeDamping", "Damping coefficient (1/ps) for "
1193  "Drude oscillators", &drudeDamping);
1194  opts.range("drudeDamping", POSITIVE);
1195  opts.optional("drude", "drudeNbtholeCut", "Nonbonded Thole interactions "
1196  "interaction radius", &drudeNbtholeCut, 5.0);
1197  opts.range("drudeNbtholeCut", POSITIVE);
1198  opts.optionalB("drude", "drudeHardWall", "Apply maximum Drude bond length "
1199  "restriction?", &drudeHardWallOn, TRUE);
1200  opts.optional("drude", "drudeBondLen", "Drude oscillator bond length "
1201  "beyond which to apply restraint", &drudeBondLen, 0.25);
1202  opts.range("drudeBondLen", POSITIVE);
1203  opts.optional("drude", "drudeBondConst", "Drude oscillator restraining "
1204  "force constant", &drudeBondConst, 40000.0);
1205  opts.range("drudeBondConst", POSITIVE);
1206 
1207  // Pair interaction calculations
1208  opts.optionalB("main", "pairInteraction",
1209  "Are pair interactions calculated?", &pairInteractionOn, FALSE);
1210  opts.optional("pairInteraction", "pairInteractionFile",
1211  "PDB files with interaction flags " "default is the input PDB file",
1212  PARSE_STRING);
1213  opts.optional("pairInteraction", "pairInteractionCol",
1214  "Column in the pairInteractionFile with the interaction flags",
1215  PARSE_STRING);
1216  opts.require("pairInteraction", "pairInteractionGroup1",
1217  "Flag for interaction group 1", &pairInteractionGroup1);
1218  opts.optional("pairInteraction", "pairInteractionGroup2",
1219  "Flag for interaction group 2", &pairInteractionGroup2, -1);
1220  opts.optionalB("pairInteraction", "pairInteractionSelf",
1221  "Compute only within-group interactions?", &pairInteractionSelf,
1222  FALSE);
1223  // Options for CG simulations
1224  opts.optionalB("main", "cosAngles", "Are some angles cosine-based?", &cosAngles, FALSE);
1225 
1226 
1227  // Dihedral angle dynamics
1228  opts.optionalB("main", "globalTest", "Should global integration (for development) be used?",
1229  &globalOn, FALSE);
1230  opts.optionalB("main", "dihedral", "Should dihedral angle dynamics be performed?",
1231  &dihedralOn, FALSE);
1232  COLDOn = FALSE;
1233  opts.optionalB("dihedral", "COLD", "Should overdamped Langevin dynamics be performed?",
1234  &COLDOn, FALSE);
1235  opts.require("COLD", "COLDTemp", "Temperature for heat bath in COLD",
1236  &COLDTemp);
1237  opts.range("COLDTemp", NOT_NEGATIVE);
1238  opts.units("COLDTemp", N_KELVIN);
1239  opts.require("COLD", "COLDRate", "Damping rate for COLD",
1240  &COLDRate, 3000.0);
1241  opts.range("COLDRate", NOT_NEGATIVE);
1242 
1243  // Get the parameters for temperature coupling
1244  opts.optionalB("main", "tcouple",
1245  "Should temperature coupling be performed?",
1246  &tCoupleOn, FALSE);
1247  opts.require("tcouple", "tCoupleTemp",
1248  "Temperature for temperature coupling", &tCoupleTemp);
1249  opts.range("tCoupleTemp", NOT_NEGATIVE);
1250  opts.units("tCoupleTemp", N_KELVIN);
1251  opts.optional("tCouple", "tCoupleFile", "PDB file with temperature "
1252  "coupling terms (B(i)) (default is the PDB input file)",
1253  PARSE_STRING);
1254  opts.optional("tCouple", "tCoupleCol", "Column in the tCoupleFile "
1255  "containing the temperature coupling term B(i);\n"
1256  "default is 'O'", PARSE_STRING);
1257 
1258  opts.optionalB("main", "stochRescale",
1259  "Should stochastic velocity rescaling be performed?",
1260  &stochRescaleOn, FALSE);
1261  opts.require("stochRescale", "stochRescaleTemp",
1262  "Temperature for stochastic velocity rescaling",
1263  &stochRescaleTemp);
1264  opts.range("stochRescaleTemp", NOT_NEGATIVE);
1265  opts.units("stochRescaleTemp", N_KELVIN);
1266  opts.require("stochRescale", "stochRescalePeriod",
1267  "Time scale for stochastic velocity rescaling (ps)",
1268  &stochRescalePeriod);
1269  opts.range("stochRescalePeriod", POSITIVE);
1270  opts.optional("stochRescale", "stochRescaleFreq",
1271  "Number of steps between stochastic rescalings",
1272  &stochRescaleFreq);
1273  opts.range("stochRescaleFreq", POSITIVE);
1274  opts.optionalB("stochRescale", "stochRescaleHeat",
1275  "Should heat transfer and work be computed?", &stochRescaleHeat, FALSE);
1276 
1277  opts.optional("main", "rescaleFreq", "Number of steps between "
1278  "velocity rescaling", &rescaleFreq);
1279  opts.range("rescaleFreq", POSITIVE);
1280  opts.optional("main", "rescaleTemp", "Target temperature for velocity rescaling",
1281  &rescaleTemp);
1282  opts.range("rescaleTemp", NOT_NEGATIVE);
1283  opts.units("rescaleTemp", N_KELVIN);
1284 
1285  opts.optional("main", "reassignFreq", "Number of steps between "
1286  "velocity reassignment", &reassignFreq);
1287  opts.range("reassignFreq", POSITIVE);
1288  opts.optional("main", "reassignTemp", "Target temperature for velocity reassignment",
1289  &reassignTemp);
1290  opts.range("reassignTemp", NOT_NEGATIVE);
1291  opts.units("reassignTemp", N_KELVIN);
1292  opts.optional("main", "reassignIncr", "Temperature increment for velocity reassignment",
1293  &reassignIncr);
1294  opts.units("reassignIncr", N_KELVIN);
1295  opts.optional("main", "reassignHold", "Final holding temperature for velocity reassignment",
1296  &reassignHold);
1297  opts.range("reassignHold", NOT_NEGATIVE);
1298  opts.units("reassignHold", N_KELVIN);
1299 
1301  opts.optionalB("main", "useGroupPressure",
1302  "Use group rather than atomic quantities for pressure control?",
1303  &useGroupPressure, FALSE);
1304 
1306  opts.optionalB("main", "useFlexibleCell",
1307  "Use anisotropic cell fluctuation for pressure control?",
1308  &useFlexibleCell, FALSE);
1309 
1310  // Fix specific cell dimensions
1311  opts.optionalB("main", "fixCellDims",
1312  "Fix some cell dimensions?",
1313  &fixCellDims, FALSE);
1314 
1315  opts.optionalB("fixCellDims", "fixCellDimX",
1316  "Fix the X dimension?",
1317  &fixCellDimX, FALSE);
1318  opts.optionalB("fixCellDims", "fixCellDimY",
1319  "Fix the Y dimension?",
1320  &fixCellDimY, FALSE);
1321  opts.optionalB("fixCellDims", "fixCellDimZ",
1322  "Fix the Z dimension?",
1323  &fixCellDimZ, FALSE);
1324 
1326  opts.optionalB("main", "useConstantRatio",
1327  "Use constant X-Y ratio for pressure control?",
1328  &useConstantRatio, FALSE);
1329 
1331  opts.optionalB("main", "useConstantArea",
1332  "Use constant area for pressure control?",
1333  &useConstantArea, FALSE);
1334 
1336  opts.optionalB("main", "excludeFromPressure",
1337  "Should some atoms be excluded from pressure rescaling?",
1338  &excludeFromPressure, FALSE);
1339  opts.optional("excludeFromPressure", "excludeFromPressureFile",
1340  "PDB file for atoms to be excluded from pressure",
1341  PARSE_STRING);
1342  opts.optional("excludeFromPressure", "excludeFromPressureCol",
1343  "Column in the excludeFromPressureFile"
1344  "containing the flags (nonzero means excluded);\n"
1345  "default is 'O'", PARSE_STRING);
1346 
1348  opts.optionalB("main", "BerendsenPressure",
1349  "Should Berendsen pressure bath coupling be performed?",
1350  &berendsenPressureOn, FALSE);
1351  opts.require("BerendsenPressure", "BerendsenPressureTarget",
1352  "Target pressure for pressure coupling",
1353  &berendsenPressureTarget);
1354  // opts.units("BerendsenPressureTarget",);
1355  opts.require("BerendsenPressure", "BerendsenPressureCompressibility",
1356  "Isothermal compressibility for pressure coupling",
1357  &berendsenPressureCompressibility);
1358  // opts.units("BerendsenPressureCompressibility",);
1359  opts.require("BerendsenPressure", "BerendsenPressureRelaxationTime",
1360  "Relaxation time for pressure coupling",
1361  &berendsenPressureRelaxationTime);
1362  opts.range("BerendsenPressureRelaxationTime", POSITIVE);
1363  opts.units("BerendsenPressureRelaxationTime", N_FSEC);
1364  opts.optional("BerendsenPressure", "BerendsenPressureFreq",
1365  "Number of steps between volume rescaling",
1366  &berendsenPressureFreq, 1);
1367  opts.range("BerendsenPressureFreq", POSITIVE);
1368 
1370  opts.optionalB("main", "LangevinPiston",
1371  "Should Langevin piston pressure control be used?",
1372  &langevinPistonOn, FALSE);
1373  opts.optionalB("LangevinPiston", "LangevinPistonBarrier",
1374  "Should Langevin piston barrier be used?",
1375  &langevinPistonBarrier, TRUE);
1376  opts.require("LangevinPiston", "LangevinPistonTarget",
1377  "Target pressure for pressure control",
1378  &langevinPistonTarget);
1379  opts.require("LangevinPiston", "LangevinPistonPeriod",
1380  "Oscillation period for pressure control",
1381  &langevinPistonPeriod);
1382  opts.range("LangevinPistonPeriod", POSITIVE);
1383  opts.units("LangevinPistonPeriod", N_FSEC);
1384  opts.require("LangevinPiston", "LangevinPistonDecay",
1385  "Decay time for pressure control",
1386  &langevinPistonDecay);
1387  opts.range("LangevinPistonDecay", POSITIVE);
1388  opts.units("LangevinPistonDecay", N_FSEC);
1389  opts.require("LangevinPiston", "LangevinPistonTemp",
1390  "Temperature for pressure control piston",
1391  &langevinPistonTemp);
1392  opts.range("LangevinPistonTemp", POSITIVE);
1393  opts.units("LangevinPistonTemp", N_KELVIN);
1394  opts.optional("LangevinPiston", "StrainRate",
1395  "Initial strain rate for pressure control (x y z)",
1396  &strainRate);
1397 
1398  // Multigrator temperature and/or pressure control
1399  opts.optionalB("main", "Multigrator",
1400  "Should multigrator temperature and/or pressure control be used?",
1401  &multigratorOn, FALSE);
1402  opts.require("Multigrator", "MultigratorPressureTarget",
1403  "Target pressure for pressure coupling",
1404  &multigratorPressureTarget);
1405  opts.require("Multigrator", "MultigratorTemperatureTarget",
1406  "Target temperature for temperature coupling",
1407  &multigratorTemperatureTarget);
1408  opts.require("Multigrator", "MultigratorPressureFreq",
1409  "Number of steps between pressure control moves",
1410  &multigratorPressureFreq);
1411  opts.range("MultigratorPressureFreq", POSITIVE);
1412  opts.optional("Multigrator", "MultigratorPressureRelaxationTime",
1413  "Relaxation time for pressure coupling is fs",
1414  &multigratorPressureRelaxationTime, 30000);
1415  opts.range("MultigratorPressureRelaxationTime", POSITIVE);
1416  opts.units("MultigratorPressureRelaxationTime", N_FSEC);
1417  opts.optional("Multigrator", "MultigratorTemperatureRelaxationTime",
1418  "Relaxation time for temperature coupling is fs",
1419  &multigratorTemperatureRelaxationTime, 1000);
1420  opts.range("MultigratorTemperatureRelaxationTime", POSITIVE);
1421  opts.units("MultigratorTemperatureRelaxationTime", N_FSEC);
1422  opts.require("Multigrator", "MultigratorTemperatureFreq",
1423  "Number of steps between temperature control moves",
1424  &multigratorTemperatureFreq);
1425  opts.range("MultigratorTemperatureFreq", POSITIVE);
1426  opts.optional("Multigrator", "MultigratorNoseHooverChainLength",
1427  "Nose-Hoover chain length",
1428  &multigratorNoseHooverChainLength, 4);
1429  opts.range("MultigratorNoseHooverChainLength", POSITIVE);
1430 
1432  opts.optional("main", "SurfaceTensionTarget",
1433  "Surface tension in the x-y plane",
1434  &surfaceTensionTarget, 0);
1435 
1437  opts.optionalB("main", "pressureprofile", "Compute pressure profile?",
1438  &pressureProfileOn, FALSE);
1439  opts.require("pressureprofile", "pressureprofileslabs",
1440  "Number of pressure profile slabs", &pressureProfileSlabs, 10);
1441  opts.optional("pressureprofile", "pressureprofilefreq",
1442  "How often to store profile data", &pressureProfileFreq, 1);
1443  opts.optional("pressureprofile", "pressureProfileAtomTypes",
1444  "Number of pressure profile atom types", &pressureProfileAtomTypes, 1);
1445  opts.range("pressureProfileAtomTypes", POSITIVE);
1446  opts.optional("pressureProfile", "pressureProfileAtomTypesFile",
1447  "PDB files with pressure profile atom types" "default is the input PDB file",
1448  PARSE_STRING);
1449  opts.optional("pressureProfile", "pressureProfileAtomTypesCol",
1450  "Column in the pressureProfileAtomTypesFile with the atom types ",
1451  PARSE_STRING);
1452  opts.optionalB("pressureProfile", "pressureProfileEwald",
1453  "Compute Ewald contribution to pressure profile",
1454  &pressureProfileEwaldOn, FALSE);
1455  opts.optional("pressureProfile", "pressureProfileEwaldX",
1456  "Ewald grid size X", &pressureProfileEwaldX, 10);
1457  opts.range("pressureProfileEwaldX", POSITIVE);
1458  opts.optional("pressureProfile", "pressureProfileEwaldY",
1459  "Ewald grid size Y", &pressureProfileEwaldY, 10);
1460  opts.range("pressureProfileEwaldY", POSITIVE);
1461  opts.optional("pressureProfile", "pressureProfileEwaldZ",
1462  "Ewald grid size Z", &pressureProfileEwaldZ, 10);
1463  opts.range("pressureProfileEwaldZ", POSITIVE);
1464 
1466  opts.optionalB("main", "accelMD", "Perform acclerated MD?", &accelMDOn, FALSE);
1467  opts.optional("accelMD", "accelMDFirstStep", "First accelMD step", &accelMDFirstStep, 0);
1468  opts.range("accelMDFirstStep", NOT_NEGATIVE);
1469  opts.optional("accelMD", "accelMDLastStep", "Last accelMD step", &accelMDLastStep, 0);
1470  opts.range("accelMDLastStep", NOT_NEGATIVE);
1471  opts.optional("accelMD", "accelMDOutFreq", "Frequency of accelMD output", &accelMDOutFreq, 1);
1472  opts.range("accelMDOutFreq", POSITIVE);
1473  opts.optionalB("accelMD", "accelMDdihe", "Apply boost to dihedral potential", &accelMDdihe, TRUE);
1474  opts.optionalB("accelMD", "accelMDDebugOn", "Debugging accelMD", &accelMDDebugOn, FALSE);
1475  opts.optional("accelMD", "accelMDE","E for AMD", &accelMDE);
1476  opts.units("accelMDE", N_KCAL);
1477  opts.optional("accelMD", "accelMDalpha","alpha for AMD", &accelMDalpha);
1478  opts.units("accelMDalpha", N_KCAL);
1479  opts.range("accelMDalpha", POSITIVE);
1480  opts.optionalB("accelMD", "accelMDdual", "Apply dual boost", &accelMDdual, FALSE);
1481  opts.optional("accelMDdual", "accelMDTE","E for total potential under accelMDdual mode", &accelMDTE);
1482  opts.units("accelMDTE", N_KCAL);
1483  opts.optional("accelMDdual", "accelMDTalpha","alpha for total potential under accelMDdual mode", &accelMDTalpha);
1484  opts.units("accelMDTalpha", N_KCAL);
1485  opts.range("accelMDTalpha", POSITIVE);
1486  // GaMD parameters
1487  opts.optionalB("accelMD", "accelMDG", "Perform Gaussian accelMD calculation?", &accelMDG, FALSE);
1488  opts.optional("accelMDG", "accelMDGiE", "Flag to set the mode iE in Gaussian accelMD", &accelMDGiE, 1);
1489  opts.optional("accelMDG", "accelMDGcMDSteps", "Number of cMD steps", &accelMDGcMDSteps, 1000000);
1490  opts.range("accelMDGcMDSteps", NOT_NEGATIVE);
1491  opts.optional("accelMDG", "accelMDGEquiSteps", "Number of equilibration steps after adding boost potential", &accelMDGEquiSteps, 1000000);
1492  opts.range("accelMDGEquiSteps", NOT_NEGATIVE);
1493  opts.require("accelMDG", "accelMDGcMDPrepSteps", "Number of preparation cMD steps", &accelMDGcMDPrepSteps, 200000);
1494  opts.range("accelMDGcMDPrepSteps", NOT_NEGATIVE);
1495  opts.require("accelMDG", "accelMDGEquiPrepSteps", "Number of preparation equilibration steps", &accelMDGEquiPrepSteps, 200000);
1496  opts.range("accelMDGEquiPrepSteps", NOT_NEGATIVE);
1497  opts.optional("accelMDG", "accelMDGStatWindow", "Number of steps to calculate avg and std", &accelMDGStatWindow, -1);
1498  opts.optional("accelMDG", "accelMDGSigma0P", "Upper limit of std of total potential", &accelMDGSigma0P, 6.0);
1499  opts.units("accelMDGSigma0P", N_KCAL);
1500  opts.range("accelMDGSigma0P", NOT_NEGATIVE);
1501  opts.optional("accelMDG", "accelMDGSigma0D", "Upper limit of std of dihedral potential", &accelMDGSigma0D, 6.0);
1502  opts.units("accelMDGSigma0D", N_KCAL);
1503  opts.range("accelMDGSigma0D", NOT_NEGATIVE);
1504  opts.optionalB("accelMDG", "accelMDGRestart", "Flag to set use restart file in Gaussian accelMD", &accelMDGRestart, FALSE);
1505  opts.require("accelMDGRestart", "accelMDGRestartFile", "Restart file name for Gaussian accelMD", accelMDGRestartFile);
1506  opts.optionalB("accelMDG", "accelMDGresetVaftercmd", "Flag to reset potential after accelMDGcMDSteps steps",
1507  &accelMDGresetVaftercmd, FALSE);
1508 
1509  // Adaptive Temperature Sampling (adaptTemp) parameters
1510  opts.optionalB("main", "adaptTempMD", "Perform adaptive temperature sampling", &adaptTempOn, FALSE);
1511  opts.optional("adaptTempMD", "adaptTempFirstStep", "First adaptTemp step", &adaptTempFirstStep, 0);
1512  opts.range("adaptTempFirstStep", NOT_NEGATIVE);
1513  opts.optional("adaptTempMD", "adaptTempLastStep", "Last adaptTemp step", &adaptTempLastStep, 0);
1514  opts.range("adaptTempLastStep", NOT_NEGATIVE);
1515  opts.optional("adaptTempMD", "adaptTempOutFreq", "Frequency of adaptTemp output", &adaptTempOutFreq, 10);
1516  opts.range("adaptTempOutFreq", POSITIVE);
1517  opts.optional("adaptTempMD", "adaptTempFreq", "Frequency of writing average energies to adaptTempOutFile", &adaptTempFreq, 10);
1518  opts.range("adaptTempFreq", POSITIVE);
1519  opts.optionalB("adaptTempMD", "adaptTempDebug", "Print debug output for adaptTemp", &adaptTempDebug, FALSE);
1520  opts.optional("adaptTempMD", "adaptTempTmin","Minimun temperature for adaptTemp", &adaptTempTmin);
1521  opts.units("adaptTempTmin", N_KELVIN);
1522  opts.range("adaptTempTmin", POSITIVE);
1523  opts.optional("adaptTempMD", "adaptTempTmax","Maximum temperature for adaptTemp", &adaptTempTmax);
1524  opts.units("adaptTempTmax", N_KELVIN);
1525  opts.range("adaptTempTmax", POSITIVE);
1526  opts.optional("adaptTempMD", "adaptTempBins","Number of bins to store average energies", &adaptTempBins,0);
1527  opts.range("adaptTempBins", NOT_NEGATIVE);
1528  opts.optional("adaptTempMD", "adaptTempDt", "Integration timestep for Temp. updates", &adaptTempDt, 0.0001);
1529  opts.units("adaptTempDt", N_FSEC);
1530  opts.range("adaptTempDt", NOT_NEGATIVE);
1531  opts.optional("adaptTempMD", "adaptTempAutoDt", "Average temperature update in percent of temperature range", &adaptTempAutoDt, 0.0);
1532  opts.range("adaptTempAutoDt", NOT_NEGATIVE);
1533  opts.optional("adaptTempMD", "adaptTempCgamma", "Adaptive bin averaging constant", &adaptTempCgamma, 0.1);
1534  opts.range("adaptTempCgamma", NOT_NEGATIVE);
1535  opts.optionalB("adaptTempMD","adaptTempLangevin","Send adaptTemp temperature to langevin thermostat",&adaptTempLangevin,TRUE);
1536  opts.optionalB("adaptTempMD","adaptTempRescaling","Send adaptTemp temperature to velocity rescaling thermostat", &adaptTempRescale,TRUE);
1537  opts.optional("adaptTempMD", "adaptTempInFile", "File containing restart information for adaptTemp", adaptTempInFile);
1538  opts.optional("adaptTempMD", "adaptTempRestartFile", "File for writing adaptTemp restart information", adaptTempRestartFile);
1539  opts.require("adaptTempRestartFile","adaptTempRestartFreq", "Frequency of writing restart file", &adaptTempRestartFreq,0);
1540  opts.range("adaptTempRestartFreq",NOT_NEGATIVE);
1541  opts.optionalB("adaptTempMD", "adaptTempRandom", "Randomly assign a temperature if we step out of range", &adaptTempRandom, FALSE);
1542 }
1543 
1544 void SimParameters::config_parser_constraints(ParseOptions &opts) {
1545 
1547  opts.optionalB("main", "fixedatoms", "Are there fixed atoms?",
1548  &fixedAtomsOn, FALSE);
1549  opts.optionalB("fixedatoms", "fixedAtomsForces",
1550  "Calculate forces between fixed atoms? (Required to unfix during run.)",
1551  &fixedAtomsForces, FALSE);
1552  opts.optional("fixedatoms", "fixedAtomsFile", "PDB file with flags for "
1553  "fixed atoms (default is the PDB input file)",
1554  PARSE_STRING);
1555  opts.optional("fixedatoms", "fixedAtomsCol", "Column in the fixedAtomsFile "
1556  "containing the flags (nonzero means fixed);\n"
1557  "default is 'O'", PARSE_STRING);
1558  opts.optional("fixedatoms", "fixedAtomListFile", "the text input file for fixed atoms "
1559  "used for parallel input IO", PARSE_STRING);
1560  opts.optionalB("fixedatoms", "fixedAtomsForceOutput",
1561  "Do we write out forces acting on fixed atoms?",
1562  &fixedAtomsForceOutput, FALSE);
1563 
1565  opts.optionalB("main", "constraints", "Are harmonic constraints active?",
1566  &constraintsOn, FALSE);
1567  opts.require("constraints", "consexp", "Exponent for harmonic potential",
1568  &constraintExp, 2);
1569  opts.range("consexp", POSITIVE);
1570 #ifndef MEM_OPT_VERSION
1571  opts.require("constraints", "consref", "PDB file containing reference "
1572  "positions",
1573  PARSE_STRING);
1574  opts.require("constraints", "conskfile", "PDB file containing force "
1575  "constaints in one of the columns", PARSE_STRING);
1576  opts.require("constraints", "conskcol", "Column of conskfile to use "
1577  "for the force constants", PARSE_STRING);
1578 #else
1579  opts.require("constraints", "consAtomListFile", "the text input file for constrained atoms "
1580  "used for parallel input IO", PARSE_STRING);
1581 #endif
1582  opts.require("constraints", "constraintScaling", "constraint scaling factor",
1583  &constraintScaling, 1.0);
1584  opts.range("constraintScaling", NOT_NEGATIVE);
1585 
1586 
1587 
1588  //****** BEGIN selective restraints (X,Y,Z) changes
1589 
1591  opts.optionalB("constraints", "selectConstraints",
1592  "Restrain only selected Cartesian components of the coordinates?",
1593  &selectConstraintsOn, FALSE);
1594  opts.optionalB("selectConstraints", "selectConstrX",
1595  "Restrain X components of coordinates ", &constrXOn, FALSE);
1596  opts.optionalB("selectConstraints", "selectConstrY",
1597  "Restrain Y components of coordinates ", &constrYOn, FALSE);
1598  opts.optionalB("selectConstraints", "selectConstrZ",
1599  "Restrain Z components of coordinates ", &constrZOn, FALSE);
1600  //****** END selective restraints (X,Y,Z) changes
1601 
1602  // spherical constraints
1603  opts.optionalB("constraints", "sphericalConstraints",
1604  "Restrain only radial spherical component of the coordinates?",
1605  &sphericalConstraintsOn, FALSE);
1606  opts.optional("sphericalConstraints", "sphericalConstrCenter",
1607  "Center of spherical constraints", &sphericalConstrCenter);
1608 
1609  //****** BEGIN moving constraints changes
1610 
1612  opts.optionalB("constraints", "movingConstraints",
1613  "Are some of the constraints moving?",
1614  &movingConstraintsOn, FALSE);
1615  opts.require("movingConstraints", "movingConsVel",
1616  "Velocity of the movement, A/timestep", &movingConsVel);
1617  //****** END moving constraints changes
1618 
1619  // BEGIN rotating constraints changes
1620  opts.optionalB("constraints", "rotConstraints",
1621  "Are the constraints rotating?",
1622  &rotConstraintsOn, FALSE);
1623  opts.require("rotConstraints", "rotConsAxis",
1624  "Axis of rotation", &rotConsAxis);
1625  opts.require("rotConstraints", "rotConsPivot",
1626  "Pivot point of rotation",
1627  &rotConsPivot);
1628  opts.require("rotConstraints", "rotConsVel",
1629  "Velocity of rotation, deg/timestep", &rotConsVel);
1630 
1631  // END rotating constraints changes
1632 
1633  // external command forces
1634  opts.optionalB("main", "extForces", "External command forces?",
1635  &extForcesOn, FALSE);
1636  opts.require("extForces", "extForcesCommand",
1637  "External forces command", extForcesCommand);
1638  opts.require("extForces", "extCoordFilename",
1639  "External forces coordinate filename", extCoordFilename);
1640  opts.require("extForces", "extForceFilename",
1641  "External forces force filename", extForceFilename);
1642 
1643 
1644  // QM/MM forces
1645  opts.optionalB("main", "QMForces", "Apply QM forces?",
1646  &qmForcesOn, FALSE);
1647  opts.require("QMForces", "QMSoftware",
1648  "software whose format will be used for input/output", qmSoftware);
1649  opts.require("QMForces", "QMExecPath",
1650  "path to executable", qmExecPath);
1651  opts.optional("QMForces", "QMChargeMode",
1652  "type of QM atom charges gathered from the QM software", qmChrgModeS);
1653  opts.require("QMForces", "QMColumn",
1654  "column defining QM and MM regions", qmColumn);
1655  opts.require("QMForces", "QMBaseDir",
1656  "base path and name for QM input and output (preferably in memory)", qmBaseDir);
1657  opts.optional("QMForces", "QMConfigLine",
1658  "Configuration line for QM (multiple inputs allowed)", PARSE_MULTIPLES);
1659  opts.optional("QMForces", "QMParamPDB",
1660  "PDB with QM parameters", qmParamPDB);
1661  opts.optional("QMForces", "QMPrepProc",
1662  "initial preparation executable", qmPrepProc);
1663  opts.optional("QMForces", "QMSecProc",
1664  "secondary executable", qmSecProc);
1665  opts.optional("QMForces", "QMCharge",
1666  "charge of the QM group", PARSE_MULTIPLES);
1667  opts.optionalB("QMForces", "QMChargeFromPSF",
1668  "gets charge of the QM group form PSF values", &qmChrgFromPSF, FALSE);
1669  opts.optional("QMForces", "QMMult",
1670  "multiplicity of the QM group", PARSE_MULTIPLES);
1671  opts.optional("QMForces", "QMLinkElement",
1672  "element of link atom", PARSE_MULTIPLES);
1673  opts.optionalB("QMForces", "QMReplaceAll",
1674  "replace all NAMD forces with QM forces", &qmReplaceAll, FALSE);
1675  opts.optional("QMForces", "QMPCStride",
1676  "frequency of selection of point charges", &qmPCSelFreq, 1);
1677  opts.range("QMPCStride", POSITIVE);
1678  opts.optionalB("QMForces", "QMNoPntChrg",
1679  "no point charges will be passed to the QM system(s)", &qmNoPC, FALSE);
1680  opts.optionalB("QMForces", "QMElecEmbed",
1681  "activates electrostatic embedding", &qmElecEmbed, TRUE);
1682  opts.optionalB("QMForces", "QMVdWParams",
1683  "use special VdW parameters for QM atoms", &qmVDW, FALSE);
1684  opts.optional("QMForces", "QMBondColumn",
1685  "column defining QM-MM bomnds", qmBondColumn);
1686  opts.optionalB("QMForces", "QMBondDist",
1687  "values in QMBondColumn defines the distance of new link atom", &qmBondDist, FALSE);
1688  opts.optional("QMForces", "QMBondValueType",
1689  "type of value in bond column: len or ratio", qmBondValueTypeS);
1690  opts.optional("QMForces", "QMBondScheme",
1691  "type of treatment given to QM-MM bonds.", qmBondSchemeS);
1692  opts.optional("QMForces", "QMenergyStride",
1693  "frequency of QM specific energy output (every x steps)", &qmEnergyOutFreq, 1);
1694  opts.optional("QMForces", "QMOutStride",
1695  "frequency of QM specific charge output (every x steps)", &qmOutFreq, 0);
1696  opts.range("QMOutStride", NOT_NEGATIVE);
1697  opts.optional("QMForces", "QMPositionOutStride",
1698  "frequency of QM specific position output (every x steps)", &qmPosOutFreq, 0);
1699  opts.range("QMPositionOutStride", NOT_NEGATIVE);
1700  opts.optional("QMForces", "QMSimsPerNode",
1701  "QM executions per node", &qmSimsPerNode, 1);
1702  opts.range("QMSimsPerNode", POSITIVE);
1703  opts.optionalB("QMForces", "QMSwitching",
1704  "apply switching to point charges.", &qmPCSwitchOn, FALSE);
1705  opts.optional("QMForces", "QMSwitchingType",
1706  "How are charges scaled down to be presented to QM groups.", qmPCSwitchTypeS);
1707  opts.optional("QMForces", "QMPointChargeScheme",
1708  "type of treatment given to the total sum of point charges.", qmPCSchemeS);
1709  opts.optionalB("QMForces", "QMCustomPCSelection",
1710  "custom and fixed selection of point charges per QM group.", &qmCustomPCSel, FALSE);
1711  opts.optional("QMForces", "QMCustomPCFile",
1712  "file with a selection of point charges for a single QM group", PARSE_MULTIPLES);
1713  opts.optionalB("QMForces", "QMLiveSolventSel",
1714  "Continuously update the selection of solvent molecules in QM groups", &qmLSSOn, FALSE);
1715  opts.optional("QMForces", "QMLSSFreq",
1716  "frequency of QM water selection update", &qmLSSFreq, 100);
1717  opts.range("QMLSSFreq", POSITIVE);
1718  opts.optional("QMForces", "QMLSSResname",
1719  "residue name for the solvent molecules (TIP3).", qmLSSResname);
1720  opts.optional("QMForces", "QMLSSMode",
1721  "mode of selection of point solvent molecules", qmLSSModeS);
1722  opts.optional("QMForces", "QMLSSRef",
1723  "for COM mode, defines reference for COM distance calculation", PARSE_MULTIPLES);
1724  opts.optionalB("QMForces", "QMCSMD",
1725  "Do we use Conditional SMD option?", &qmCSMD, FALSE);
1726  opts.optional("QMForces", "QMCSMDFile",
1727  "File for Conditional SMD information",qmCSMDFile);
1728 
1729  //print which bad contacts are being moved downhill
1730  opts.optionalB("main", "printBadContacts", "Print atoms with huge forces?",
1731  &printBadContacts, FALSE);
1732 
1733  /* GBIS generalized born implicit solvent*/
1734 
1735  opts.optionalB("main", "GBIS", "Use GB implicit solvent?",
1736  &GBISOn, FALSE);
1737  opts.optionalB("main", "GBISSer", "Use GB implicit solvent?",
1738  &GBISserOn, FALSE);
1739 
1740  opts.optional("GBIS", "solventDielectric",
1741  "Solvent Dielectric", &solvent_dielectric, 78.5);
1742  opts.optional("GBIS", "intrinsicRadiusOffset",
1743  "Coulomb Radius Offset", &coulomb_radius_offset, 0.09);
1744  opts.optional("GBIS", "ionConcentration",
1745  "Ion Concentration", &ion_concentration, 0.2); //0.2 mol/L
1746  opts.optional("GBIS", "GBISDelta",
1747  "delta from GBOBC", &gbis_delta, 1.0); //0.8 or 1.0
1748  opts.optional("GBIS", "GBISBeta",
1749  "beta from GBOBC", &gbis_beta, 0.8); //0.0 or 0.8
1750  opts.optional("GBIS", "GBISGamma",
1751  "gamma from GBOBC", &gbis_gamma, 4.85);//2.290912 or 4.85
1752  opts.optional("GBIS", "alphaCutoff",
1753  "cutoff for calculating effective born radius", &alpha_cutoff, 15);
1754  opts.optional("GBIS", "alphaMax",
1755  "maximum allowable born radius", &alpha_max, 30);
1756  opts.optional("GBIS", "fsMax",
1757  "maximum screened intrinsic radius", &fsMax, 1.728);
1758 
1759  opts.optionalB("main", "SASA", "Use Linear Combination of Pairwise Overlaps (LCPO) for calculating SASA",
1760  &LCPOOn, FALSE);
1761  opts.optional("SASA", "surfaceTension",
1762  "Surfce Tension for SASA (kcal/mol/Ang^2)", &surface_tension, 0.005);
1763 
1764  //****** BEGIN SMD constraints changes
1765 
1766  // SMD constraints
1767  opts.optionalB("main", "SMD",
1768  "Do we use SMD option?",
1769  &SMDOn, FALSE);
1770  opts.require("SMD", "SMDVel",
1771  "Velocity of the movement, A/timestep", &SMDVel);
1772  opts.range("SMDVel", NOT_NEGATIVE);
1773  opts.require("SMD", "SMDDir",
1774  "Direction of movement", &SMDDir);
1775  opts.require("SMD", "SMDk",
1776  "Elastic constant for SMD", &SMDk);
1777  opts.optional("SMD", "SMDk2",
1778  "Transverse elastic constant for SMD", &SMDk2, 0);
1779  opts.range("SMDk", NOT_NEGATIVE);
1780  opts.range("SMDk2", NOT_NEGATIVE);
1781  opts.require("SMD", "SMDFile",
1782  "File for SMD information",
1783  SMDFile);
1784  opts.optional("SMD", "SMDOutputFreq",
1785  "Frequency of output",
1786  &SMDOutputFreq, 1);
1787  opts.range("SMDOutputFreq", POSITIVE);
1788 
1789  //****** END SMD constraints changes
1790 
1791  //****** BEGIN tabulated energies section
1792  opts.optionalB("main", "tabulatedEnergies", "Do we get energies from a table?", &tabulatedEnergies, FALSE);
1793 // opts.require("tabulatedEnergies", "tableNumTypes","Number of types for energy tabulation", &tableNumTypes);
1794  opts.require("tabulatedEnergies", "tabulatedEnergiesFile", "File containing energy table", tabulatedEnergiesFile);
1795  opts.require("tabulatedEnergies", "tableInterpType", "Cubic or linear interpolation", tableInterpType);
1796 
1797  // TMD parameters
1798  opts.optionalB("main", "TMD", "Perform Targeted MD?", &TMDOn, FALSE);
1799  opts.optional("TMD", "TMDk", "Elastic constant for TMD", &TMDk, 0);
1800  opts.range("TMDk", NOT_NEGATIVE);
1801  opts.require("TMD", "TMDFile", "File for TMD information", TMDFile);
1802  opts.optionalB("TMD", "TMDDiffRMSD", "Restrain Difference between the RMSD from two structures", &TMDDiffRMSD, FALSE);
1803  opts.require("TMDDiffRMSD", "TMDFile2", "Second file for TMD information", TMDFile2);
1804 
1805  opts.optional("TMD", "TMDOutputFreq", "Frequency of TMD output",
1806  &TMDOutputFreq, 1);
1807  opts.range("TMDOutputFreq", POSITIVE);
1808  opts.require("TMD", "TMDLastStep", "Last TMD timestep", &TMDLastStep);
1809  opts.range("TMDLastStep", POSITIVE);
1810  opts.optional("TMD", "TMDFirstStep", "First TMD step (default 0)", &TMDFirstStep, 0);
1811  opts.optional("TMD", "TMDInitialRMSD", "Target RMSD at first TMD step (default -1 to use initial coordinates)", &TMDInitialRMSD);
1812  TMDInitialRMSD = -1;
1813  opts.optional("TMD", "TMDFinalRMSD", "Target RMSD at last TMD step (default 0 )", &TMDFinalRMSD, 0);
1814  opts.range("TMDInitialRMSD", NOT_NEGATIVE);
1815  // End of TMD parameters
1816 
1817  // Symmetry restraint parameters
1818  opts.optionalB("main", "symmetryRestraints", "Enable symmetry restraints?", &symmetryOn, FALSE);
1819  opts.optional("symmetryRestraints", "symmetryk", "Elastic constant for symmetry restraints", &symmetryk, 0);
1820  opts.range("symmetryk", NOT_NEGATIVE);
1821  opts.optional("symmetryRestraints", "symmetrykfile", "PDB file specifying force contants on a per-atom basis", PARSE_MULTIPLES);
1822  opts.optionalB("symmetryRestraints", "symmetryScaleForces", "Scale applied forces over time?", &symmetryScaleForces, FALSE);
1823  opts.require("symmetryRestraints", "symmetryFile", "File for symmetry information", PARSE_MULTIPLES);
1824  opts.optional("symmetryRestraints", "symmetryMatrixFile", "File(s) for transfromation matrices", PARSE_MULTIPLES);
1825  opts.optional("symmetryRestraints", "symmetryLastStep", "Last symmetry timestep", &symmetryLastStep, -1);
1826  opts.optional("symmetryRestraints", "symmetryFirstStep", "First symmetry step (default 0)", &symmetryFirstStep, 0);
1827  opts.optional("symmetryRestraints", "symmetryLastFullStep", "Last full force symmetry timestep (default symmetryLastStep)", &symmetryLastFullStep, symmetryLastStep);
1828  opts.optional("symmetryRestraints", "symmetryFirstFullStep", "First full force symmetry step (default symmetryFirstStep)", &symmetryFirstFullStep, symmetryFirstStep);
1829  //End of symmetry restraint parameters.
1830 
1832  opts.optionalB("main", "tclForces", "Are Tcl global forces active?",
1833  &tclForcesOn, FALSE);
1834  opts.require("tclForces", "tclForcesScript",
1835  "Tcl script for global forces", PARSE_MULTIPLES);
1836 
1838  opts.optionalB("main", "tclBC", "Are Tcl boundary forces active?",
1839  &tclBCOn, FALSE);
1840  opts.require("tclBC", "tclBCScript",
1841  "Tcl script defining calcforces for boundary forces", PARSE_STRING);
1842  tclBCScript = 0;
1843  opts.optional("tclBC", "tclBCArgs", "Extra args for calcforces command",
1844  tclBCArgs);
1845  tclBCArgs[0] = 0;
1846 
1848  opts.optionalB("main", "miscForces", "Are misc global forces active?",
1849  &miscForcesOn, FALSE);
1850  opts.optional("miscForces", "miscForcesScript",
1851  "script for misc forces", PARSE_MULTIPLES);
1852 
1854  opts.optionalB("main", "freeEnergy", "Perform free energy perturbation?",
1855  &freeEnergyOn, FALSE);
1856  opts.require("freeEnergy", "freeEnergyConfig",
1857  "Configuration file for free energy perturbation", PARSE_MULTIPLES);
1858 
1860  opts.optionalB("main", "constantforce", "Apply constant force?",
1861  &consForceOn, FALSE);
1862  opts.optional("constantforce", "consForceFile",
1863  "Configuration file for constant forces", PARSE_STRING);
1864  opts.require("constantforce", "consForceScaling",
1865  "Scaling factor for constant forces", &consForceScaling, 1.0);
1866 
1868  opts.optionalB("main", "colvars", "Is the colvars module enabled?",
1869  &colvarsOn, FALSE);
1870  opts.optional("colvars", "colvarsConfig",
1871  "configuration for the collective variables", PARSE_STRING);
1872  opts.optional("colvars", "colvarsInput",
1873  "input restart file for the collective variables", PARSE_STRING);
1874 
1875 }
1876 
1877 #ifdef OPENATOM_VERSION
1878 void SimParameters::config_parser_openatom(ParseOptions &opts) {
1879  opts.optionalB("main", "openatom", "OpenAtom active?", &openatomOn, FALSE);
1880  opts.require("openatom", "openatomDriverFile", "What config file specifies openatom input parameters", PARSE_STRING);
1881  opts.require("openatom", "openatomPhysicsFile", "What structure file specifies openatom input system", PARSE_STRING);
1882  opts.require("openatom", "openatomPdbFile", "NAMD input file defining QM and MM regions", PARSE_STRING);
1883  opts.optional("openatom", "openatomCol", "Column in the openatomPdb with the QM/MM flag", PARSE_STRING);
1884 }
1885 #endif // OPENATOM_VERSION
1886 
1887 /* BEGIN gf */
1888 void SimParameters::config_parser_mgridforce(ParseOptions &opts) {
1890  opts.optionalB("main", "mgridforce", "Is Multiple gridforce active?",
1891  &mgridforceOn, FALSE);
1892  opts.optional("mgridforce", "mgridforcevolts", "Is Gridforce using Volts/eV as units?",
1893  PARSE_MULTIPLES);
1894  opts.require("mgridforce", "mgridforcescale", "Scale factor by which to multiply "
1895  "grid forces", PARSE_MULTIPLES);
1896  opts.require("mgridforce", "mgridforcefile", "PDB file containing force "
1897  "multipliers in one of the columns", PARSE_MULTIPLES);
1898  opts.require("mgridforce", "mgridforcecol", "Column of gridforcefile to "
1899  "use for force multiplier", PARSE_MULTIPLES);
1900  opts.optional("mgridforce", "mgridforcechargecol", "Column of gridforcefile to "
1901  "use for charge", PARSE_MULTIPLES);
1902  opts.require("mgridforce", "mgridforcepotfile", "Gridforce potential file",
1903  PARSE_MULTIPLES);
1904  opts.optional("mgridforce", "mgridforcecont1", "Use continuous grid "
1905  "in K1 direction?", PARSE_MULTIPLES);
1906  opts.optional("mgridforce", "mgridforcecont2", "Use continuous grid "
1907  "in K2 direction?", PARSE_MULTIPLES);
1908  opts.optional("mgridforce", "mgridforcecont3", "Use continuous grid "
1909  "in K3 direction?", PARSE_MULTIPLES);
1910  opts.optional("mgridforce", "mgridforcevoff", "Gridforce potential offsets",
1911  PARSE_MULTIPLES);
1912  opts.optional("mgridforce", "mgridforcelite", "Use Gridforce Lite?",
1913  PARSE_MULTIPLES);
1914  opts.optional("mgridforce", "mgridforcechecksize", "Check if grid exceeds PBC cell dimensions?", PARSE_MULTIPLES);
1915 }
1916 
1917 void SimParameters::config_parser_gridforce(ParseOptions &opts) {
1919  opts.optionalB("main", "gridforce", "Is Gridforce active?",
1920  &gridforceOn, FALSE);
1921  opts.optionalB("gridforce", "gridforcevolts", "Is Gridforce using Volts/eV as units?",
1922  &gridforceVolts, FALSE);
1923  opts.require("gridforce", "gridforcescale", "Scale factor by which to multiply "
1924  "grid forces", &gridforceScale);
1925  opts.require("gridforce", "gridforcefile", "PDB file containing force "
1926  "multipliers in one of the columns", PARSE_STRING);
1927  opts.require("gridforce", "gridforcecol", "Column of gridforcefile to "
1928  "use for force multiplier", PARSE_STRING);
1929  opts.optional("gridforce", "gridforcechargecol", "Column of gridforcefile to "
1930  "use for charge", PARSE_STRING);
1931  opts.require("gridforce", "gridforcepotfile", "Gridforce potential file",
1932  PARSE_STRING);
1933  opts.optionalB("gridforce", "gridforcecont1", "Use continuous grid "
1934  "in A1 direction?", &gridforceContA1, FALSE);
1935  opts.optionalB("gridforce", "gridforcecont2", "Use continuous grid "
1936  "in A2 direction?", &gridforceContA2, FALSE);
1937  opts.optionalB("gridforce", "gridforcecont3", "Use continuous grid "
1938  "in A3 direction?", &gridforceContA3, FALSE);
1939  opts.optional("gridforce", "gridforcevoff", "Gridforce potential offsets",
1940  &gridforceVOffset);
1941  opts.optionalB("gridforce", "gridforcelite", "Use Gridforce Lite?",
1942  &gridforceLite, FALSE);
1943  opts.optionalB("gridforce", "gridforcechecksize", "Check if grid exceeds PBC cell dimensions?",
1944  &gridforcechecksize, TRUE);
1945 }
1946 /* END gf */
1947 
1948 void SimParameters::config_parser_movdrag(ParseOptions &opts) {
1950  opts.optionalB("main", "movDragOn", "Do we apply moving drag?",
1951  &movDragOn, FALSE);
1952  opts.require("movDragOn", "movDragFile",
1953  "Main moving drag PDB file", movDragFile);
1954  opts.require("movDragOn", "movDragCol",
1955  "Main moving drag PDB column", PARSE_STRING);
1956  opts.require("movDragOn", "movDragGlobVel",
1957  "Global moving drag velocity (A/step)", &movDragGlobVel);
1958  opts.require("movDragOn", "movDragVelFile",
1959  "Moving drag linear velocity file", movDragVelFile);
1960 }
1961 
1962 void SimParameters::config_parser_rotdrag(ParseOptions &opts) {
1964  opts.optionalB("main", "rotDragOn", "Do we apply rotating drag?",
1965  &rotDragOn, FALSE);
1966  opts.require("rotDragOn", "rotDragFile",
1967  "Main rotating drag PDB file", rotDragFile);
1968  opts.require("rotDragOn", "rotDragCol",
1969  "Main rotating drag PDB column", PARSE_STRING);
1970  opts.require("rotDragOn", "rotDragAxisFile",
1971  "Rotating drag axis file", rotDragAxisFile);
1972  opts.require("rotDragOn", "rotDragPivotFile",
1973  "Rotating drag pivot point file", rotDragPivotFile);
1974  opts.require("rotDragOn", "rotDragGlobVel",
1975  "Global rotating drag angular velocity (deg/step)", &rotDragGlobVel);
1976  opts.require("rotDragOn", "rotDragVelFile",
1977  "Rotating drag angular velocity file", rotDragVelFile);
1978  opts.require("rotDragOn", "rotDragVelCol",
1979  "Rotating drag angular velocity column", PARSE_STRING);
1980 }
1981 
1982 void SimParameters::config_parser_constorque(ParseOptions &opts) {
1984  opts.optionalB("main", "consTorqueOn", "Do we apply \"constant\" torque?",
1985  &consTorqueOn, FALSE);
1986  opts.require("consTorqueOn", "consTorqueFile",
1987  "Main \"constant\" torque PDB file", consTorqueFile);
1988  opts.require("consTorqueOn", "consTorqueCol",
1989  "Main \"constant\" torque PDB column", PARSE_STRING);
1990  opts.require("consTorqueOn", "consTorqueAxisFile",
1991  "\"Constant\" torque axis file", consTorqueAxisFile);
1992  opts.require("consTorqueOn", "consTorquePivotFile",
1993  "\"Constant\" torque pivot point file", consTorquePivotFile);
1994  opts.require("consTorqueOn", "consTorqueGlobVal",
1995  "Global \"constant\" torque value (Kcal/(mol*A^2))", &consTorqueGlobVal);
1996  opts.require("consTorqueOn", "consTorqueValFile",
1997  "\"constant\" torque factors file", consTorqueValFile);
1998  opts.require("consTorqueOn", "consTorqueValCol",
1999  "\"constant\" torque factors column", PARSE_STRING);
2000 }
2001 
2002 void SimParameters::config_parser_boundary(ParseOptions &opts) {
2003 
2005  opts.optionalB("main", "sphericalBC", "Are spherical boundary counditions "
2006  "active?", &sphericalBCOn, FALSE);
2007  opts.require("sphericalBC", "sphericalBCCenter",
2008  "Center of spherical boundaries", &sphericalCenter);
2009  opts.require("sphericalBC", "sphericalBCr1", "Radius for first sphere "
2010  "potential", &sphericalBCr1);
2011  opts.range("sphericalBCr1", POSITIVE);
2012  opts.units("sphericalBCr1", N_ANGSTROM);
2013  opts.require("sphericalBC", "sphericalBCk1", "Force constant for first "
2014  "sphere potential (+ is an inward force, - outward)",
2015  &sphericalBCk1);
2016  opts.units("sphericalBCk1", N_KCAL);
2017  opts.optional("sphericalBC", "sphericalBCexp1", "Exponent for first "
2018  "sphere potential", &sphericalBCexp1, 2);
2019  opts.range("sphericalBCexp1", POSITIVE);
2020 
2021  opts.optional("sphericalBCr1", "sphericalBCr2", "Radius for second sphere "
2022  "potential", &sphericalBCr2);
2023  opts.range("sphericalBCr2", POSITIVE);
2024  opts.units("sphericalBCr2", N_ANGSTROM);
2025  opts.require("sphericalBCr2", "sphericalBCk2", "Force constant for second "
2026  "sphere potential (+ is an inward force, - outward)",
2027  &sphericalBCk2);
2028  opts.units("sphericalBCk2", N_KCAL);
2029  opts.optional("sphericalBCr2", "sphericalBCexp2", "Exponent for second "
2030  "sphere potential", &sphericalBCexp2, 2);
2031  opts.range("sphericalBCexp2", POSITIVE);
2032 
2034  opts.optionalB("main", "cylindricalBC", "Are cylindrical boundary counditions "
2035  "active?", &cylindricalBCOn, FALSE);
2036  opts.require("cylindricalBC", "cylindricalBCr1", "Radius for first cylinder "
2037  "potential", &cylindricalBCr1);
2038  opts.range("cylindricalBCr1", POSITIVE);
2039  opts.units("cylindricalBCr1", N_ANGSTROM);
2040  opts.require("cylindricalBC", "cylindricalBCk1", "Force constant for first "
2041  "cylinder potential (+ is an inward force, - outward)",
2042  &cylindricalBCk1);
2043  opts.units("cylindricalBCk1", N_KCAL);
2044  opts.optional("cylindricalBC", "cylindricalBCexp1", "Exponent for first "
2045  "cylinder potential", &cylindricalBCexp1, 2);
2046  opts.range("cylindricalBCexp1", POSITIVE);
2047 
2048 
2049 // additions beyond those already found in spherical parameters JJU
2050  opts.optional("cylindricalBC", "cylindricalBCAxis", "Cylinder axis (defaults to x)",
2051  PARSE_STRING);
2052  opts.require("cylindricalBC", "cylindricalBCCenter",
2053  "Center of cylindrical boundaries", &cylindricalCenter);
2054  opts.require ("cylindricalBC", "cylindricalBCl1", "Length of first cylinder",
2055  &cylindricalBCl1);
2056  opts.range("cylindricalBCl1", POSITIVE);
2057  opts.units("cylindricalBCl1", N_ANGSTROM);
2058  opts.optional ("cylindricalBCl1", "cylindricalBCl2", "Length of second cylinder",
2059  &cylindricalBCl2);
2060  opts.range ("cylindricalBCl2", POSITIVE);
2061  opts.units ("cylindricalBCl2", N_ANGSTROM);
2062 // end additions
2063 
2064  opts.optional("cylindricalBCr1", "cylindricalBCr2", "Radius for second cylinder "
2065  "potential", &cylindricalBCr2);
2066  opts.range("cylindricalBCr2", POSITIVE);
2067  opts.units("cylindricalBCr2", N_ANGSTROM);
2068  opts.require("cylindricalBCr2", "cylindricalBCk2", "Force constant for second "
2069  "cylinder potential (+ is an inward force, - outward)",
2070  &cylindricalBCk2);
2071  opts.units("cylindricalBCk2", N_KCAL);
2072  opts.optional("cylindricalBCr2", "cylindricalBCexp2", "Exponent for second "
2073  "cylinder potential", &cylindricalBCexp2, 2);
2074  opts.range("cylindricalBCexp2", POSITIVE);
2075 
2077  opts.optionalB("main", "eFieldOn", "Should an electric field be applied",
2078  &eFieldOn, FALSE);
2079  opts.optionalB("eFieldOn", "eFieldNormalized", "Is eField vector scaled by cell basis vectors?",
2080  &eFieldNormalized, FALSE);
2081  opts.require("eFieldOn", "eField", "Electric field vector", &eField);
2082  opts.optional("eFieldOn", "eFieldFreq", "Electric field frequency", &eFieldFreq);
2083  opts.optional("eFieldOn", "eFieldPhase", "Electric field phase", &eFieldPhase);
2084 
2086  opts.optionalB("main", "stirOn", "Should stirring torque be applied",
2087  &stirOn, FALSE);
2088  opts.optional("stirOn", "stirFilename", "PDB file with flags for "
2089  "stirred atoms (default is the PDB input file)",
2090  PARSE_STRING);
2091  opts.optional("stirOn", "stirredAtomsCol", "Column in the stirredAtomsFile "
2092  "containing the flags (nonzero means fixed);\n"
2093  "default is 'O'", PARSE_STRING);
2094  opts.require("stirOn", "stirStartingTheta", "Stir starting theta offset", &stirStartingTheta);
2095  opts.require("stirOn", "stirK", "Stir force harmonic spring constant", &stirK);
2096  //should make this optional, compute from firsttimestep * stirVel
2097  opts.require("stirOn", "stirVel", "Stir angular velocity (deg/timestep)", &stirVel);
2098  opts.require("stirOn", "stirAxis", "Stir axis (direction vector)", &stirAxis);
2099  opts.require("stirOn", "stirPivot", "Stir pivot point (coordinate)", &stirPivot);
2100 
2102  opts.optionalB("main", "extraBonds",
2103  "Should extra bonded forces be applied",
2104  &extraBondsOn, FALSE);
2105  opts.optional("extraBonds", "extraBondsFile",
2106  "file with list of extra bonds",
2107  PARSE_MULTIPLES);
2108  opts.optionalB("extraBonds", "extraBondsCosAngles",
2109  "Should extra angles be cosine-based to match ancient bug",
2110  &extraBondsCosAngles, TRUE);
2111 
2112 }
2113 
2114 void SimParameters::config_parser_misc(ParseOptions &opts) {
2115 
2117  opts.optional("main", "ldBalancer", "Load balancer",
2118  loadBalancer);
2119  opts.optional("main", "ldbStrategy", "Load balancing strategy",
2120  loadStrategy);
2121  opts.optional("main", "ldbPeriod", "steps between load balancing",
2122  &ldbPeriod);
2123  opts.range("ldbPeriod", POSITIVE);
2124  opts.optional("main", "firstLdbStep", "when to start load balancing",
2125  &firstLdbStep);
2126  opts.range("firstLdbStep", POSITIVE);
2127  opts.optional("main", "lastLdbStep", "when to stop load balancing",
2128  &lastLdbStep);
2129  opts.range("lastLdbStep", POSITIVE);
2130  opts.optional("main", "hybridGroupSize", "Hybrid load balancing group size",
2131  &hybridGroupSize);
2132  opts.optional("main", "ldbBackgroundScaling",
2133  "background load scaling", &ldbBackgroundScaling);
2134  opts.range("ldbBackgroundScaling", NOT_NEGATIVE);
2135  opts.optional("main", "ldbPMEBackgroundScaling",
2136  "PME node background load scaling", &ldbPMEBackgroundScaling);
2137  opts.range("ldbPMEBackgroundScaling", NOT_NEGATIVE);
2138  opts.optional("main", "ldbHomeBackgroundScaling",
2139  "home node background load scaling", &ldbHomeBackgroundScaling);
2140  opts.range("ldbHomeBackgroundScaling", NOT_NEGATIVE);
2141  opts.optional("main", "ldbRelativeGrainsize",
2142  "fraction of average load per compute", &ldbRelativeGrainsize, 0.);
2143  opts.range("ldbRelativeGrainsize", NOT_NEGATIVE);
2144 
2145  opts.optional("main", "traceStartStep", "when to start tracing", &traceStartStep);
2146  opts.range("traceStartStep", POSITIVE);
2147  opts.optional("main", "numTraceSteps", "the number of timesteps to be traced", &numTraceSteps);
2148  opts.range("numTraceSteps", POSITIVE);
2149 
2150 #ifdef MEASURE_NAMD_WITH_PAPI
2151  opts.optionalB("main", "papiMeasure", "whether use PAPI to measure performacne", &papiMeasure, FALSE);
2152  opts.optional("main", "papiMeasureStartStep", "when to measure performacne using PAPI", &papiMeasureStartStep);
2153  opts.range("papiMeasureStartStep", POSITIVE);
2154  opts.optional("main", "numPapiMeasureSteps", "the number of timesteps to be measured using PAPI", &numPapiMeasureSteps);
2155  opts.range("numPapiMeasureSteps", POSITIVE);
2156 #endif
2157 
2158  opts.optionalB("main", "outputMaps", "whether to dump compute map and patch map for analysis just before load balancing", &outputMaps, FALSE);
2159  opts.optionalB("main", "benchTimestep", "whether to do benchmarking timestep in which case final file output is disabled", &benchTimestep, FALSE);
2160  opts.optional("main", "useCkLoop", "whether to use CkLoop library to parallelize a loop in a function like OpenMP", &useCkLoop,
2161  #if CMK_SMP && USE_CKLOOP
2162  ( CkNumPes() < 2 * CkNumNodes() ? 0 : CKLOOP_CTRL_PME_FORWARDFFT ) );
2163  #else
2164  0);
2165  #endif
2166  opts.range("useCkLoop", NOT_NEGATIVE);
2167 
2168  opts.optionalB("main", "simulateInitialMapping", "whether to study the initial mapping scheme", &simulateInitialMapping, FALSE);
2169  opts.optional("main", "simulatedPEs", "the number of PEs to be used for studying initial mapping", &simulatedPEs);
2170  opts.range("simulatedPEs", POSITIVE);
2171  opts.optional("main", "simulatedNodeSize", "the node size to be used for studying initial mapping", &simulatedNodeSize);
2172  opts.range("simulatedNodeSize", POSITIVE);
2173  opts.optionalB("main", "disableTopology", "ignore torus information during patch placement", &disableTopology, FALSE);
2174  opts.optionalB("main", "verboseTopology", "print torus information during patch placement", &verboseTopology, FALSE);
2175 
2176  opts.optionalB("main", "ldbUnloadPME", "no load on PME nodes",
2177  &ldbUnloadPME, FALSE);
2178  opts.optionalB("main", "ldbUnloadZero", "no load on pe zero",
2179  &ldbUnloadZero, FALSE);
2180  opts.optionalB("main", "ldbUnloadOne", "no load on pe one",
2181  &ldbUnloadOne, FALSE);
2182  opts.optionalB("main", "ldbUnloadOutputPEs", "no load on output PEs",
2183  &ldbUnloadOutputPEs, FALSE);
2184  opts.optionalB("main", "noPatchesOnZero", "no patches on pe zero",
2185  &noPatchesOnZero, FALSE);
2186  opts.optionalB("main", "noPatchesOnOutputPEs", "no patches on Output PEs",
2187  &noPatchesOnOutputPEs, FALSE);
2188  opts.optionalB("main", "noPatchesOnOne", "no patches on pe one",
2189  &noPatchesOnOne, FALSE);
2190  opts.optionalB("main", "useCompressedPsf", "The structure file psf is in the compressed format",
2191  &useCompressedPsf, FALSE);
2192  opts.optionalB("main", "genCompressedPsf", "Generate the compressed version of the psf file",
2193  &genCompressedPsf, FALSE);
2194  opts.optionalB("main", "usePluginIO", "Use the plugin I/O to load the molecule system",
2195  &usePluginIO, FALSE);
2196  opts.optionalB("main", "mallocTest", "test how much memory all PEs can allocate",
2197  &mallocTest, FALSE);
2198  opts.optionalB("main", "printExclusions", "print exclusion lists to stdout",
2199  &printExclusions, FALSE);
2200  opts.optional("main", "proxySendSpanningTree", "using spanning tree to send proxies",
2201  &proxySendSpanningTree, -1);
2202  opts.optional("main", "proxyRecvSpanningTree", "using spanning tree to receive proxies",
2203  &proxyRecvSpanningTree, 0); // default off due to memory leak -1);
2204  opts.optional("main", "proxyTreeBranchFactor", "the branch factor when building a spanning tree",
2205  &proxyTreeBranchFactor, 0); // actual default in ProxyMgr.C
2206  opts.optionalB("main", "twoAwayX", "half-size patches in 1st dimension",
2207  &twoAwayX, -1);
2208  opts.optionalB("main", "twoAwayY", "half-size patches in 2nd dimension",
2209  &twoAwayY, -1);
2210  opts.optionalB("main", "twoAwayZ", "half-size patches in 3rd dimension",
2211  &twoAwayZ, -1);
2212  opts.optional("main", "maxPatches", "maximum patch count", &maxPatches, -1);
2213 
2215  opts.optional("main", "firsttimestep", "Timestep to start simulation at",
2216  &firstTimestep, 0);
2217  opts.range("firsttimestep", NOT_NEGATIVE);
2218 
2220  opts.optionalB("main", "test", "Perform self-tests rather than simulation",
2221  &testOn, FALSE);
2222  opts.optionalB("main", "commOnly", "Do not evaluate forces or integrate",
2223  &commOnly, FALSE);
2224 
2225  opts.optionalB("main", "statsOn", "counters in machine layer",
2226  &statsOn, FALSE);
2228  opts.optionalB("main", "hbonds", "Use explicit hydrogen bond term",
2229  &HydrogenBonds, FALSE);
2230  opts.optionalB("hbonds","hbAntecedents","Include Antecedent in hbond term",
2231  &useAntecedent, TRUE);
2232  opts.optional("hbonds","hbAAexp","Hbond AA-A-H angle cos exponential",
2233  &aaAngleExp, 2);
2234  opts.optional("hbonds","hbHAexp","Hbond D-H-A angle cos exponential",
2235  &haAngleExp, 4);
2236  opts.optional("hbonds","hbDistAexp","Hbond A-D dist attractive exponential",
2237  &distAttExp, 4);
2238  opts.optional("hbonds","hbDistRexp","Hbond A-D dist repulstive exponential",
2239  &distRepExp, 6);
2240  opts.optional("hbonds","hbCutoffAngle","Hbond D-H-A cutoff angle",
2241  &dhaCutoffAngle, 100.0);
2242  opts.range("hbCutoffAngle", NOT_NEGATIVE);
2243  opts.optional("hbonds","hbOnAngle","Hbond D-H-A switch function on angle",
2244  &dhaOnAngle, 60.0);
2245  opts.range("hbOnAngle", NOT_NEGATIVE);
2246  opts.optional("hbonds","hbOffAngle","Hbond D-H-A switch function off angle",
2247  &dhaOffAngle, 80.0);
2248  opts.range("hbOffAngle", NOT_NEGATIVE);
2249  opts.optional("hbonds","hbCutoffDist","Hbond A-D cutoff distance",
2250  &daCutoffDist, 7.5);
2251  opts.range("hbCutoffDist", POSITIVE);
2252  opts.units("hbCutoffDist", N_ANGSTROM);
2253  opts.optional("hbonds","hbOnDist","Hbond A-D switch function on distance",
2254  &daOnDist, 5.5);
2255  opts.range("hbOnDist", POSITIVE);
2256  opts.units("hbOnDist", N_ANGSTROM);
2257  opts.optional("hbonds","hbOffDist","Hbond A-D switch function off distance",
2258  &daOffDist, 6.5);
2259  opts.range("hbOffDist", POSITIVE);
2260  opts.units("hbOffDist", N_ANGSTROM);
2261 
2262  // IMD options
2263  opts.optionalB("main","IMDon","Connect using IMD?",&IMDon, FALSE);
2264  opts.require("IMDon","IMDport", "Port to which to bind", &IMDport);
2265  opts.range("IMDport",POSITIVE);
2266  opts.require("IMDon","IMDfreq", "Frequency at which to report", &IMDfreq);
2267  opts.range("IMDfreq",POSITIVE);
2268  opts.optionalB("IMDon","IMDwait","Pause until IMD connection?",&IMDwait,
2269  FALSE);
2270  opts.optionalB("IMDon","IMDignore","Ignore any user input?",&IMDignore,
2271  FALSE);
2272  opts.optionalB("IMDon","IMDignoreForces","Ignore forces ONLY?",&IMDignoreForces,
2273  FALSE);
2274  // Maximum Partition options
2275  opts.optional("ldBalancer", "maxSelfPart",
2276  "maximum number of self partitions in one patch", &maxSelfPart, 20);
2277  opts.range("maxSelfPart",POSITIVE);
2278  opts.optional("ldBalancer", "maxPairPart",
2279  "maximum number of pair partitions in one patch", &maxPairPart, 8);
2280  opts.range("maxPairPart",POSITIVE);
2281  opts.optional("ldBalancer", "numAtomsSelf",
2282  "maximum number of atoms in one self compute distribution",
2283  &numAtomsSelf, 154);
2284  opts.range("numAtomsSelf",NOT_NEGATIVE);
2285 
2286  opts.optional("ldBalancer", "numAtomsSelf2",
2287  "maximum number of atoms in one self compute distribution",
2288  &numAtomsSelf2, 154);
2289  opts.range("numAtomsSelf2",NOT_NEGATIVE);
2290 
2291  opts.optional("ldBalancer", "numAtomsPair",
2292  "maximum number of atoms in one pair compute distribution",
2293  &numAtomsPair, 318);
2294  opts.range("numAtomsPair",NOT_NEGATIVE);
2295  opts.optional("ldBalancer", "numAtomsPair2",
2296  "maximum number of atoms in one pair compute distribution",
2297  &numAtomsPair2, 637);
2298  opts.range("numAtomsPair2",NOT_NEGATIVE);
2299  opts.optional("main", "minAtomsPerPatch",
2300  "minimum average atoms per patch",
2301  &minAtomsPerPatch, 40);
2302  opts.range("minAtomsPerPatch",NOT_NEGATIVE);
2303 
2304  // Set number added to patch atom count during initial node assignment
2305  opts.optional("main", "emptyPatchLoad",
2306  "load generated by empty patch, in atoms",
2307  &emptyPatchLoad, 40);
2308  opts.range("emptyPatchLoad",POSITIVE);
2309 
2310  // Maximum exclusion flags per atom
2311  opts.optional("main", "maxExclusionFlags",
2312  "maximum number of exclusion flags per atom", &maxExclusionFlags, 256);
2313  opts.range("maxExclusionFlags",POSITIVE);
2314 
2315  // Bonded interactions on GPU
2316  opts.optional("main", "bondedCUDA", "Bitmask for calculating bonded interactions on GPU", &bondedCUDA, 255);
2317 
2318  // Automatically disable individual CUDA kernels that are
2319  // incompatible with simulation options.
2320  // Set FALSE to manually control kernel use for development.
2321  opts.optionalB("main", "useCUDAdisable", "Disable kernels to maintain feature compatibility with CUDA", &useCUDAdisable, TRUE);
2322 
2323  // MIC specific parameters
2324  opts.optional("main", "mic_unloadMICPEs", "Indicates whether or not the load balancer should unload PEs driving Xeon Phi cards", &mic_unloadMICPEs, 1);
2325  opts.optional("main", "mic_singleKernel", "Set to non-zero to have all MIC work to be placed in a single kernel", &mic_singleKernel, 1);
2326  opts.optional("main", "mic_deviceThreshold", "Threshold to use for directing computes to Xeon Phi devices", &mic_deviceThreshold, -1);
2327  opts.optional("main", "mic_hostSplit", "DMK - reserved", &mic_hostSplit, -1);
2328  opts.optional("main", "mic_numParts_self_p1", "MIC-Specific NumParts SELF Parameter 1", &mic_numParts_self_p1, -1);
2329  opts.optional("main", "mic_numParts_pair_p1", "MIC-Specific NumParts PAIR Parameter 1", &mic_numParts_pair_p1, -1);
2330  opts.optional("main", "mic_numParts_pair_p2", "MIC-Specific NumParts PAIR Parameter 2", &mic_numParts_pair_p2, -1);
2331  opts.range("mic_unloadMICPEs", NOT_NEGATIVE);
2332  opts.range("mic_singleKernel", NOT_NEGATIVE);
2333 }
2334 
2335 void SimParameters::readExtendedSystem(const char *filename, Lattice *latptr) {
2336 
2337  if ( ! latptr ) {
2338  iout << iINFO << "EXTENDED SYSTEM FILE " << filename << "\n" << endi;
2339  }
2340 
2341  ifstream xscFile(filename);
2342  if ( ! xscFile ) NAMD_die("Unable to open extended system file.\n");
2343 
2344  char labels[1024];
2345  do {
2346  if ( ! xscFile ) NAMD_die("Error reading extended system file.\n");
2347  xscFile.getline(labels,1023);
2348  } while ( strncmp(labels,"#$LABELS ",9) );
2349 
2350  int a_x, a_y, a_z, b_x, b_y, b_z, c_x, c_y, c_z;
2351  a_x = a_y = a_z = b_x = b_y = b_z = c_x = c_y = c_z = -1;
2352  int o_x, o_y, o_z, s_u, s_v, s_w, s_x, s_y, s_z;
2353  o_x = o_y = o_z = s_u = s_v = s_w = s_x = s_y = s_z = -1;
2354 
2355  int pos = 0;
2356  char *l_i = labels + 8;
2357  while ( *l_i ) {
2358  if ( *l_i == ' ' ) { ++l_i; continue; }
2359  char *l_i2;
2360  for ( l_i2 = l_i; *l_i2 && *l_i2 != ' '; ++l_i2 );
2361  if ( (l_i2 - l_i) == 3 && (l_i[1] == '_') ) {
2362  if (l_i[0] == 'a' && l_i[2] == 'x') a_x = pos;
2363  if (l_i[0] == 'a' && l_i[2] == 'y') a_y = pos;
2364  if (l_i[0] == 'a' && l_i[2] == 'z') a_z = pos;
2365  if (l_i[0] == 'b' && l_i[2] == 'x') b_x = pos;
2366  if (l_i[0] == 'b' && l_i[2] == 'y') b_y = pos;
2367  if (l_i[0] == 'b' && l_i[2] == 'z') b_z = pos;
2368  if (l_i[0] == 'c' && l_i[2] == 'x') c_x = pos;
2369  if (l_i[0] == 'c' && l_i[2] == 'y') c_y = pos;
2370  if (l_i[0] == 'c' && l_i[2] == 'z') c_z = pos;
2371  if (l_i[0] == 'o' && l_i[2] == 'x') o_x = pos;
2372  if (l_i[0] == 'o' && l_i[2] == 'y') o_y = pos;
2373  if (l_i[0] == 'o' && l_i[2] == 'z') o_z = pos;
2374  if (l_i[0] == 's' && l_i[2] == 'u') s_u = pos;
2375  if (l_i[0] == 's' && l_i[2] == 'v') s_v = pos;
2376  if (l_i[0] == 's' && l_i[2] == 'w') s_w = pos;
2377  if (l_i[0] == 's' && l_i[2] == 'x') s_x = pos;
2378  if (l_i[0] == 's' && l_i[2] == 'y') s_y = pos;
2379  if (l_i[0] == 's' && l_i[2] == 'z') s_z = pos;
2380  }
2381  ++pos;
2382  l_i = l_i2;
2383  }
2384  int numpos = pos;
2385 
2386  for ( pos = 0; pos < numpos; ++pos ) {
2387  double tmp;
2388  xscFile >> tmp;
2389  if ( ! xscFile ) NAMD_die("Error reading extended system file.\n");
2390  if ( pos == a_x ) cellBasisVector1.x = tmp;
2391  if ( pos == a_y ) cellBasisVector1.y = tmp;
2392  if ( pos == a_z ) cellBasisVector1.z = tmp;
2393  if ( pos == b_x ) cellBasisVector2.x = tmp;
2394  if ( pos == b_y ) cellBasisVector2.y = tmp;
2395  if ( pos == b_z ) cellBasisVector2.z = tmp;
2396  if ( pos == c_x ) cellBasisVector3.x = tmp;
2397  if ( pos == c_y ) cellBasisVector3.y = tmp;
2398  if ( pos == c_z ) cellBasisVector3.z = tmp;
2399  if ( pos == o_x ) cellOrigin.x = tmp;
2400  if ( pos == o_y ) cellOrigin.y = tmp;
2401  if ( pos == o_z ) cellOrigin.z = tmp;
2402  if ( pos == s_u ) strainRate2.x = tmp;
2403  if ( pos == s_v ) strainRate2.y = tmp;
2404  if ( pos == s_w ) strainRate2.z = tmp;
2405  if ( pos == s_x ) strainRate.x = tmp;
2406  if ( pos == s_y ) strainRate.y = tmp;
2407  if ( pos == s_z ) strainRate.z = tmp;
2408  }
2409 
2410  if ( latptr ) {
2411  Lattice test;
2412  test.set(cellBasisVector1,cellBasisVector2,cellBasisVector3,cellOrigin);
2413 
2414  if ( test.a_p() && ! lattice.a_p() ) {
2415  NAMD_die("cellBasisVector1 added during atom reinitialization");
2416  }
2417  if ( lattice.a_p() && ! test.a_p() ) {
2418  NAMD_die("cellBasisVector1 dropped during atom reinitialization");
2419  }
2420  if ( test.b_p() && ! lattice.b_p() ) {
2421  NAMD_die("cellBasisVector2 added during atom reinitialization");
2422  }
2423  if ( lattice.b_p() && ! test.b_p() ) {
2424  NAMD_die("cellBasisVector2 dropped during atom reinitialization");
2425  }
2426  if ( test.c_p() && ! lattice.c_p() ) {
2427  NAMD_die("cellBasisVector3 added during atom reinitialization");
2428  }
2429  if ( lattice.c_p() && ! test.c_p() ) {
2430  NAMD_die("cellBasisVector3 dropped during atom reinitialization");
2431  }
2432 
2433  latptr->set(cellBasisVector1,cellBasisVector2,cellBasisVector3,cellOrigin);
2434  }
2435 
2436 }
2437 
2438 #ifdef MEM_OPT_VERSION
2439 //This global var is defined in mainfunc.C
2440 extern char *gWorkDir;
2441 #endif
2442 
2443 void SimParameters::check_config(ParseOptions &opts, ConfigList *config, char *&cwd) {
2444 
2445  int len; // String length
2446  StringList *current; // Pointer to config option list
2447 
2448 #ifdef MEM_OPT_VERSION
2449  char *namdWorkDir = NULL;
2450 #endif
2451 
2452  if ( opts.defined("obsolete") ) {
2453  iout << iWARN <<
2454  "\"obsolete\" defined, silently ignoring obsolete options\n" << endi;
2455  }
2456 
2457  // Take care of cwd processing
2458  if (opts.defined("cwd"))
2459  {
2460  // First allocate and get the cwd value
2461  current = config->find("cwd");
2462 
2463  len = strlen(current->data);
2464 
2465  if ( CHDIR(current->data) )
2466  {
2467  NAMD_die("chdir() to given cwd failed!");
2468  } else {
2469  iout << iINFO << "Changed directory to " << current->data << "\n" << endi;
2470  }
2471 
2472  if (current->data[len-1] != PATHSEP)
2473  len++;
2474 
2475  cwd = new char[len+1];
2476 
2477  strcpy(cwd, current->data);
2478 
2479  if (current->data[strlen(current->data)-1] != PATHSEP)
2480  strcat(cwd, PATHSEPSTR);
2481  }
2482 
2483 #ifdef MEM_OPT_VERSION
2484  if(cwd!=NULL)namdWorkDir = cwd;
2485  else namdWorkDir = gWorkDir;
2486  int dirlen = strlen(namdWorkDir);
2487  //only support the path representation on UNIX-like platforms
2488  char *tmpDir;
2489  if(namdWorkDir[dirlen-1]=='/'){
2490  tmpDir = new char[dirlen+1];
2491  tmpDir[dirlen] = 0;
2492  }else{
2493  tmpDir = new char[dirlen+2];
2494  tmpDir[dirlen]='/';
2495  tmpDir[dirlen+1]=0;
2496  }
2497  memcpy(tmpDir, namdWorkDir, dirlen);
2498  namdWorkDir = tmpDir;
2499  //finished recording the per atom files, free the space for gWorkDir
2500  delete [] gWorkDir;
2501 #endif
2502 
2503 
2504  // Don't try to specify coordinates with pluginIO
2505  if ( usePluginIO && opts.defined("coordinates") ) {
2506  NAMD_die("Separate coordinates file not allowed with plugin IO, coordinates will be taken from structure file.");
2507  }
2508 
2509  // If it's not AMBER||GROMACS, then "coordinates", "structure"
2510  // and "parameters" must be specified.
2511  if (!amberOn && !gromacsOn) {
2512 #ifndef MEM_OPT_VERSION
2513  if (useCompressedPsf)
2514  NAMD_die("useCompressedPsf requires memory-optimized build!");
2515  if (!usePluginIO && !genCompressedPsf && !opts.defined("coordinates"))
2516  NAMD_die("coordinates not found in the configuration file!");
2517 #else
2518  if(!usePluginIO && !opts.defined("bincoordinates")) {
2519  NAMD_die("bincoordinates not found in the configuration file for the memory optimized version!");
2520  }
2521  if(!usePluginIO && opts.defined("coordinates")) {
2522  NAMD_die("coordinates not allowed in the configuration file for the memory optimized version!");
2523  }
2524 #endif
2525  if (!opts.defined("structure"))
2526  NAMD_die("structure not found in the configuration file!");
2527  if (!opts.defined("parameters"))
2528  NAMD_die("parameters not found in the configuration file!");
2529  }
2530 
2531  // In any case, there should be either "coordinates" or
2532  // "ambercoor", but not both
2533  if (opts.defined("coordinates") && opts.defined("ambercoor"))
2534  NAMD_die("Cannot specify both coordinates and ambercoor!");
2535 #ifndef MEM_OPT_VERSION
2536  if (!genCompressedPsf && !opts.defined("coordinates") && !opts.defined("ambercoor")
2537  && !opts.defined("grocoorfile") && !usePluginIO)
2538  NAMD_die("Coordinate file not found!");
2539 #endif
2540 
2541  // Make sure that both a temperature and a velocity PDB were
2542  // specified
2543  if (opts.defined("temperature") &&
2544  (opts.defined("velocities") || opts.defined("binvelocities")) )
2545  {
2546  NAMD_die("Cannot specify both an initial temperature and a velocity file");
2547  }
2548 
2549 #ifdef MEM_OPT_VERSION
2550 //record the absolute file name for binAtomFile, binCoorFile and binVelFile etc.
2551  binAtomFile = NULL;
2552  binCoorFile = NULL;
2553  binVelFile = NULL;
2554  binRefFile = NULL;
2555 
2556  char *curfile = NULL;
2557  dirlen = strlen(namdWorkDir);
2558  current = config->find("structure");;
2559  curfile = current->data;
2560  int filelen = strlen(curfile);
2561  if(*curfile == '/' || *curfile=='~') {
2562  //check whether it is an absolute path
2563  //WARNING: Only works on Unix-like platforms!
2564  //Needs to fix on Windows platform.
2565  //-Chao Mei
2566  //adding 5 because of ".bin"+"\0"
2567  binAtomFile = new char[filelen+5];
2568  memcpy(binAtomFile, curfile, filelen);
2569  memcpy(binAtomFile+filelen, ".bin", 4);
2570  binAtomFile[filelen+4] = 0;
2571  }else{
2572  binAtomFile = new char[dirlen+filelen+5];
2573  memcpy(binAtomFile, namdWorkDir, dirlen);
2574  memcpy(binAtomFile+dirlen, curfile, filelen);
2575  memcpy(binAtomFile+dirlen+filelen, ".bin", 4);
2576  binAtomFile[dirlen+filelen+4] = 0;
2577  }
2578 
2579  current = config->find("bincoordinates");
2580  curfile = current->data;
2581  filelen = strlen(curfile);
2582  if(*curfile == '/' || *curfile=='~') {
2583  binCoorFile = new char[filelen+1];
2584  memcpy(binCoorFile, curfile, filelen);
2585  binCoorFile[filelen] = 0;
2586  }else{
2587  binCoorFile = new char[dirlen+filelen+1];
2588  memcpy(binCoorFile, namdWorkDir, dirlen);
2589  memcpy(binCoorFile+dirlen, curfile, filelen);
2590  binCoorFile[dirlen+filelen] = 0;
2591  }
2592 
2593  if(opts.defined("binvelocities")){
2594  current = config->find("binvelocities");
2595  curfile = current->data;
2596  filelen = strlen(curfile);
2597  if(*curfile == '/' || *curfile=='~') {
2598  binVelFile = new char[filelen+1];
2599  memcpy(binVelFile, curfile, filelen);
2600  binVelFile[filelen] = 0;
2601  }else{
2602  binVelFile = new char[dirlen+filelen+1];
2603  memcpy(binVelFile, namdWorkDir, dirlen);
2604  memcpy(binVelFile+dirlen, curfile, filelen);
2605  binVelFile[dirlen+filelen] = 0;
2606  }
2607  }
2608 
2609  if(opts.defined("binrefcoords")){
2610  current = config->find("binrefcoords");
2611  curfile = current->data;
2612  filelen = strlen(curfile);
2613  if(*curfile == '/' || *curfile=='~') {
2614  binRefFile = new char[filelen+1];
2615  memcpy(binRefFile, curfile, filelen);
2616  binRefFile[filelen] = 0;
2617  }else{
2618  binRefFile = new char[dirlen+filelen+1];
2619  memcpy(binRefFile, namdWorkDir, dirlen);
2620  memcpy(binRefFile+dirlen, curfile, filelen);
2621  binRefFile[dirlen+filelen] = 0;
2622  }
2623  }
2624 
2625  //deal with output file name to make it absolute path for parallel output
2626  if(outputFilename[0] != '/' && outputFilename[0]!='~') {
2627  filelen = strlen(outputFilename);
2628  char *tmpout = new char[filelen];
2629  memcpy(tmpout, outputFilename, filelen);
2630  CmiAssert(filelen+dirlen <= 120); //leave 8 chars for file suffix
2631  memcpy(outputFilename, namdWorkDir, dirlen);
2632  memcpy(outputFilename+dirlen, tmpout, filelen);
2633  outputFilename[filelen+dirlen] = 0;
2634  delete [] tmpout;
2635  }
2636 
2637  if ( dcdFrequency && opts.defined("dcdfile") &&
2638  dcdFilename[0] != '/' && dcdFilename[0]!='~' ) {
2639  filelen = strlen(dcdFilename);
2640  char *tmpout = new char[filelen];
2641  memcpy(tmpout, dcdFilename, filelen);
2642  CmiAssert(filelen+dirlen <= 120); //leave 8 chars for file suffix
2643  memcpy(dcdFilename, namdWorkDir, dirlen);
2644  memcpy(dcdFilename+dirlen, tmpout, filelen);
2645  dcdFilename[filelen+dirlen] = 0;
2646  delete [] tmpout;
2647  }
2648 
2649  if ( velDcdFrequency && opts.defined("veldcdfile") &&
2650  velDcdFilename[0] != '/' && velDcdFilename[0]!='~' ) {
2651  filelen = strlen(velDcdFilename);
2652  char *tmpout = new char[filelen];
2653  memcpy(tmpout, velDcdFilename, filelen);
2654  CmiAssert(filelen+dirlen <= 120); //leave 8 chars for file suffix
2655  memcpy(velDcdFilename, namdWorkDir, dirlen);
2656  memcpy(velDcdFilename+dirlen, tmpout, filelen);
2657  velDcdFilename[filelen+dirlen] = 0;
2658  delete [] tmpout;
2659  }
2660 
2661  if ( forceDcdFrequency && opts.defined("forcedcdfile") &&
2662  forceDcdFilename[0] != '/' && forceDcdFilename[0]!='~' ) {
2663  filelen = strlen(forceDcdFilename);
2664  char *tmpout = new char[filelen];
2665  memcpy(tmpout, forceDcdFilename, filelen);
2666  CmiAssert(filelen+dirlen <= 120); //leave 8 chars for file suffix
2667  memcpy(forceDcdFilename, namdWorkDir, dirlen);
2668  memcpy(forceDcdFilename+dirlen, tmpout, filelen);
2669  forceDcdFilename[filelen+dirlen] = 0;
2670  delete [] tmpout;
2671  }
2672 
2673  if ( restartFrequency && opts.defined("restartname") &&
2674  restartFilename[0] != '/' && restartFilename[0]!='~' ) {
2675  filelen = strlen(restartFilename);
2676  char *tmpout = new char[filelen];
2677  memcpy(tmpout, restartFilename, filelen);
2678  CmiAssert(filelen+dirlen <= 120); //leave 8 chars for file suffix
2679  memcpy(restartFilename, namdWorkDir, dirlen);
2680  memcpy(restartFilename+dirlen, tmpout, filelen);
2681  restartFilename[filelen+dirlen] = 0;
2682  delete [] tmpout;
2683  }
2684 
2685  delete [] namdWorkDir;
2686 
2687  if (opts.defined("numinputprocs")) {
2688  if(numinputprocs > CkNumPes()) {
2689  iout << iWARN << "The number of input processors exceeds the total number of processors. Resetting to half of the number of total processors.\n" << endi;
2690  numinputprocs = (CkNumPes()>>1)+(CkNumPes()&1);
2691  }
2692  }
2693 
2694  if (opts.defined("numoutputprocs")) {
2695  if(numoutputprocs > CkNumPes()) {
2696  iout << iWARN << "The number of output processors exceeds the total number of processors. Resetting to half of the number of total processors.\n" << endi;
2697  numoutputprocs = (CkNumPes()>>1)+(CkNumPes()&1);
2698  }
2699  }
2700 
2701 #ifndef OUTPUT_SINGLE_FILE
2702 #error OUTPUT_SINGLE_FILE not defined!
2703 #endif
2704 
2705  #if !OUTPUT_SINGLE_FILE
2706  //create directories for multi-file output scheme
2707  create_output_directories("coor");
2708  create_output_directories("vel");
2709  if(dcdFrequency) {
2710  create_output_directories("dcd");
2711  if(opts.defined("dcdfile")){
2712  iout << iWARN << "The dcd file output has been changed to directory: " << outputFilename << ".\n" << endi;
2713  }
2714  }
2715  if (velDcdFrequency) {
2716  create_output_directories("veldcd");
2717  if(opts.defined("veldcdfile")){
2718  iout << iWARN << "The veldcd file output has been changed to directory: " << outputFilename << ".\n" << endi;
2719  }
2720  }
2721  if (forceDcdFrequency) {
2722  create_output_directories("forcedcd");
2723  if(opts.defined("forcedcdfile")){
2724  iout << iWARN << "The forcedcd file output has been changed to directory: " << outputFilename << ".\n" << endi;
2725  }
2726  }
2727  #endif
2728 #endif
2729 
2730  if (! opts.defined("auxFile")) {
2731  strcpy(auxFilename,outputFilename);
2732  strcat(auxFilename,".aux");
2733  }
2734 
2735  // Check for frequencies
2736  if (dcdFrequency) {
2737  if (! opts.defined("dcdfile")) {
2738  strcpy(dcdFilename,outputFilename);
2739  strcat(dcdFilename,".dcd");
2740  }
2741  } else {
2742  dcdFilename[0] = STRINGNULL;
2743  }
2744 
2745  if (velDcdFrequency) {
2746  if (! opts.defined("veldcdfile")) {
2747  strcpy(velDcdFilename,outputFilename);
2748  strcat(velDcdFilename,".veldcd");
2749  }
2750  } else {
2751  velDcdFilename[0] = STRINGNULL;
2752  }
2753 
2754  if (forceDcdFrequency) {
2755  if (! opts.defined("forcedcdfile")) {
2756  strcpy(forceDcdFilename,outputFilename);
2757  strcat(forceDcdFilename,".forcedcd");
2758  }
2759  } else {
2760  forceDcdFilename[0] = STRINGNULL;
2761  }
2762 
2763  if (xstFrequency) {
2764  if (! opts.defined("xstfile")) {
2765  strcpy(xstFilename,outputFilename);
2766  strcat(xstFilename,".xst");
2767  }
2768  } else {
2769  xstFilename[0] = STRINGNULL;
2770  }
2771 
2772  if (restartFrequency) {
2773  if (! opts.defined("restartname")) {
2774  strcpy(restartFilename,outputFilename);
2775  if ( ! restartSave ) strcat(restartFilename,".restart");
2776  }
2777  } else {
2778  restartFilename[0] = STRINGNULL;
2779  restartSave = FALSE;
2780  binaryRestart = FALSE;
2781  }
2782 
2783  if (storeComputeMap || loadComputeMap) {
2784  if (! opts.defined("computeMapFile")) {
2785  strcpy(computeMapFilename,"computeMapFile");
2786  strcat(computeMapFilename,".txt");
2787  }
2788  }
2789 
2790 
2791  if (!amberOn)
2792  { //****** BEGIN CHARMM/XPLOR type changes
2794  if (!paraTypeXplorOn && !paraTypeCharmmOn)
2795  {
2796  paraTypeXplorOn = TRUE;
2797  }
2799  if (paraTypeXplorOn && paraTypeCharmmOn)
2800  {
2801  NAMD_die("Please specify either XPLOR or CHARMM format for parameters!");
2802  }
2803  //****** END CHARMM/XPLOR type changes
2804  }
2805 
2806 
2807  // If minimization isn't on, must have a temp or velocity
2808  if (!(minimizeOn||minimizeCGOn) && !opts.defined("temperature") &&
2809  !opts.defined("velocities") && !opts.defined("binvelocities") )
2810  {
2811  NAMD_die("Must have either an initial temperature or a velocity file");
2812  }
2813 
2814  if (minimizeOn||minimizeCGOn) { initialTemp = 0.0; }
2815  if (opts.defined("velocities") || opts.defined("binvelocities") )
2816  {
2817  initialTemp = -1.0;
2818  }
2819 
2821 
2822  if ( opts.defined("extendedSystem") ) readExtendedSystem(config->find("extendedSystem")->data);
2823 
2824 #ifdef MEM_OPT_VERSION
2825  if ( LJcorrection ) {
2826  NAMD_die("LJ tail corrections not yet available for memory optimized builds");
2827  }
2828 #endif
2829 
2830  if ( LJcorrection && ! cellBasisVector3.length2() ) {
2831  NAMD_die("Can't use LJ tail corrections without periodic boundary conditions!");
2832  }
2833 
2834  if ( cellBasisVector3.length2() && ! cellBasisVector2.length2() ) {
2835  NAMD_die("Used cellBasisVector3 without cellBasisVector2!");
2836  }
2837 
2838  if ( cellBasisVector2.length2() && ! cellBasisVector1.length2() ) {
2839  NAMD_die("Used cellBasisVector2 without cellBasisVector1!");
2840  }
2841 
2842  if ( cellOrigin.length2() && ! cellBasisVector1.length2() ) {
2843  NAMD_die("Used cellOrigin without cellBasisVector1!");
2844  }
2845 
2846  lattice.set(cellBasisVector1,cellBasisVector2,cellBasisVector3,cellOrigin);
2847 
2848  if (! opts.defined("DCDunitcell")) {
2849  dcdUnitCell = lattice.a_p() && lattice.b_p() && lattice.c_p();
2850  }
2851 
2852  char s[129];
2853 
2855  if ( ! opts.defined("cylindricalBCAxis") )
2856  {
2857  cylindricalBCAxis = 'x';
2858  }
2859  else
2860  {
2861  opts.get("cylindricalBCAxis", s);
2862 
2863  if (!strcasecmp(s, "x"))
2864  {
2865  cylindricalBCAxis = 'x';
2866  }
2867  else if (!strcasecmp(s, "y"))
2868  {
2869  cylindricalBCAxis = 'y';
2870  }
2871  else if (!strcasecmp(s, "z"))
2872  {
2873  cylindricalBCAxis = 'z';
2874  }
2875  else
2876  {
2877  char err_msg[128];
2878 
2879  sprintf(err_msg, "Illegal value '%s' for 'cylindricalBCAxis' in configuration file", s);
2880  NAMD_die(err_msg);
2881  }
2882  }
2883 
2884  if (!opts.defined("splitPatch"))
2885  {
2886  splitPatch = SPLIT_PATCH_HYDROGEN;
2887  }
2888  else
2889  {
2890  opts.get("splitPatch", s);
2891  if (!strcasecmp(s, "position"))
2892  splitPatch = SPLIT_PATCH_POSITION;
2893  else if (!strcasecmp(s,"hydrogen"))
2894  splitPatch = SPLIT_PATCH_HYDROGEN;
2895  else
2896  {
2897  char err_msg[129];
2898  sprintf(err_msg,
2899  "Illegal value '%s' for 'splitPatch' in configuration file",
2900  s);
2901  NAMD_die(err_msg);
2902  }
2903  }
2904 
2906  opts.get("exclude", s);
2907 
2908  if (!strcasecmp(s, "none"))
2909  {
2910  exclude = NONE;
2911  splitPatch = SPLIT_PATCH_POSITION;
2912  }
2913  else if (!strcasecmp(s, "1-2"))
2914  {
2915  exclude = ONETWO;
2916  splitPatch = SPLIT_PATCH_POSITION;
2917  }
2918  else if (!strcasecmp(s, "1-3"))
2919  {
2920  exclude = ONETHREE;
2921  }
2922  else if (!strcasecmp(s, "1-4"))
2923  {
2924  exclude = ONEFOUR;
2925  }
2926  else if (!strcasecmp(s, "scaled1-4"))
2927  {
2928  exclude = SCALED14;
2929  }
2930  else
2931  {
2932  char err_msg[128];
2933 
2934  sprintf(err_msg, "Illegal value '%s' for 'exclude' in configuration file",
2935  s);
2936  NAMD_die(err_msg);
2937  }
2938 
2939  if (scale14 != 1.0 && exclude != SCALED14)
2940  {
2941  iout << iWARN << "Exclude is not scaled1-4; 1-4scaling ignored.\n" << endi;
2942  }
2943 
2944  // water model stuff
2945  if (!opts.defined("waterModel")) {
2946  watmodel = WAT_TIP3;
2947  } else {
2948  opts.get("waterModel", s);
2949  if (!strncasecmp(s, "tip4", 4)) {
2950  iout << iINFO << "Using TIP4P water model.\n" << endi;
2951  watmodel = WAT_TIP4;
2952  } else if (!strncasecmp(s, "tip3", 4)) {
2953  iout << iINFO << "Using TIP3P water model.\n" << endi;
2954  watmodel = WAT_TIP3;
2955  } else if (!strncasecmp(s, "swm4", 4)) {
2956  iout << iINFO << "Using SWM4-DP water model.\n" << endi;
2957  watmodel = WAT_SWM4;
2958  } else {
2959  char err_msg[128];
2960  sprintf(err_msg,
2961  "Illegal value %s for 'waterModel' in configuration file", s);
2962  NAMD_die(err_msg);
2963  }
2964  }
2965  if (watmodel == WAT_SWM4 && !drudeOn) {
2966  NAMD_die("Must have 'drudeOn' enabled to use SWM4-DP water model.");
2967  }
2968  if (drudeOn && watmodel != WAT_SWM4) {
2969  watmodel = WAT_SWM4;
2970  iout << iWARN
2971  << "Setting water model to 'swm4' (SWM4-DP) for Drude polarization.\n"
2972  << endi;
2973  }
2974 
2975  // Drude water model uses "lonepairs"
2976  if (watmodel == WAT_SWM4) {
2977  lonepairs = TRUE;
2978  }
2979 
2980  // Added by JLai -- 8.2.11 -- Checks if Go method is defined
2981  if (goForcesOn) {
2982  iout << iINFO << "Go forces are on\n" << endi;
2983  // Added by JLai -- 6.3.11 -- Checks if Go method is defined
2984  int * gomethod = &goMethod;
2985  if (!opts.defined("GoMethod")) {
2986  *gomethod = 0;
2987  // printf("GO METHOD IS NOT DEFINED SO WE'LL SET IT TO SOME WEIRD VALUE\n");
2988  } else {
2989  opts.get("GoMethod",s);
2990  // printf("GO METHOD IS DEFINED SO WE'LL PRINT IT OUT: %s\n",s);
2991  *gomethod = atoi(s);
2992  }
2993  if (!strcasecmp(s, "matrix")) {
2994  goMethod = 1;
2995  //GoMethod = GO_MATRIX;
2996  } else if (!strcasecmp(s, "faster")) {
2997  goMethod = 2;
2998  //GoMethod = GO_FASTER;
2999  } else if (!strcasecmp(s, "lowmem")) {
3000  goMethod = 3;
3001  //GoMethod = GO_LOWMEM;
3002  }
3003  else {
3004  char err_msg[129];
3005  sprintf(err_msg,
3006  "Illegal value '%s' for 'GoMethod' in configuration file",
3007  s);
3008  NAMD_die(err_msg);
3009  }
3010  } // End of NAMD code to check goMethod
3011  // End of Port -- JL
3012 
3013  // Get multiple timestep integration scheme
3014  if (!opts.defined("MTSAlgorithm"))
3015  {
3016  MTSAlgorithm = VERLETI;
3017  }
3018  else
3019  {
3020  opts.get("MTSAlgorithm", s);
3021 
3022  if (!strcasecmp(s, "naive"))
3023  {
3024  MTSAlgorithm = NAIVE;
3025  }
3026  else if (!strcasecmp(s, "constant"))
3027  {
3028  MTSAlgorithm = NAIVE;
3029  }
3030  else if (!strcasecmp(s, "impulse"))
3031  {
3032  MTSAlgorithm = VERLETI;
3033  }
3034  else if (!strcasecmp(s, "verleti"))
3035  {
3036  MTSAlgorithm = VERLETI;
3037  }
3038  else
3039  {
3040  char err_msg[129];
3041 
3042  sprintf(err_msg,
3043  "Illegal value '%s' for 'MTSAlgorithm' in configuration file",
3044  s);
3045  NAMD_die(err_msg);
3046  }
3047  }
3048 
3049  // Get the long range force splitting specification
3050  if (!opts.defined("longSplitting"))
3051  {
3052  longSplitting = C1;
3053  }
3054  else
3055  {
3056  opts.get("longSplitting", s);
3057  if (!strcasecmp(s, "sharp"))
3058  longSplitting = SHARP;
3059  else if (!strcasecmp(s, "xplor"))
3060  longSplitting = XPLOR;
3061  else if (!strcasecmp(s, "c1"))
3062  longSplitting = C1;
3063  else if (!strcasecmp(s, "c2"))
3064  longSplitting = C2;
3065  else
3066  {
3067  char err_msg[129];
3068 
3069  sprintf(err_msg,
3070  "Illegal value '%s' for 'longSplitting' in configuration file",
3071  s);
3072  NAMD_die(err_msg);
3073  }
3074  }
3075 
3076  // take care of rigid bond options
3077  if (!opts.defined("rigidBonds"))
3078  {
3079  rigidBonds = RIGID_NONE;
3080  }
3081  else
3082  {
3083  opts.get("rigidBonds", s);
3084  if (!strcasecmp(s, "all"))
3085  {
3086  rigidBonds = RIGID_ALL;
3087  }
3088  else if (!strcasecmp(s, "water"))
3089  {
3090  rigidBonds = RIGID_WATER;
3091  }
3092  else if (!strcasecmp(s, "none"))
3093  {
3094  rigidBonds = RIGID_NONE;
3095  }
3096  else
3097  {
3098  char err_msg[256];
3099  sprintf(err_msg,
3100  "Illegal value '%s' for 'rigidBonds' in configuration file", s);
3101  NAMD_die(err_msg);
3102  }
3103  }
3104 
3105  // TIP4P and SWM4-DP water models require rigid water
3106  if ((watmodel == WAT_TIP4 || watmodel == WAT_SWM4)
3107  && rigidBonds == RIGID_NONE) {
3108  char err_msg[256];
3109  sprintf(err_msg,
3110  "Water model %s requires rigidBonds set to \"all\" or \"water\"",
3111  (watmodel == WAT_TIP4 ? "TIP4P" : "SWM4-DP"));
3112  NAMD_die(err_msg);
3113  }
3114 
3115  // Take care of switching stuff
3116  if (switchingActive)
3117  {
3118 
3119  if (!opts.defined("switchDist")) {
3120  NAMD_die("switchDist must be defined when switching is enabled");
3121  }
3122 
3123  if ( (switchingDist>cutoff) || (switchingDist<0) )
3124  {
3125  char err_msg[129];
3126 
3127  sprintf(err_msg,
3128  "switchDist muct be between 0 and cutoff, which is %f", cutoff);
3129  NAMD_die(err_msg);
3130  }
3131 
3132  }
3133 
3134  if ( martiniSwitching )
3135  {
3136  if ( ! switchingActive )
3137  {
3138  NAMD_die("martiniSwitching requires switching");
3139  }
3140  if ( vdwForceSwitching )
3141  {
3142  NAMD_die("martiniSwitching and vdwForceSwitching are exclusive to one another. Select only one.");
3143  }
3144  if ( dielectric != 15.0 && ! martiniDielAllow )
3145  {
3146  iout << iWARN << "USE DIELECTRIC OF 15.0 WITH MARTINI.\n";
3147  iout << iWARN << "SETTING dielectric 15.0\n";
3148  iout << iWARN << "FOR NON-STANDARD DIELECTRIC WITH MARTINI, SET: martiniDielAllow on\n";
3149  dielectric = 15.0;
3150  }
3151  if ( ! cosAngles )
3152  {
3153  iout << iWARN << "USE COSINE BASED ANGLES WITH MARTINI.\n";
3154  iout << iWARN << "SETTING cosAngles on\n";
3155  cosAngles = TRUE;
3156  }
3157  if ( PMEOn )
3158  {
3159  NAMD_die("Do not use Particle Mesh Ewald with Martini. Set: PME off");
3160  }
3161  if ( MSMOn )
3162  {
3163  NAMD_die("Do not use Multilevel Summation Method with Martini. Set: MSM off");
3164  }
3165  if ( FMMOn )
3166  {
3167  NAMD_die("Do not use Fast Multipole Method with Martini. Set: FMM off");
3168  }
3169 
3170  }
3171 
3172 
3173  if (!opts.defined("pairlistDist"))
3174  {
3175  pairlistDist = cutoff;
3176  }
3177  else if (pairlistDist < cutoff)
3178  {
3179  NAMD_die("pairlistDist must be >= cutoff distance");
3180  }
3181 
3182  patchDimension = pairlistDist;
3183 
3184  if ( splitPatch == SPLIT_PATCH_HYDROGEN ) {
3185  patchDimension += hgroupCutoff;
3186  }
3187 
3188  BigReal defaultMargin = 0.0;
3189  if (berendsenPressureOn || langevinPistonOn) {
3190  defaultMargin = ( useFlexibleCell ? 0.06 : 0.03 ) * patchDimension;
3191  }
3192  if ( margin == XXXBIGREAL ) {
3193  margin = defaultMargin;
3194  }
3195  if ( defaultMargin != 0.0 && margin == 0.0 ) {
3196  margin = defaultMargin;
3197  iout << iWARN << "ALWAYS USE NON-ZERO MARGIN WITH CONSTANT PRESSURE!\n";
3198  iout << iWARN << "CHANGING MARGIN FROM 0 to " << margin << "\n" << endi;
3199  }
3200 
3201  patchDimension += margin;
3202 
3203  //ensure patch can handle alpha_cutoff for gbis
3204  if (GBISOn) {
3205  //Check compatibility
3206  if (fullDirectOn) {
3207  NAMD_die("GBIS not compatible with FullDirect");
3208  }
3209  if (PMEOn) {
3210  NAMD_die("GBIS not compatible with PME");
3211  }
3212  if (MSMOn) {
3213  NAMD_die("GBIS not compatible with MSM");
3214  }
3215  if (FMMOn) {
3216  NAMD_die("GBIS not compatible with FMM");
3217  }
3218  if (alchOn) {
3219  NAMD_die("GBIS not compatible with Alchemical Transformations");
3220  }
3221  if (lesOn) {
3222  NAMD_die("GBIS not compatible with Locally Enhanced Sampling");
3223  }
3224  if (FMAOn) {
3225  NAMD_die("GBIS not compatible with FMA");
3226  }
3227  if (drudeOn) {
3228  NAMD_die("GBIS not compatible with Drude Polarization");
3229  }
3230 
3231  if (alpha_cutoff > patchDimension) {
3232  patchDimension = alpha_cutoff;
3233  }
3234  //calculate kappa
3235  BigReal tmp = (initialTemp > 0) ? initialTemp : 300;
3236  kappa = 50.29216*sqrt(ion_concentration/solvent_dielectric/tmp);
3237  /*magic number = 1/sqrt(eps0*kB/(2*nA*e^2*1000))*/
3238  } // GBISOn
3239 
3240  if (LCPOOn) {
3241 #ifdef MEM_OPT_VERSION
3242  NAMD_die("SASA not yet available for memory optimized builds");
3243 #endif
3244  if ( lattice.volume() > 0 ) {
3245  NAMD_die("SASA does not yet support periodic boundary conditions.");
3246  }
3247  //LCPO requires patches to be at least 16.2Ang in each dimension
3248  // twoAway[XYZ} is ignored for now
3249  }
3250 
3251  // Turn on global integration if not explicitly specified
3252 
3253  if ( dihedralOn ) globalOn = TRUE;
3254 
3255 #if defined(NAMD_CUDA) || defined(NAMD_HIP)
3256  if (loweAndersenOn) {
3257  NAMD_die("Lowe-Andersen dynamics not compatible with CUDA at this time");
3258  }
3259 #endif
3260  // END LA
3261 
3262  // BEGIN LA
3263  if (loweAndersenOn && (langevinOn || tCoupleOn))
3264  {
3265  NAMD_die("Lowe-Andersen dynamics, Langevin dynamics and temperature coupling are mutually exclusive dynamics modes");
3266  }
3267  // END LA
3268 
3269  if (tCoupleOn && opts.defined("rescaleFreq") )
3270  {
3271  NAMD_die("Temperature coupling and temperature rescaling are mutually exclusive");
3272  }
3273 
3274  if (globalOn && CkNumPes() > 1)
3275  {
3276  NAMD_die("Global integration does not run in parallel (yet).");
3277  }
3278 
3279  if (COLDOn && langevinOn)
3280  {
3281  NAMD_die("COLD and Langevin dynamics are mutually exclusive dynamics modes");
3282  }
3283  if (COLDOn && minimizeOn)
3284  {
3285  NAMD_die("COLD and minimization are mutually exclusive dynamics modes");
3286  }
3287  if (COLDOn && tCoupleOn)
3288  {
3289  NAMD_die("COLD and temperature coupling are mutually exclusive dynamics modes");
3290  }
3291  if (COLDOn && opts.defined("rescaleFreq"))
3292  {
3293  NAMD_die("COLD and velocity rescaling are mutually exclusive dynamics modes");
3294  }
3295 
3296  if (splitPatch == SPLIT_PATCH_POSITION && mollyOn )
3297  {
3298  NAMD_die("splitPatch hydrogen is required for MOLLY");
3299  }
3300 
3301  if (splitPatch == SPLIT_PATCH_POSITION && rigidBonds != RIGID_NONE)
3302  {
3303  NAMD_die("splitPatch hydrogen is required for rigidBonds");
3304  }
3305 
3306  if (accelMDOn) {
3307  if(accelMDG){
3308  char msg[128];
3309  if(accelMDGiE < 1 || accelMDGiE > 2){
3310  sprintf(msg, "accelMDGiE was set to %d but it should be 1 or 2", accelMDGiE);
3311  NAMD_die(msg);
3312  }
3313  if(accelMDGStatWindow > 0){
3314  if(accelMDGcMDPrepSteps % accelMDGStatWindow != 0)
3315  NAMD_die("'accelMDGcMDPrepSteps' has to be a multiple of 'accelMDGStatWindow'");
3316  if(accelMDGcMDSteps % accelMDGStatWindow != 0)
3317  NAMD_die("'accelMDGcMDSteps' has to be a multiple of 'accelMDGStatWindow'");
3318  if(accelMDGEquiPrepSteps % accelMDGStatWindow != 0)
3319  NAMD_die("'accelMDGEquiPrepSteps' has to be a multiple of 'accelMDGStatWindow'");
3320  if(accelMDGEquiSteps % accelMDGStatWindow != 0)
3321  NAMD_die("'accelMDGEquiSteps' has to be a multiple of 'accelMDGStatWindow'");
3322  }
3323  if(accelMDGRestart && accelMDGcMDSteps == 0)
3324  accelMDGcMDPrepSteps = 0;
3325  else if(accelMDGcMDSteps - accelMDGcMDPrepSteps < 2)
3326  NAMD_die("'accelMDGcMDSteps' should be larger than 'accelMDGcMDPrepSteps'");
3327 
3328  if(accelMDGEquiSteps == 0)
3329  accelMDGEquiPrepSteps = 0;
3330  else if(accelMDGresetVaftercmd){
3331  if(accelMDGEquiPrepSteps <= 0)
3332  NAMD_die("'accelMDGEquiPrepSteps' should be non-zero");
3333  if(accelMDGEquiSteps - accelMDGEquiPrepSteps < 1)
3334  NAMD_die("'accelMDGEquiSteps' should be larger than 'accelMDGEquiPrepSteps'");
3335  }
3336 
3337  //warn user that accelMD params will be ignored
3338  if(opts.defined("accelMDE"))
3339  iout << iWARN << "accelMDE will be ignored with accelMDG on.\n" << endi;
3340  if(opts.defined("accelMDalpha"))
3341  iout << iWARN << "accelMDalpha will be ignored with accelMDG on.\n" << endi;
3342  if(opts.defined("accelMDTE"))
3343  iout << iWARN << "accelMDTE will be ignored with accelMDG on.\n" << endi;
3344  if(opts.defined("accelMDTalpha"))
3345  iout << iWARN << "accelMDTalpha will be ignored with accelMDG on.\n" << endi;
3346  }
3347  else{
3348  if(!opts.defined("accelMDE") || !opts.defined("accelMDalpha"))
3349  NAMD_die("accelMDE and accelMDalpha are required for accelMD with accelMDG off");
3350 
3351  if(accelMDdual && (!opts.defined("accelMDTE") || !opts.defined("accelMDTalpha"))){
3352  NAMD_die("accelMDTE and accelMDTalpha are required for accelMDdual with accelMDG off");
3353  }
3354  }
3355  }
3356 
3357  // Set the default value for the maximum movement parameter
3358  // for minimization
3359  if (minimizeOn && (maximumMove == 0.0))
3360  {
3361  maximumMove = 0.75 * pairlistDist/stepsPerCycle;
3362  }
3363  if (adaptTempOn) {
3364  if (!adaptTempRescale && !adaptTempLangevin)
3365  NAMD_die("Adaptive tempering needs to be coupled to either the Langevin thermostat or velocity rescaling.");
3366  if (opts.defined("adaptTempInFile") && (opts.defined("adaptTempTmin") ||
3367  opts.defined("adaptTempTmax") ||
3368  adaptTempBins != 0))
3369  NAMD_die("cannot simultaneously specify adaptTempInFile and any of {adaptTempTmin, adaptTempTmax,adaptTempBins} as these are read from the input file");
3370  if (!opts.defined("adaptTempInFile") && !(opts.defined("adaptTempTmin") &&
3371  opts.defined("adaptTempTmax") &&
3372  adaptTempBins != 0 ))
3373  NAMD_die("Need to specify either adaptTempInFile or all of {adaptTempTmin, adaptTempTmax,adaptTempBins} if adaptTempMD is on.");
3374  }
3375 
3376  langevinOnAtStartup = langevinOn;
3377  if (langevinOn) {
3378  if ( ! opts.defined("langevinDamping") ) langevinDamping = 0.0;
3379  if ( ! opts.defined("langevinHydrogen") ) langevinHydrogen = TRUE;
3380  if ( (opts.defined("langevinDamping") || opts.defined("langevinHydrogen"))
3381  && (opts.defined("langevinFile") || opts.defined("langevinCol")) )
3382  NAMD_die("To specify Langevin dynamics parameters, use either langevinDamping and langevinHydrogen or langevinFile and langevinCol. Do not combine them.");
3383  if ( opts.defined("langevinHydrogen") && langevinDamping == 0.0 )
3384  NAMD_die("langevinHydrogen requires langevinDamping to be set.");
3385  }
3386 
3387  // BEGIN LA
3388  if (loweAndersenOn) {
3389  if (!opts.defined("loweAndersenRate")) loweAndersenRate = 100;
3390  if (!opts.defined("loweAndersenCutoff")) loweAndersenCutoff = 2.7;
3391  }
3392  // END LA
3393 
3394  // BKR - stochastic velocity rescaling
3395  if (stochRescaleOn) {
3396  if (langevinOn || loweAndersenOn || tCoupleOn ||
3397  opts.defined("rescaleFreq") || opts.defined("reassignFreq"))
3398  NAMD_die("Stochastic velocity rescaling is incompatible with other temperature control methods");
3399  // This is largely the same default used in GROMACS.
3400  if (!opts.defined("stochRescaleFreq")) stochRescaleFreq = stepsPerCycle;
3401  }
3402 
3403  if (opts.defined("rescaleFreq"))
3404  {
3405  if (!opts.defined("rescaleTemp"))
3406  {
3407  if (opts.defined("temperature"))
3408  {
3409  rescaleTemp = initialTemp;
3410  }
3411  else
3412  {
3413  NAMD_die("Must give a rescale temperature if rescaleFreq is defined");
3414  }
3415  }
3416  }
3417  else
3418  {
3419  rescaleFreq = -1;
3420  rescaleTemp = 0.0;
3421  }
3422 
3423  if (opts.defined("rescaleTemp"))
3424  {
3425  if (!opts.defined("rescaleFreq"))
3426  {
3427  NAMD_die("Must give a rescale freqency if rescaleTemp is given");
3428  }
3429  }
3430 
3431  if (opts.defined("reassignFreq"))
3432  {
3433  if (!opts.defined("reassignTemp"))
3434  {
3435  if (opts.defined("temperature"))
3436  {
3437  reassignTemp = initialTemp;
3438  }
3439  else
3440  {
3441  NAMD_die("Must give a reassign temperature if reassignFreq is defined");
3442  }
3443  }
3444  }
3445  else
3446  {
3447  reassignFreq = -1;
3448  reassignTemp = 0.0;
3449  }
3450 
3451  if (opts.defined("reassignTemp"))
3452  {
3453  if (!opts.defined("reassignFreq"))
3454  {
3455  NAMD_die("Must give a reassignment freqency if reassignTemp is given");
3456  }
3457  }
3458 
3459  if (opts.defined("reassignIncr"))
3460  {
3461  if (!opts.defined("reassignFreq"))
3462  {
3463  NAMD_die("Must give a reassignment freqency if reassignIncr is given");
3464  }
3465  }
3466  else
3467  {
3468  reassignIncr = 0.0;
3469  }
3470 
3471  if (opts.defined("reassignHold"))
3472  {
3473  if (!opts.defined("reassignIncr"))
3474  {
3475  NAMD_die("Must give a reassignment increment if reassignHold is given");
3476  }
3477  }
3478  else
3479  {
3480  reassignHold = 0.0;
3481  }
3482 
3483  if (!opts.defined("seed"))
3484  {
3485  randomSeed = (unsigned int) time(NULL) + 31530001 * CmiMyPartition();
3486  }
3487 
3488  // REST2
3489  if (opts.defined("soluteScaling")) {
3490  // Parameters soluteScalingFactorCharge and soluteScalingFactorVdw
3491  // allow independent scaling of electrostatics and van der Waals.
3492  // Initialize with soluteScalingFactor if either is not already set.
3493  if ( ! opts.defined("soluteScalingFactorCharge") ) {
3494  soluteScalingFactorCharge = soluteScalingFactor;
3495  }
3496  if ( ! opts.defined("soluteScalingFactorVdw") ) {
3497  soluteScalingFactorVdw = soluteScalingFactor;
3498  }
3499  }
3500 
3501 //fepb
3502  alchFepOnAtStartup = alchFepOn = FALSE;
3503  alchThermIntOnAtStartup = alchThermIntOn = FALSE;
3504  alchOnAtStartup = alchOn;
3505 
3506  if (alchOn) {
3507  if (martiniSwitching) {
3508  iout << iWARN << "Martini switching disabled for alchemical "
3509  "interactions.\n" << endi;
3510  }
3511 
3512  if (!opts.defined("alchType")) {
3513  NAMD_die("Must define type of alchemical simulation: fep or ti\n");
3514  }
3515  else {
3516  opts.get("alchType",s);
3517  if (!strcasecmp(s, "fep")) {
3518  alchFepOnAtStartup = alchFepOn = TRUE;
3519  }
3520  else if (!strcasecmp(s, "ti")) {
3521  alchThermIntOnAtStartup = alchThermIntOn = TRUE;
3522  }
3523  else {
3524  NAMD_die("Unknown type of alchemical simulation; choices are fep or ti\n");
3525  }
3526  }
3527 
3528  if (rescaleFreq > 0) alchTemp = rescaleTemp;
3529  else if (reassignFreq > 0) alchTemp = reassignTemp;
3530  else if (langevinOn) alchTemp = langevinTemp;
3531  else if (stochRescaleOn) alchTemp = stochRescaleTemp;
3532  else if (tCoupleOn) alchTemp = tCoupleTemp;
3533  else NAMD_die("Alchemical FEP can be performed only in constant temperature simulations\n");
3534 
3535  if (reassignFreq > 0 && reassignIncr != 0)
3536  NAMD_die("reassignIncr cannot be used in alchemical simulations\n");
3537 
3538  if (alchLambda < 0.0 || alchLambda > 1.0)
3539  NAMD_die("alchLambda values should be in the range [0.0, 1.0]\n");
3540 
3541  if (alchVdwLambdaEnd > 1.0)
3542  NAMD_die("Gosh tiny Elvis, you kicked soft-core in the van der Waals! alchVdwLambdaEnd should be in the range [0.0, 1.0]\n");
3543 
3544  if (alchBondLambdaEnd > 1.0)
3545  NAMD_die("alchBondLambdaEnd should be in the range [0.0, 1.0]\n");
3546 
3547  if (alchElecLambdaStart > 1.0)
3548  NAMD_die("alchElecLambdaStart should be in the range [0.0, 1.0]\n");
3549 
3550  if (alchWCAOn) {
3551  if (alchRepLambdaEnd > 1.0)
3552  NAMD_die("alchRepLambdaEnd should be in the range [0.0, 1.0]\n");
3553  if (alchVdwLambdaEnd < alchRepLambdaEnd)
3554  NAMD_die("alchVdwLambdaEnd should be greater than alchRepLambdaEnd\n");
3555  if (alchVdwShiftCoeff > 0.0) {
3556  iout << iWARN << "alchVdwShiftCoeff is non-zero but not used when WCA"
3557  << " is active. Setting it to zero now.\n" << endi;
3558  alchVdwShiftCoeff = 0.0;
3559  }
3560  if (alchThermIntOn) {
3561  NAMD_die("alchWCA is not currently compatible with TI");
3562  }
3563 #if defined(NAMD_CUDA) || defined(NAMD_HIP)
3564  NAMD_die("alchWCA is not currently available with CUDA");
3565 #endif
3566  }
3567 
3568  if (alchFepOn) {
3569  if (alchLambda2 < 0.0 || alchLambda2 > 1.0)
3570  NAMD_die("alchLambda2 values should be in the range [0.0, 1.0]\n");
3571 
3572  setupIDWS(); // setup IDWS if it was activated.
3573 
3574  if (!opts.defined("alchoutfile")) {
3575  strcpy(alchOutFile, outputFilename);
3576  strcat(alchOutFile, ".fep");
3577  }
3578 
3579  if (!opts.defined("alchLambda") || !opts.defined("alchLambda2")) {
3580  NAMD_die("alchFepOn is on, but alchLambda or alchLambda2 is not set.");
3581  }
3582  }
3583  else if (alchThermIntOn) {
3584  // alchLambda2 is only needed for nonequilibrium switching
3585  if (alchLambdaFreq && (alchLambda2 < 0.0 || alchLambda2 > 1.0))
3586  NAMD_die("alchLambda2 values should be in the range [0.0, 1.0]\n");
3587 
3588  if (!opts.defined("alchoutfile")) {
3589  strcpy(alchOutFile, outputFilename);
3590  strcat(alchOutFile, ".ti");
3591  }
3592  }
3593  }
3594 
3595 //fepe
3596 
3597  if ( alchOn && alchFepOn && alchThermIntOn )
3598  NAMD_die("Sorry, combined TI and FEP is not implemented.\n");
3599  if ( alchOn && lesOn )
3600  NAMD_die("Sorry, combined LES with FEP or TI is not implemented.\n");
3601  if ( alchOn && alchThermIntOn && lesOn )
3602  NAMD_die("Sorry, combined LES and TI is not implemented.\n");
3603  if ( alchWCAOn && !alchOn ) {
3604  iout << iWARN << "Alchemical WCA decomposition was requested but \
3605  alchemical free energy calculation is not active. Setting \
3606  alchWCA to off.\n" << endi;
3607  alchWCAOn = FALSE;
3608  }
3609  if ( alchDecouple && !alchOn ) {
3610  iout << iWARN << "Alchemical decoupling was requested but \
3611  alchemical free energy calculation is not active. Setting \
3612  alchDecouple to off.\n" << endi;
3613  alchDecouple = FALSE;
3614  }
3615  if ( alchBondDecouple && !alchOn ) {
3616  iout << iWARN << "Alchemical bond decoupling was requested but \
3617  alchemical free energy calculation is not active. Setting \
3618  alchBondDecouple to off.\n" << endi;
3619  alchBondDecouple = FALSE;
3620  }
3621 
3622  if ( lesOn && ( lesFactor < 1 || lesFactor > 255 ) ) {
3623  NAMD_die("lesFactor must be positive and less than 256");
3624  }
3625  if ((pairInteractionOn && alchOn) || (pairInteractionOn && lesOn))
3626  NAMD_die("Sorry, pair interactions may not be calculated when LES, FEP or TI is enabled.");
3627 
3628  // Drude model
3629  if (drudeOn) {
3630  if ( ! langevinOn ) {
3631  NAMD_die("Drude model requires use of Langevin thermostat.");
3632  }
3633  if ( ! opts.defined("drudeDamping")) {
3634  drudeDamping = langevinDamping;
3635  iout << iWARN << "Undefined 'drudeDamping' will be set to "
3636  "value of 'langevinDamping'\n" << endi;
3637  }
3638  if ( alchOn ) {
3639  NAMD_die("Drude implementation is incompatible with alchemical "
3640  "free energy calculation.");
3641  }
3642  }
3643 
3644  // Set up load balancing variables
3645  if (opts.defined("ldBalancer")) {
3646  if (strcasecmp(loadBalancer, "none") == 0)
3647  ldBalancer = LDBAL_NONE;
3648  else if (strcasecmp(loadBalancer, "hybrid") == 0)
3649  ldBalancer = LDBAL_HYBRID;
3650  else
3651  NAMD_die("Unknown ldBalancer selected");
3652  } else {
3653  ldBalancer = LDBAL_CENTRALIZED;
3654 #ifdef MEM_OPT_VERSION
3655  if ( CkNumPes() > 1400 ) ldBalancer = LDBAL_HYBRID;
3656 #endif
3657  }
3658 
3659  if (opts.defined("ldbStrategy")) {
3660  // Assign the load balancing strategy
3661  if (strcasecmp(loadStrategy, "comprehensive") == 0)
3662  ldbStrategy = LDBSTRAT_COMPREHENSIVE;
3663  else if (strcasecmp(loadStrategy, "refineonly") == 0)
3664  ldbStrategy = LDBSTRAT_REFINEONLY;
3665  else if (strcasecmp(loadStrategy, "old") == 0)
3666  ldbStrategy = LDBSTRAT_OLD;
3667  else
3668  NAMD_die("Unknown ldbStrategy selected");
3669  } else {
3670  ldbStrategy = LDBSTRAT_DEFAULT;
3671  }
3672 
3673  if (!opts.defined("ldbPeriod")) {
3674  ldbPeriod=200*stepsPerCycle;
3675  }
3676 
3677  // Set default values
3678  if (!opts.defined("firstLdbStep")) {
3679  firstLdbStep=5*stepsPerCycle;
3680  }
3681 
3682  if (ldbPeriod <= firstLdbStep) {
3683  NAMD_die("ldbPeriod must greater than firstLdbStep.");
3684  }
3685 
3686  if (!opts.defined("lastLdbStep")) {
3687  lastLdbStep = -1;
3688  }
3689 
3690  if (!opts.defined("hybridGroupSize")) {
3691  hybridGroupSize = 512;
3692  }
3693  if ( hybridGroupSize < CkNumPes() ) {
3694  // match load balancer boundaries to physical nodes if possible
3695  int groupsize = hybridGroupSize;
3696  int *rpelist;
3697  int nodesize;
3698  CmiGetPesOnPhysicalNode(CmiPhysicalNodeID(0), &rpelist, &nodesize);
3699  if ( CkNumPes() % nodesize ) nodesize = CmiNodeSize(CmiNodeOf(0));
3700  if ( CkNumPes() % nodesize ) nodesize = 1;
3701  groupsize += nodesize - 1;
3702  while ( 2 * groupsize > CkNumPes() ) --groupsize;
3703  if ( groupsize < nodesize ) groupsize = nodesize;
3704  while ( groupsize % nodesize ) --groupsize;
3705  while ( groupsize && CkNumPes() % groupsize ) groupsize -= nodesize;
3706  if ( 2 * groupsize < hybridGroupSize ) {
3707  groupsize += nodesize;
3708  while ( CkNumPes() % groupsize ) groupsize += nodesize;
3709  }
3710  if ( 2 * groupsize <= CkNumPes() ) hybridGroupSize = groupsize;
3711  }
3712 
3713  // tracing will be done if trace is available and user says +traceOff
3714  // in that case we set nice values for some functions
3715  bool specialTracing = traceAvailable() && (traceIsOn() == 0);
3716 
3717  if(!opts.defined("traceStartStep")) {
3718  traceStartStep = 4 * firstLdbStep + 2 * ldbPeriod;
3719  }
3720  if(!opts.defined("numTraceSteps")) {
3721  numTraceSteps = 100;
3722  }
3723 
3724  if(specialTracing) {
3725  if (!opts.defined("firstLdbStep")) firstLdbStep = 20;
3726  if (!opts.defined("ldbPeriod")) ldbPeriod = 100;
3727 
3728  if(!opts.defined("traceStartStep")) {
3729  traceStartStep = 4 * firstLdbStep + 2 * ldbPeriod; // 380
3730  }
3731 
3732  if(!opts.defined("numTraceSteps")) {
3733  numTraceSteps = 80;
3734  }
3735  }
3736 
3737 #ifdef MEASURE_NAMD_WITH_PAPI
3738  if(papiMeasure){
3739  if(!opts.defined("papiMeasureStartStep")) {
3740  papiMeasureStartStep = 3 * firstLdbStep;
3741  }
3742  if(!opts.defined("numPapiMeasureSteps")) {
3743  numPapiMeasureSteps = 8; //including two pme steps
3744  }
3745  }
3746 #endif
3747 
3748  if(simulateInitialMapping) {
3749  if(!opts.defined("simulatedPEs")){
3750  simulatedPEs = CkNumPes();
3751  }
3752  if(!opts.defined("simulatedNodeSize")){
3753  simulatedNodeSize = CkMyNodeSize();
3754  }
3755  }
3756 
3757 #ifdef MEM_OPT_VERSION
3758  //Some constraints on the values of load balancing parameters.
3759  //The reason is related to communication schemes used in sending proxy
3760  //data. If the step immediately after the load balancing is not a step
3761  //for atom migration, then it's possible there are some necessary information
3762  // missing inside the ProxyPatch which will crash the program. Therefore,
3763  // It's better that the step immediately after the load balancing be a step
3764  // for atom migration so that the some overhead in Proxy msgs are removed.
3765  // --Chao Mei
3766  if(ldbPeriod%stepsPerCycle!=0 || firstLdbStep%stepsPerCycle!=0) {
3767  iout << iWARN << "In memory optimized version, the ldbPeriod parameter or firstLdbStep parameter is better set to be a multiple of stepsPerCycle parameter!\n";
3768  }
3769 #endif
3770 
3771  if (N < firstTimestep) { N = firstTimestep; }
3772 
3773  if ( (firstTimestep%stepsPerCycle) != 0)
3774  {
3775  NAMD_die("First timestep must be a multiple of stepsPerCycle!!");
3776  }
3777 
3778  // Make sure only one full electrostatics algorithm is selected
3779  {
3780  int i = 0;
3781  if ( FMAOn ) ++i;
3782  if ( PMEOn ) ++i;
3783  if ( MSMOn ) ++i;
3784  if ( FMMOn ) ++i;
3785  if ( fullDirectOn ) ++i;
3786  if ( i > 1 )
3787  NAMD_die("More than one full electrostatics algorithm selected!!!");
3788  }
3789 
3790  if (!opts.defined("ldbBackgroundScaling")) {
3791  ldbBackgroundScaling = 1.0;
3792  }
3793  if (!opts.defined("ldbPMEBackgroundScaling")) {
3794  ldbPMEBackgroundScaling = ldbBackgroundScaling;
3795  }
3796  if (!opts.defined("ldbHomeBackgroundScaling")) {
3797  ldbHomeBackgroundScaling = ldbBackgroundScaling;
3798  }
3799 
3800  // Check on PME parameters
3801  if (PMEOn) { // idiot checking
3802  if ( lattice.volume() == 0. ) {
3803  NAMD_die("PME requires periodic boundary conditions.");
3804  }
3805  if ( PMEGridSpacing == 0. ) {
3806  if ( PMEGridSizeX * PMEGridSizeY * PMEGridSizeZ == 0 )
3807  NAMD_die("Either PMEGridSpacing or PMEGridSizeX, PMEGridSizeY, and PMEGridSizeZ must be specified.");
3808  else PMEGridSpacing = 1.5; // only exit in very bad cases
3809  }
3810 #ifndef TEST_PME_GRID
3811  for ( int idim = 0; idim < 3; ++idim ) {
3812  int *gridSize;
3813  BigReal cellLength;
3814  const char *direction;
3815  switch ( idim ) {
3816  case 0: direction = "X";
3817  gridSize = &PMEGridSizeX; cellLength = lattice.a().length();
3818  break;
3819  case 1: direction = "Y";
3820  gridSize = &PMEGridSizeY; cellLength = lattice.b().length();
3821  break;
3822  case 2: direction = "Z";
3823  gridSize = &PMEGridSizeZ; cellLength = lattice.c().length();
3824  break;
3825  }
3826  int minSize = (int) ceil(cellLength/PMEGridSpacing);
3827 #else
3828  for ( int minSize = 1; minSize < 300; ++minSize ) {
3829 #endif
3830  int bestSize = 10 * (minSize + 10); // make sure it's big
3831  int max2, max3, ts;
3832  for ( max2=2, ts=1; ts < minSize; ++max2 ) ts *= 2;
3833  for ( max3=2, ts=1; ts < minSize; ++max3 ) ts *= 3;
3834  int max5 = 2;
3835  int max7 = 1;
3836  int max11 = 1;
3837  for ( int i2 = 0; i2 <= max2; ++i2 ) {
3838  for ( int i3 = 0; i3 <= max3; ++i3 ) {
3839  for ( int i5 = 0; i5 <= max5; ++i5 ) {
3840  for ( int i7 = 0; i7 <= max7; ++i7 ) {
3841  for ( int i11 = 0; i11 <= max11; ++i11 ) {
3842  if ( i5 + i7 + i11 > i2 ) continue;
3843  int testSize = 2; // must be even
3844  for ( int j2 = 0; j2 < i2; ++j2 ) testSize *= 2;
3845  if ( testSize > bestSize ) continue;
3846  for ( int j3 = 0; j3 < i3; ++j3 ) testSize *= 3;
3847  if ( testSize > bestSize ) continue;
3848  for ( int j5 = 0; j5 < i5; ++j5 ) testSize *= 5;
3849  if ( testSize > bestSize ) continue;
3850  for ( int j7 = 0; j7 < i7; ++j7 ) testSize *= 7;
3851  if ( testSize > bestSize ) continue;
3852  for ( int j11 = 0; j11 < i11; ++j11 ) testSize *= 11;
3853  if ( testSize > bestSize ) continue;
3854  if ( testSize >= minSize ) bestSize = testSize;
3855  } } } } }
3856 #ifdef TEST_PME_GRID
3857  iout << minSize << " " << bestSize << "\n" << endi;
3858 #else
3859  if ( ! *gridSize ) { // set it
3860  *gridSize = bestSize;
3861  }
3862  if ( *gridSize * PMEGridSpacing < cellLength ) {
3863  char errmsg[512];
3864  sprintf(errmsg, "PMEGridSize%s %d is too small for cell length %f and PMEGridSpacing %f\n",
3865  direction, *gridSize, cellLength, PMEGridSpacing);
3866  NAMD_die(errmsg);
3867  }
3868 #endif
3869  }
3870  if ( PMEGridSizeX < 5 ) {
3871  NAMD_die("PMEGridSizeX (number of grid points) is very small.");
3872  }
3873  if ( PMEGridSizeY < 5 ) {
3874  NAMD_die("PMEGridSizeY (number of grid points) is very small.");
3875  }
3876  if ( PMEGridSizeZ < 5 ) {
3877  NAMD_die("PMEGridSizeZ (number of grid points) is very small.");
3878  }
3879  BigReal tolerance = PMETolerance;
3880  BigReal ewaldcof = 1.0;
3881  while ( erfc(ewaldcof*cutoff)/cutoff >= tolerance ) ewaldcof *= 2.0;
3882  BigReal ewaldcof_lo = 0.;
3883  BigReal ewaldcof_hi = ewaldcof;
3884  for ( int i = 0; i < 100; ++i ) {
3885  ewaldcof = 0.5 * ( ewaldcof_lo + ewaldcof_hi );
3886  if ( erfc(ewaldcof*cutoff)/cutoff >= tolerance ) {
3887  ewaldcof_lo = ewaldcof;
3888  } else {
3889  ewaldcof_hi = ewaldcof;
3890  }
3891  }
3892  PMEEwaldCoefficient = ewaldcof;
3893 
3894 #if defined(NAMD_CUDA) || defined(NAMD_HIP)
3895  bool one_device_per_node = deviceCUDA->one_device_per_node(); // only checks node 0
3896  if ( ! opts.defined("PMEOffload") ) {
3897  PMEOffload = ( (PMEInterpOrder > 4) && one_device_per_node );
3898  if ( PMEOffload ) iout << iINFO << "Enabling PMEOffload because PMEInterpOrder > 4.\n" << endi;
3899  } else if ( PMEOffload && ! one_device_per_node ) {
3900  PMEOffload = 0;
3901  iout << iWARN << "Disabling PMEOffload because multiple CUDA devices per process are not supported.\n" << endi;
3902  }
3903  if ( usePMECUDA && ! ( useCUDA2 || one_device_per_node ) ) {
3904  usePMECUDA = 0;
3905  iout << iWARN << "Disabling usePMECUDA because multiple CUDA devices per process requires useCUDA2.\n" << endi;
3906  }
3907 #else
3908  PMEOffload = 0;
3909 #endif
3910  } else { // initialize anyway
3911  useDPME = 0;
3912  PMEGridSizeX = 0;
3913  PMEGridSizeY = 0;
3914  PMEGridSizeZ = 0;
3915  PMEGridSpacing = 1000.;
3916  PMEEwaldCoefficient = 0;
3917  PMEOffload = 0;
3918  }
3919 
3920  // Take care of initializing FMA values to something if FMA is not
3921  // active
3922  if (!FMAOn)
3923  {
3924  FMALevels = 0;
3925  FMAMp = 0;
3926  FMAFFTOn = FALSE;
3927  FMAFFTBlock = 0;
3928  }
3929  else
3930  {
3931  // idiot checking: frm bug reported by Tom Bishop.
3932  // DPMTA requires: (#terms in fma)/(fft blocking factor) = integer.
3933  if (FMAFFTBlock != 4)
3934  NAMD_die("FMAFFTBlock: Block length must be 4 for short FFT's");
3935  if (FMAMp % FMAFFTBlock != 0)
3936  NAMD_die("FMAMp: multipole term must be multiple of block length (FMAFFTBlock)");
3937  }
3938 
3939  if ( (nonbondedFrequency > stepsPerCycle) || ( (stepsPerCycle % nonbondedFrequency) != 0) )
3940  {
3941  NAMD_die("stepsPerCycle must be a multiple of nonbondedFreq");
3942  }
3943 
3944  if (!LCPOOn && !GBISOn && !GBISserOn && !FMAOn && !PMEOn && !MSMOn && !fullDirectOn && !FMMOn)
3945  {
3946  fullElectFrequency = 0;
3947  }
3948  else
3949  {
3950  if (!opts.defined("fullElectFrequency"))
3951  {
3952  if (opts.defined("fmaFrequency")) {
3953  iout << iWARN << "The parameter fmaFrequency has been renamed fullElectFrequency.\n" << endi;
3954  fullElectFrequency = fmaFrequency;
3955  } else {
3956  iout << iWARN << "The parameter fullElectFrequency now defaults to nonbondedFreq (" << nonbondedFrequency << ") rather than stepsPerCycle.\n" << endi;
3957  fullElectFrequency = nonbondedFrequency;
3958  }
3959  }
3960  else
3961  {
3962  if (opts.defined("fmaFrequency")) {
3963  iout << iWARN << "Ignoring redundant parameter fmaFrequency in favor of fullElectFrequency.\n" << endi;
3964  }
3965  if ( (fullElectFrequency > stepsPerCycle) || ( (stepsPerCycle % fullElectFrequency) != 0) )
3966  {
3967  NAMD_die("stepsPerCycle must be a multiple of fullElectFrequency");
3968  }
3969  }
3970 
3971  if ( (nonbondedFrequency > fullElectFrequency) || ( (fullElectFrequency % nonbondedFrequency) != 0) )
3972  {
3973  NAMD_die("fullElectFrequency must be a multiple of nonbondedFreq");
3974  }
3975 
3976  if (singleTopology && fullElectFrequency > 1) NAMD_die("Single topology free energy calculation discourages multiple timesteps to assure accuracy!");
3977  if (singleTopology && alchDecouple) NAMD_die("Single topology free energy calculation can NOT work with alchDecouple on");
3978 
3979  if (multigratorOn) {
3980  if ( (multigratorTemperatureFreq > multigratorPressureFreq) || ( (multigratorPressureFreq % multigratorTemperatureFreq) != 0) )
3981  {
3982  NAMD_die("multigratorTemperatureFreq must be a multiple of multigratorPressureFreq");
3983  }
3984  if ( (fullElectFrequency > multigratorTemperatureFreq) || ( (multigratorTemperatureFreq % fullElectFrequency) != 0) )
3985  {
3986  NAMD_die("fullElectFrequency must be a multiple of multigratorTemperatureFreq");
3987  }
3988  if (multigratorNoseHooverChainLength <= 2) {
3989  NAMD_die("multigratorNoseHooverChainLength must be greater than 2");
3990  }
3991  }
3992 
3993  if (!opts.defined("fmaTheta"))
3994  fmaTheta=0.715; /* Suggested by Duke developers */
3995  }
3996 
3997  if ( lesOn && ( FMAOn || useDPME || fullDirectOn ) ) {
3998  NAMD_die("Sorry, LES is only implemented for PME full electrostatics.");
3999  }
4000  if ( alchFepOn && ( FMAOn || useDPME || fullDirectOn ) ) {
4001  NAMD_die("Sorry, FEP is only implemented for PME full electrostatics.");
4002  }
4003  if ( alchThermIntOn && ( FMAOn || useDPME || fullDirectOn ) ) {
4004  NAMD_die("Sorry, TI is only implemented for PME full electrostatics.");
4005  }
4006  if ( pairInteractionOn && FMAOn ) {
4007  NAMD_die("Sorry, pairInteraction not implemented for FMA.");
4008  }
4009  if ( pairInteractionOn && useDPME ) {
4010  NAMD_die("Sorry, pairInteraction not implemented for DPME.");
4011  }
4012  if ( pairInteractionOn && fullDirectOn ) {
4013  NAMD_die("Sorry, pairInteraction not implemented for full direct electrostatics.");
4014  }
4015  if ( ! pairInteractionOn ) {
4016  pairInteractionSelf = 0;
4017  }
4018  if ( pairInteractionOn && !pairInteractionSelf && !config->find("pairInteractionGroup2"))
4019  NAMD_die("pairInteractionGroup2 must be specified");
4020 
4021  if ( ! fixedAtomsOn ) {
4022  fixedAtomsForces = 0;
4023  }
4024 
4025  if ( gridforceOn || mgridforceOn ) {
4026  parse_mgrid_params(config);
4027  }
4028 
4029  if ( extraBondsOn ) {
4030  extraBondsCosAnglesSetByUser = ! ! config->find("extraBondsCosAngles");
4031  } else {
4032  extraBondsCosAnglesSetByUser = false;
4033  }
4034 
4035  if (!opts.defined("constraints"))
4036  {
4037  constraintExp = 0;
4038  constraintScaling = 1.0;
4039 
4040  //****** BEGIN selective restraints (X,Y,Z) changes
4041  selectConstraintsOn = FALSE;
4042  //****** END selective restraints (X,Y,Z) changes
4043 
4044  //****** BEGIN moving constraints changes
4045  movingConstraintsOn = FALSE;
4046  //****** END moving constraints changes
4047  //****** BEGIN rotating constraints changes
4048  rotConstraintsOn = FALSE;
4049  //****** END rotating constraints changes
4050  }
4051  //****** BEGIN rotating constraints changes
4052  else {
4053  if (rotConstraintsOn) {
4054  rotConsAxis = rotConsAxis.unit();
4055  }
4056  }
4057  if(opts.defined("rotConstraints")
4058  && opts.defined("movingConstraints")) {
4059  NAMD_die("Rotating and moving constraints are mutually exclusive!");
4060  }
4061  //****** END rotating constraints changes
4062 
4063  //****** BEGIN selective restraints (X,Y,Z) changes
4064  if(opts.defined("selectConstraints") && !opts.defined("selectConstrX")
4065  && !opts.defined("selectConstrY") && !opts.defined("selectConstrZ")) {
4066  NAMD_die("selectConstraints was specified, but no Cartesian components were defined!");
4067  }
4068  if (!opts.defined("selectConstraints")) {
4069  constrXOn = FALSE;
4070  constrYOn = FALSE;
4071  constrZOn = FALSE;
4072  }
4073  //****** END selective restraints (X,Y,Z) changes
4074 
4075 
4076  //****** BEGIN SMD constraints changes
4077 
4078  if (!opts.defined("SMD")) {
4079  SMDOn = FALSE;
4080  }
4081 
4082  if (SMDOn) {
4083  // normalize direction
4084  if (SMDDir.length2() == 0) {
4085  NAMD_die("SMD direction vector must be non-zero");
4086  }
4087  else {
4088  SMDDir = SMDDir.unit();
4089  }
4090 
4091  if (SMDOutputFreq > 0 && SMDOutputFreq < stepsPerCycle
4092  || SMDOutputFreq % stepsPerCycle != 0) {
4093  NAMD_die("SMDOutputFreq must be a multiple of stepsPerCycle");
4094  }
4095  }
4096 
4097  //****** END SMD constraints changes
4098 
4099  if (!sphericalBCOn)
4100  {
4101  sphericalBCr1 = 0.0;
4102  sphericalBCk1 = 0.0;
4103  sphericalBCexp1 = 0;
4104  sphericalBCr2 = 0.0;
4105  sphericalBCk2 = 0.0;
4106  sphericalBCexp2 = 0;
4107  }
4108  else if (!opts.defined("sphericalBCr2"))
4109  {
4110  sphericalBCr2 = -1.0;
4111  sphericalBCk2 = 0.0;
4112  sphericalBCexp2 = 0;
4113  }
4114 
4115  if (!cylindricalBCOn)
4116  {
4117  cylindricalBCr1 = 0.0;
4118  cylindricalBCk1 = 0.0;
4119  cylindricalBCexp1 = 0;
4120  cylindricalBCr2 = 0.0;
4121  cylindricalBCk2 = 0.0;
4122  cylindricalBCexp2 = 0;
4123  cylindricalBCl1 = 0.0;
4124  cylindricalBCl2 = 0.0;
4125  }
4126  else if (!opts.defined("cylindricalBCr2"))
4127  {
4128  cylindricalBCr2 = -1.0;
4129  cylindricalBCk2 = 0.0;
4130  cylindricalBCexp2 = 0;
4131  cylindricalBCl2 = 0.0;
4132  }
4133 
4134  if (!eFieldOn)
4135  {
4136  eField.x = 0.0;
4137  eField.y = 0.0;
4138  eField.z = 0.0;
4139  eFieldFreq = 0.0;
4140  eFieldPhase = 0.0;
4141  }
4142  else
4143  {
4144  if (!opts.defined("eFieldFreq")) eFieldFreq = 0.0;
4145  if (!opts.defined("eFieldPhase")) eFieldPhase = 0.0;
4146  }
4147 
4148  if (!stirOn)
4149  {
4150  stirFilename[0] = STRINGNULL;
4151  stirStartingTheta = 0.0;
4152  stirVel = 0.0;
4153  stirK = 0.0;
4154  stirAxis.x = 0.0;
4155  stirAxis.y = 0.0;
4156  stirAxis.z = 0.0;
4157  stirPivot.x = 0.0;
4158  stirPivot.y = 0.0;
4159  stirPivot.z = 0.0;
4160  }
4161 
4162  if (!opts.defined("langevin"))
4163  {
4164  langevinTemp = 0.0;
4165  }
4166 
4167  // BEGIN LA
4168  if (!opts.defined("loweAndersen"))
4169  {
4170  loweAndersenTemp = 0.0;
4171  }
4172  // END LA
4173 
4174  if (!opts.defined("tcouple"))
4175  {
4176  tCoupleTemp = 0.0;
4177  }
4178 
4179  if (HydrogenBonds)
4180  {
4181  if (daCutoffDist > pairlistDist)
4182  NAMD_die("Hydrogen bond cutoff distance must be <= pairlist distance");
4183  }
4184 
4185  // If we're doing pair interaction, set
4186  // outputEnergies to 1 to make NAMD not die (the other nonbonded code paths
4187  // aren't defined when these options are enabled), and set nonbondedFreq to
4188  // 1 to avoid getting erroneous output. Warn the user of what we're doing.
4189  if (pairInteractionOn) {
4190  if (outputEnergies != 1) {
4191  iout << iWARN << "Setting outputEnergies to 1 due to\n";
4192  iout << iWARN << "pairInteraction calculations\n" << endi;
4193  outputEnergies = 1;
4194  }
4195  }
4196  if (pairInteractionOn || pressureProfileOn) {
4197  if (nonbondedFrequency != 1) {
4198  iout << iWARN << "Setting nonbondedFreq to 1 due to\n";
4199  iout << iWARN << "pairInteraction or pressure profile calculations\n" << endi;
4200  }
4201  }
4202 
4203  // print timing at a reasonable interval by default
4204  if (!opts.defined("outputTiming"))
4205  {
4206  outputTiming = firstLdbStep;
4207  int ot2 = 10 * outputEnergies;
4208  if ( outputTiming < ot2 ) outputTiming = ot2;
4209  }
4210 
4211  // Checks if a secondary process was added in the configuration, and sets
4212  // the appropriated variable
4213  if(qmForcesOn){
4214 
4215  if (opts.defined("QMSecProc")){
4216  qmSecProcOn = true;
4217  }
4218  else {
4219  qmSecProcOn = false;
4220  }
4221 
4222  if (opts.defined("qmPrepProc")){
4223  qmPrepProcOn = true;
4224  }
4225  else {
4226  qmPrepProcOn = false;
4227  }
4228 
4229  if (opts.defined("QMParamPDB")){
4230  qmParamPDBDefined = true;
4231  }
4232  else {
4233  qmParamPDBDefined = false;
4234  }
4235 
4236  if (opts.defined("QMBondColumn")){
4237  qmBondOn = true;
4238  }
4239  else {
4240  qmBondOn = false;
4241  }
4242 
4243  if ( strcasecmp(qmSoftware,"orca") != 0 &&
4244  strcasecmp(qmSoftware,"mopac") != 0 &&
4245  strcasecmp(qmSoftware,"custom") != 0 ) {
4246  NAMD_die("Available QM software options are \'mopac\', \'orca\', or \'custom\'.");
4247  }
4248  else {
4249  if ( strcasecmp(qmSoftware,"orca") == 0 )
4250  qmFormat = QMFormatORCA;
4251  if ( strcasecmp(qmSoftware,"mopac") == 0 )
4252  qmFormat = QMFormatMOPAC;
4253  if ( strcasecmp(qmSoftware,"custom") == 0 )
4254  qmFormat = QMFormatUSR;
4255 
4256  if (qmFormat == QMFormatORCA || qmFormat == QMFormatMOPAC) {
4257 
4258  if (! opts.defined("QMConfigLine"))
4259  NAMD_die("If the selected QM software is \'mopac\' or \'orca\'\
4260 , QMConfigLine needs to be defined.");
4261 
4262  }
4263  }
4264 
4265  qmChrgMode = QMCHRGMULLIKEN;
4266  if (opts.defined("QMChargeMode")) {
4267  if ( strcasecmp(qmChrgModeS,"none") != 0 &&
4268  strcasecmp(qmChrgModeS,"mulliken") != 0 &&
4269  strcasecmp(qmChrgModeS,"chelpg") != 0) {
4270  NAMD_die("Available charge options are \'none\', \'mulliken\' or \'chelpg\'.");
4271  }
4272  else {
4273  if ( strcasecmp(qmChrgModeS,"none") == 0 )
4274  qmChrgMode = QMCHRGNONE;
4275  if ( strcasecmp(qmChrgModeS,"mulliken") == 0 )
4276  qmChrgMode = QMCHRGMULLIKEN;
4277  if ( strcasecmp(qmChrgModeS,"chelpg") == 0 )
4278  qmChrgMode = QMCHRGCHELPG;
4279  }
4280  }
4281 
4282  if (qmFormat == QMFormatMOPAC && qmChrgMode == QMCHRGCHELPG)
4283  NAMD_die("Available charge options for MOPAC are \'none\' and \'mulliken\'.");
4284 
4285  if (qmFormat == QMFormatUSR && qmChrgMode == QMCHRGCHELPG)
4286  NAMD_die("Available charge options for MOPAC are \'none\' and \'mulliken\'.");
4287 
4288  if (qmBondOn && (opts.defined("QMBondValueType"))) {
4289  if ( strcasecmp(qmBondValueTypeS,"len") != 0 &&
4290  strcasecmp(qmBondValueTypeS,"ratio") != 0 ) {
4291  NAMD_die("Available QM bond value type options are \'len\' or \'ratio\'.");
4292  }
4293  else {
4294  // #define QMLENTYPE 1
4295  // #define QMRATIOTYPE 2
4296  if ( strcasecmp(qmBondValueTypeS,"len") == 0 )
4297  qmBondValType = 1;
4298  if ( strcasecmp(qmBondValueTypeS,"ratio") == 0 )
4299  qmBondValType = 2;
4300  }
4301  }
4302  else if (qmBondOn && ! (opts.defined("QMBondValueType")))
4303  qmBondValType = 1;
4304 
4305  if ( strcmp(qmColumn,"beta") != 0 &&
4306  strcmp(qmColumn,"occ") != 0 ) {
4307  NAMD_die("Available column options are \'beta\' and \'occ\'.");
4308  }
4309 
4310  if (qmBondOn) {
4311  if ( strcmp(qmBondColumn,"beta") != 0 &&
4312  strcmp(qmBondColumn,"occ") != 0 ) {
4313  NAMD_die("Available column options are \'beta\' and \'occ\'.");
4314  }
4315 
4316  if (strcmp(qmBondColumn,qmColumn) == 0)
4317  NAMD_die("QM column and bond-column must be different!");
4318  }
4319 
4320  qmBondScheme = 1;
4321  if (opts.defined("QMBondScheme")) {
4322  if ( strcasecmp(qmBondSchemeS,"CS") == 0 )
4323  qmBondScheme = QMSCHEMECS;
4324  if ( strcasecmp(qmBondSchemeS,"RCD") == 0 )
4325  qmBondScheme = QMSCHEMERCD;
4326  if ( strcasecmp(qmBondSchemeS,"Z1") == 0 )
4327  qmBondScheme = QMSCHEMEZ1;
4328  if ( strcasecmp(qmBondSchemeS,"Z2") == 0 )
4329  qmBondScheme = QMSCHEMEZ2;
4330  if ( strcasecmp(qmBondSchemeS,"Z3") == 0 )
4331  qmBondScheme = QMSCHEMEZ3;
4332  }
4333 
4334 // #define QMPCSCHEMENONE 1
4335 // #define QMPCSCHEMEROUND 2
4336 // #define QMPCSCHEMEZERO 3
4337  qmPCScheme = 1;
4338  if (opts.defined("QMPointChargeScheme") && qmPCSwitchOn) {
4339  if ( strcasecmp(qmPCSchemeS,"none") == 0 )
4340  qmPCScheme = 1;
4341 
4342  if ( strcasecmp(qmPCSchemeS,"round") == 0 )
4343  qmPCScheme = 2;
4344  if ( strcasecmp(qmPCSchemeS,"zero") == 0 )
4345  qmPCScheme = 3;
4346 
4347  if ( qmPCScheme > 1 && ! qmPCSwitchOn)
4348  NAMD_die("QM Charge Schemes \'round\' or \'zero\' can only be applied with QMswitching set to \'on\'!");
4349  }
4350 
4351  // Redundant option to deprecate "qmNoPC" option.
4352  if (qmElecEmbed)
4353  qmNoPC = FALSE;
4354 
4355 // #define QMLSSMODEDIST 1
4356 // #define QMLSSMODECOM 2
4357  if (qmLSSOn) {
4358 
4359  if (qmNoPC)
4360  NAMD_die("QM Live Solvent Selection cannot be done with QMNoPntChrg set to \'on\'!") ;
4361 
4362  if (rigidBonds != RIGID_NONE)
4363  NAMD_die("QM Live Solvent Selection cannot be done with fixed bonds!") ;
4364 
4365  if (qmLSSFreq % qmPCSelFreq != 0)
4366  NAMD_die("Frequency of QM solvent update must be a multiple of frequency of point charge selection.");
4367 
4368  if (qmLSSFreq % stepsPerCycle != 0)
4369  NAMD_die("Frequency of QM solvent update must be a multiple of steps per cycle.");
4370 
4371  if (opts.defined("QMLSSMode") ) {
4372  if ( strcasecmp(qmLSSModeS,"dist") != 0 &&
4373  strcasecmp(qmLSSModeS,"COM") != 0 ) {
4374  NAMD_die("Available LSS mode options are \'dist\' and \'COM\'.");
4375  }
4376  if ( strcasecmp(qmLSSModeS,"dist") == 0 )
4377  qmLSSMode = 1;
4378  else if ( strcasecmp(qmLSSModeS,"COM") == 0 )
4379  qmLSSMode = 2;
4380  }
4381  else
4382  qmLSSMode = 1;
4383  }
4384 
4385 // #define QMPCSCALESHIFT 1
4386 // #define QMPCSCALESWITCH 2
4387  if (qmPCSwitchOn) {
4388 
4389  if (opts.defined("QMSwitchingType") ) {
4390  if ( strcasecmp(qmPCSwitchTypeS,"shift") != 0 &&
4391  strcasecmp(qmPCSwitchTypeS,"switch") != 0 ) {
4392  NAMD_die("Available scaling options are \'shift\' and \'switch\'.");
4393  }
4394  if ( strcasecmp(qmPCSwitchTypeS,"shift") == 0 )
4395  qmPCSwitchType = 1;
4396  else if ( strcasecmp(qmPCSwitchTypeS,"switch") == 0 )
4397  qmPCSwitchType = 2;
4398  }
4399  else
4400  qmPCSwitchType = 1;
4401  }
4402 
4403  if (qmNoPC && qmPCSelFreq > 1) {
4404  iout << iWARN << "QMPCStride being IGNORED since QMNoPntChrg is set to \'on\'!\n" << endi;
4405  qmPCSelFreq = 1;
4406  }
4407 
4408  if (qmNoPC && qmPCSwitchOn)
4409  NAMD_die("QM PC switching can only be applied with QMNoPntChrg set to \'off\'!");
4410 
4411 // if (qmNoPC && qmBondOn)
4412 // NAMD_die("QM-MM bonds can only be applied with QMNoPntChrg set to \'off\'!");
4413 
4414  if (qmPCSelFreq <= 0)
4415  NAMD_die("QMPCFreq can only be a positive number! For static point charge selection, see QMCutomPC.");
4416 
4417  if (qmCustomPCSel && qmNoPC)
4418  NAMD_die("QM Custom PC Selection is incompatible with QMNoPntChrg!");
4419 
4420 // if (qmCustomPCSel && qmPCSwitchOn)
4421 // NAMD_die("QM Custom PC Selection is incompatible with QMSwitching!");
4422 
4423  if (qmCustomPCSel && qmPCSelFreq > 1)
4424  NAMD_die("QM Custom PC Selection is incompatible with QMPCStride > 1!");
4425 
4426  if (qmCSMD && (! opts.defined("QMCSMDFile") ))
4427  NAMD_die("QM Conditional SMD is ON, but no CSMD configuration file was profided!");
4428  }
4429 
4430 #if defined(NAMD_CUDA) || defined(NAMD_HIP)
4431  // Disable various CUDA kernels if they do not fully support
4432  // or are otherwise incompatible with simulation options.
4433  if ( useCUDAdisable ) {
4434  if ( drudeOn && useCUDA2 && (bondedCUDA & 0x0001) ) {
4435  // disable CUDA kernels for spring bonds
4436  bondedCUDA &= ~0x0001;
4437  iout << iWARN << "Disabling CUDA kernel for bonds due to incompatibility with Drude oscillators.\n";
4438  }
4439  if ( accelMDOn && (accelMDdihe || accelMDdual) && useCUDA2 && (bondedCUDA & (0x0004 | 0x0020)) ) {
4440  // disable CUDA kernels for dihedrals and crossterms
4441  bondedCUDA &= ~(0x0004 | 0x0020);
4442  iout << iWARN << "Disabling CUDA kernels for dihedrals and crossterms due to incompatibility with accelerated MD options.\n";
4443  }
4444  }
4445 #endif
4446 
4447 #ifdef NAMD_AVXTILES
4448  if (avxTilesCommandLineDisable) useAVXTiles = FALSE;
4449  if (useAVXTiles) {
4450  if (alchOn || lesOn || tabulatedEnergies || drudeOn || goForcesOn ||
4451  pressureProfileOn || qmForcesOn) {
4452  useAVXTiles = FALSE;
4453  iout << iWARN << "Disabling AVX tiles optimizations due to "
4454  << "incompatible simulation params.\n";
4455  }
4456  }
4457 #else
4458  useAVXTiles = FALSE;
4459 #endif
4460 
4461 } // check_config()
4462 
4463 
4464 void SimParameters::print_config(ParseOptions &opts, ConfigList *config, char *&cwd) {
4465 
4466  StringList *current; // Pointer to config option list
4467 
4468  // Now that we have read everything, print it out so that
4469  // the user knows what is going on
4470  iout << iINFO << "SIMULATION PARAMETERS:\n";
4471  iout << iINFO << "TIMESTEP " << dt << "\n" << endi;
4472  iout << iINFO << "NUMBER OF STEPS " << N << "\n";
4473  iout << iINFO << "STEPS PER CYCLE " << stepsPerCycle << "\n";
4474  iout << endi;
4475 
4476  if ( lattice.a_p() || lattice.b_p() || lattice.c_p() ) {
4477  if ( lattice.a_p() )
4478  iout << iINFO << "PERIODIC CELL BASIS 1 " << lattice.a() << "\n";
4479  if ( lattice.b_p() )
4480  iout << iINFO << "PERIODIC CELL BASIS 2 " << lattice.b() << "\n";
4481  if ( lattice.c_p() )
4482  iout << iINFO << "PERIODIC CELL BASIS 3 " << lattice.c() << "\n";
4483  iout << iINFO << "PERIODIC CELL CENTER " << lattice.origin() << "\n";
4484  if (wrapWater) {
4485  iout << iINFO << "WRAPPING WATERS AROUND PERIODIC BOUNDARIES ON OUTPUT.\n";
4486  }
4487  if (wrapAll) {
4488  iout << iINFO << "WRAPPING ALL CLUSTERS AROUND PERIODIC BOUNDARIES ON OUTPUT.\n";
4489  }
4490  if (wrapNearest) {
4491  iout << iINFO << "WRAPPING TO IMAGE NEAREST TO PERIODIC CELL CENTER.\n";
4492  }
4493  iout << endi;
4494  }
4495 
4496  if ( CkNumPes() > 512 ) ldbUnloadOne = TRUE;
4497  if ( ldbUnloadOne || CkNumPes() > 128 ) ldbUnloadZero = TRUE;
4498 
4499  if (ldBalancer == LDBAL_NONE) {
4500  iout << iINFO << "LOAD BALANCER None\n" << endi;
4501  } else {
4502  if (ldBalancer == LDBAL_CENTRALIZED) {
4503  iout << iINFO << "LOAD BALANCER Centralized\n" << endi;
4504  } else if (ldBalancer == LDBAL_HYBRID) {
4505  iout << iINFO << "LOAD BALANCER Hybrid\n" << endi;
4506  }
4507 
4508  if (ldbStrategy == LDBSTRAT_DEFAULT) {
4509  iout << iINFO << "LOAD BALANCING STRATEGY New Load Balancers -- DEFAULT\n";
4510  } else if (ldbStrategy == LDBSTRAT_REFINEONLY) {
4511  iout << iINFO << "LOAD BALANCING STRATEGY Refinement Only\n";
4512  } else if (ldbStrategy == LDBSTRAT_COMPREHENSIVE) {
4513  iout << iINFO << "LOAD BALANCING STRATEGY Comprehensive\n";
4514  } else if (ldbStrategy == LDBSTRAT_OLD) {
4515  iout << iINFO << "LOAD BALANCING STRATEGY Old Load Balancers\n";
4516  }
4517 
4518  iout << iINFO << "LDB PERIOD " << ldbPeriod << " steps\n";
4519  iout << iINFO << "FIRST LDB TIMESTEP " << firstLdbStep << "\n";
4520  if (ldBalancer == LDBAL_HYBRID)
4521  iout << iINFO << "HYBRIDLB GROUP SIZE " << hybridGroupSize << "\n";
4522  iout << iINFO << "LAST LDB TIMESTEP " << lastLdbStep << "\n";
4523  if ( ldbRelativeGrainsize > 0. )
4524  iout << iINFO << "LDB RELATIVE GRAINSIZE " << ldbRelativeGrainsize << "\n";
4525  iout << iINFO << "LDB BACKGROUND SCALING " << ldbBackgroundScaling << "\n";
4526  iout << iINFO << "HOM BACKGROUND SCALING " << ldbHomeBackgroundScaling << "\n";
4527  if ( PMEOn ) {
4528  iout << iINFO << "PME BACKGROUND SCALING "
4529  << ldbPMEBackgroundScaling << "\n";
4530  if ( ldbUnloadPME )
4531  iout << iINFO << "REMOVING LOAD FROM PME NODES" << "\n";
4532  }
4533  if ( ldbUnloadZero ) iout << iINFO << "REMOVING LOAD FROM NODE 0\n";
4534  if ( ldbUnloadOne ) iout << iINFO << "REMOVING LOAD FROM NODE 1\n";
4535  if ( ldbUnloadOutputPEs ) iout << iINFO << "REMOVING LOAD FROM OUTPUT PES\n";
4536  iout << endi;
4537  }
4538 
4539  if ( ldbUnloadOne || CkNumPes() > 256 ) noPatchesOnOne = TRUE;
4540  if ( ldbUnloadZero || noPatchesOnOne ||
4541  CkNumPes() > 64 || ( IMDon && CkNumPes() > 8 ) ) {
4542  noPatchesOnZero = TRUE;
4543  }
4544  if ( noPatchesOnZero ) iout << iINFO << "REMOVING PATCHES FROM PROCESSOR 0\n";
4545  if ( noPatchesOnOne ) iout << iINFO << "REMOVING PATCHES FROM PROCESSOR 1\n";
4546  iout << endi;
4547 
4548 #if defined(NAMD_CUDA) || defined(NAMD_HIP) || defined(NAMD_MIC)
4549  maxSelfPart = maxPairPart = 1;
4550 #endif
4551 
4552  if (ldBalancer == LDBAL_HYBRID) {
4553  iout << iINFO << "MAX SELF PARTITIONS " << maxSelfPart << "\n"
4554  << iINFO << "MAX PAIR PARTITIONS " << maxPairPart << "\n"
4555  << iINFO << "SELF PARTITION ATOMS " << numAtomsSelf << "\n"
4556  << iINFO << "SELF2 PARTITION ATOMS " << numAtomsSelf2 << "\n"
4557  << iINFO << "PAIR PARTITION ATOMS " << numAtomsPair << "\n"
4558  << iINFO << "PAIR2 PARTITION ATOMS " << numAtomsPair2 << "\n";
4559  }
4560  iout << iINFO << "MIN ATOMS PER PATCH " << minAtomsPerPatch << "\n"
4561  << iINFO << "EMPTY PATCH LOAD " << emptyPatchLoad << " ATOMS\n"
4562  << endi;
4563 
4564  if (initialTemp < 0)
4565  {
4566  current = config->find("velocities");
4567 
4568  if (current == NULL)
4569  {
4570  current = config->find("binvelocities");
4571  }
4572 
4573  iout << iINFO << "VELOCITY FILE " << current->data << "\n";
4574  }
4575  else
4576  {
4577  iout << iINFO << "INITIAL TEMPERATURE "
4578  << initialTemp << "\n";
4579  }
4580  iout << endi;
4581 
4582  iout << iINFO << "CENTER OF MASS MOVING INITIALLY? ";
4583 
4584  if (comMove)
4585  {
4586  iout << "YES\n";
4587  }
4588  else
4589  {
4590  iout << "NO\n";
4591  }
4592  iout << endi;
4593 
4594  if ( zeroMomentum ) {
4595  iout << iINFO << "REMOVING CENTER OF MASS DRIFT DURING SIMULATION";
4596  if ( zeroMomentumAlt ) iout << " (ALT METHOD)";
4597  iout << "\n" << endi;
4598  }
4599 
4600  iout << iINFO << "DIELECTRIC "
4601  << dielectric << "\n";
4602 
4603  if ( nonbondedScaling != 1.0 )
4604  {
4605  iout << iINFO << "NONBONDED SCALING " << nonbondedScaling << "\n" << endi;
4606  }
4607  iout << iINFO << "EXCLUDE ";
4608 
4609  switch (exclude)
4610  {
4611  case NONE:
4612  iout << "NONE\n";
4613  break;
4614  case ONETWO:
4615  iout << "ONETWO\n";
4616  break;
4617  case ONETHREE:
4618  iout << "ONETHREE\n";
4619  break;
4620  case ONEFOUR:
4621  iout << "ONE-FOUR\n";
4622  break;
4623  default:
4624  iout << "SCALED ONE-FOUR\n";
4625  break;
4626  }
4627  iout << endi;
4628 
4629  if (exclude == SCALED14)
4630  {
4631  iout << iINFO << "1-4 ELECTROSTATICS SCALED BY " << scale14 << "\n";
4632  iout << iINFO << "MODIFIED 1-4 VDW PARAMETERS WILL BE USED\n" << endi;
4633  } else {
4634  iout << iWARN << "MODIFIED 1-4 VDW PARAMETERS WILL BE IGNORED\n" << endi;
4635  }
4636 
4637 #ifdef SPEC_DISABLED_VERSION
4638  if (dcdFrequency > 0) {
4639  dcdFrequency = 0;
4640  iout << iWARN << "DCD TRAJECTORY OUTPUT IS DISABLED IN SPEC RELEASE\n";
4641  }
4642 #endif
4643 
4644  if (dcdFrequency > 0)
4645  {
4646  iout << iINFO << "DCD FILENAME "
4647  << dcdFilename << "\n";
4648  iout << iINFO << "DCD FREQUENCY "
4649  << dcdFrequency << "\n";
4650  iout << iINFO << "DCD FIRST STEP "
4651  << ( ((firstTimestep + dcdFrequency)/dcdFrequency)*dcdFrequency ) << "\n";
4652  if ( dcdUnitCell ) {
4653  iout << iINFO << "DCD FILE WILL CONTAIN UNIT CELL DATA\n";
4654  }
4655  }
4656  else
4657  {
4658  iout << iINFO << "NO DCD TRAJECTORY OUTPUT\n";
4659  }
4660  iout << endi;
4661 
4662  if (xstFrequency > 0)
4663  {
4664  iout << iINFO << "XST FILENAME "
4665  << xstFilename << "\n";
4666  iout << iINFO << "XST FREQUENCY "
4667  << xstFrequency << "\n";
4668  }
4669  else
4670  {
4671  iout << iINFO << "NO EXTENDED SYSTEM TRAJECTORY OUTPUT\n";
4672  }
4673  iout << endi;
4674 
4675  if (velDcdFrequency > 0)
4676  {
4677  iout << iINFO << "VELOCITY DCD FILENAME "
4678  << velDcdFilename << "\n";
4679  iout << iINFO << "VELOCITY DCD FREQUENCY "
4680  << velDcdFrequency << "\n";
4681  iout << iINFO << "VELOCITY DCD FIRST STEP "
4682  << ( ((firstTimestep + velDcdFrequency)/velDcdFrequency)*velDcdFrequency ) << "\n";
4683  }
4684  else
4685  {
4686  iout << iINFO << "NO VELOCITY DCD OUTPUT\n";
4687  }
4688  iout << endi;
4689 
4690  if (forceDcdFrequency > 0)
4691  {
4692  iout << iINFO << "FORCE DCD FILENAME "
4693  << forceDcdFilename << "\n";
4694  iout << iINFO << "FORCE DCD FREQUENCY "
4695  << forceDcdFrequency << "\n";
4696  iout << iINFO << "FORCE DCD FIRST STEP "
4697  << ( ((firstTimestep + forceDcdFrequency)/forceDcdFrequency)*forceDcdFrequency ) << "\n";
4698  }
4699  else
4700  {
4701  iout << iINFO << "NO FORCE DCD OUTPUT\n";
4702  }
4703  iout << endi;
4704 
4705  iout << iINFO << "OUTPUT FILENAME "
4706  << outputFilename << "\n" << endi;
4707  if (binaryOutput)
4708  {
4709  iout << iINFO << "BINARY OUTPUT FILES WILL BE USED\n" << endi;
4710  }
4711 #ifdef MEM_OPT_VERSION
4712  if(!binaryOutput){
4713  iout << iWARN <<"SINCE MEMORY OPTIMIZED VERSION IS USED, OUTPUT IN TEXT FORMAT IS DISABLED!\n" << endi;
4714  binaryOutput = TRUE;
4715  }
4716 #endif
4717 
4718  if (! restartFrequency)
4719  {
4720  iout << iINFO << "NO RESTART FILE\n";
4721  }
4722  else
4723  {
4724  iout << iINFO << "RESTART FILENAME "
4725  << restartFilename << "\n";
4726  iout << iINFO << "RESTART FREQUENCY "
4727  << restartFrequency << "\n";
4728  if (restartSave) {
4729  iout << iINFO << "RESTART FILES WILL NOT BE OVERWRITTEN\n";
4730  }
4731  if (restartSaveDcd) {
4732  iout << iINFO << "DCD FILE WILL BE SPLIT WHEN RESTART FILES ARE WRITTEN\n";
4733  }
4734 
4735  if (binaryRestart)
4736  {
4737  iout << iINFO << "BINARY RESTART FILES WILL BE USED\n";
4738  }
4739  }
4740  iout << endi;
4741 
4742  if (switchingActive)
4743  {
4744  iout << iINFO << "SWITCHING ACTIVE\n";
4745  if ( vdwForceSwitching ) {
4746  iout << iINFO << "VDW FORCE SWITCHING ACTIVE\n";
4747  }
4748  if ( martiniSwitching ) {
4749  iout << iINFO << "MARTINI RESIDUE-BASED COARSE-GRAIN SWITCHING ACTIVE\n";
4750  }
4751  iout << iINFO << "SWITCHING ON "
4752  << switchingDist << "\n";
4753  iout << iINFO << "SWITCHING OFF "
4754  << cutoff << "\n";
4755  }
4756  else
4757  {
4758  iout << iINFO << "CUTOFF "
4759  << cutoff << "\n";
4760  }
4761 
4762  iout << iINFO << "PAIRLIST DISTANCE " << pairlistDist << "\n";
4763  iout << iINFO << "PAIRLIST SHRINK RATE " << pairlistShrink << "\n";
4764  iout << iINFO << "PAIRLIST GROW RATE " << pairlistGrow << "\n";
4765  iout << iINFO << "PAIRLIST TRIGGER " << pairlistTrigger << "\n";
4766  iout << iINFO << "PAIRLISTS PER CYCLE " << pairlistsPerCycle << "\n";
4767  if ( outputPairlists )
4768  iout << iINFO << "PAIRLIST OUTPUT STEPS " << outputPairlists << "\n";
4769  iout << endi;
4770 
4771  if ( pairlistMinProcs > 1 )
4772  iout << iINFO << "REQUIRING " << pairlistMinProcs << " PROCESSORS FOR PAIRLISTS\n";
4773  usePairlists = ( CkNumPes() >= pairlistMinProcs );
4774 
4775 #ifdef OPENATOM_VERSION
4776 if ( openatomOn )
4777 {
4778  iout << iINFO << "OPENATOM QM/MM CAR-PARINELLO ACTIVE\n";
4779  iout << iINFO << "OPENATOM CONFIG FILE: " << openatomConfig << "\n";
4780  iout << iINFO << "OPENATOM STRUCT FILE: " << openatomStruct << "\n";
4781  iout << iINFO << "OPENATOM PDB FILE: " << openatomPDB << "\n";
4782 }
4783 #endif // OPENATOM_VERSION
4784 
4785  // FB - FEP and TI are now dependent on pairlists - disallow usePairlists=0
4786  if ( (alchOn) && (!usePairlists)) {
4787  NAMD_die("Sorry, Alchemical simulations require pairlists to be enabled\n");
4788  }
4789 #if defined(NAMD_CUDA) || defined(NAMD_HIP)
4790  if ( ! usePairlists ) {
4791  usePairlists = 1;
4792  iout << iINFO << "CUDA ACCELERATION REQUIRES PAIRLISTS\n";
4793  }
4794 #endif
4795 
4796  iout << iINFO << "PAIRLISTS " << ( usePairlists ? "ENABLED" : "DISABLED" )
4797  << "\n" << endi;
4798 
4799  iout << iINFO << "MARGIN " << margin << "\n";
4800  if ( margin > 4.0 ) {
4801  iout << iWARN << "MARGIN IS UNUSUALLY LARGE AND WILL LOWER PERFORMANCE\n";
4802  BigReal f = patchDimension/(patchDimension-margin);
4803  f *= f*f;
4804  iout << iWARN << "MARGIN INCREASED PATCH VOLUME BY A FACTOR OF " << f << "\n";
4805  }
4806 
4807  if ( splitPatch == SPLIT_PATCH_HYDROGEN ) {
4808  iout << iINFO << "HYDROGEN GROUP CUTOFF " << hgroupCutoff << "\n";
4809  }
4810 
4811  iout << iINFO << "PATCH DIMENSION "
4812  << patchDimension << "\n";
4813 
4814  iout << endi;
4815 
4816  if (outputEnergies != 1)
4817  {
4818  iout << iINFO << "ENERGY OUTPUT STEPS "
4819  << outputEnergies << "\n";
4820  iout << endi;
4821  }
4822 
4823  if (mergeCrossterms) {
4824  iout << iINFO << "CROSSTERM ENERGY INCLUDED IN DIHEDRAL\n" << endi;
4825  }
4826 
4827  if (outputMomenta != 0)
4828  {
4829  iout << iINFO << "MOMENTUM OUTPUT STEPS "
4830  << outputMomenta << "\n";
4831  iout << endi;
4832  }
4833 
4834  if (outputTiming != 0)
4835  {
4836  iout << iINFO << "TIMING OUTPUT STEPS "
4837  << outputTiming << "\n";
4838  iout << endi;
4839  }
4840 
4841  if (outputCudaTiming != 0)
4842  {
4843  iout << iINFO << "CUDA TIMING OUTPUT STEPS "
4844  << outputCudaTiming << "\n";
4845  iout << endi;
4846  }
4847 
4848  if (outputPressure != 0)
4849  {
4850  iout << iINFO << "PRESSURE OUTPUT STEPS "
4851  << outputPressure << "\n";
4852  iout << endi;
4853  }
4854 
4855  if (fixedAtomsOn)
4856  {
4857  iout << iINFO << "FIXED ATOMS ACTIVE\n";
4858  if ( fixedAtomsForces )
4859  iout << iINFO << "FORCES BETWEEN FIXED ATOMS ARE CALCULATED\n";
4860  iout << endi;
4861  }
4862 
4863  if (constraintsOn)
4864  {
4865  iout << iINFO << "HARMONIC CONSTRAINTS ACTIVE\n";
4866 
4867  iout << iINFO << "HARMONIC CONS EXP "
4868  << constraintExp << "\n";
4869 
4870  if (constraintScaling != 1.0) {
4871  iout << iINFO << "HARMONIC CONS SCALING "
4872  << constraintScaling << "\n";
4873  }
4874 
4875  //****** BEGIN selective restraints (X,Y,Z) changes
4876 
4877  if (selectConstraintsOn) {
4878  iout << iINFO << "SELECTED CARTESIAN COMPONENTS OF HARMONIC RESTRAINTS ACTIVE\n";
4879 
4880  if (constrXOn)
4881  iout << iINFO << "RESTRAINING X-COMPONENTS OF CARTESIAN COORDINATES!\n";
4882 
4883  if (constrYOn)
4884  iout << iINFO << "RESTRAINING Y-COMPONENTS OF CARTESIAN COORDINATES!\n";
4885 
4886  if (constrZOn)
4887  iout << iINFO << "RESTRAINING Z-COMPONENTS OF CARTESIAN COORDINATES!\n";
4888  }
4889  //****** END selective restraints (X,Y,Z) changes
4890 
4891  if (sphericalConstraintsOn) {
4892  iout << iINFO << "SPHERICAL HARMONIC CONSTRAINTS ACTIVE\n";
4893  iout << iINFO << "RESTRAINING DISTANCE TO " << sphericalConstrCenter <<"\n";
4894  }
4895  iout << endi;
4896 
4897  //****** BEGIN moving constraints changes
4898 
4899  if (movingConstraintsOn) {
4900  iout << iINFO << "MOVING HARMONIC CONSTRAINTS ACTIVE\n";
4901 
4902  iout << iINFO << "MOVING CONSTRAINT VELOCITY "
4903  << movingConsVel << " ANGSTROM/TIMESTEP\n";
4904 
4905  iout << iINFO << "ALL CONSTRAINED ATOMS WILL MOVE\n";
4906  }
4907  //****** END moving constraints changes
4908  iout << endi;
4909 
4910  //****** BEGIN rotating constraints changes
4911 
4912  if (rotConstraintsOn) {
4913  iout << iINFO << "ROTATING HARMONIC CONSTRAINTS ACTIVE\n";
4914 
4915  iout << iINFO << "AXIS OF ROTATION "
4916  << rotConsAxis << "\n";
4917 
4918  iout << iINFO << "PIVOT OF ROTATION "
4919  << rotConsPivot << "\n";
4920 
4921  iout << iINFO << "ROTATING CONSTRAINT VELOCITY "
4922  << rotConsVel << " DEGREES/TIMESTEP\n";
4923  }
4924  iout << endi;
4925  //****** END rotating constraints changes
4926  }
4927 
4928  // moving drag
4929  if (movDragOn) {
4930  iout << iINFO << "MOVING DRAG ACTIVE.\n";
4931 
4932  iout << iINFO << "MOVING DRAG MAIN PDB FILE "
4933  << movDragFile << "\n";
4934 
4935  iout << iINFO << "MOVING DRAG GLOBAL VELOCITY (A/step) "
4936  << movDragGlobVel << "\n";
4937 
4938  iout << iINFO << "MOVING DRAG LINEAR VELOCITY FILE "
4939  << movDragVelFile << "\n";
4940 
4941  iout << endi;
4942  }
4943 
4944  // rotating drag
4945  if (rotDragOn) {
4946  iout << iINFO << "ROTATING DRAG ACTIVE.\n";
4947 
4948  iout << iINFO << "ROTATING DRAG MAIN PDB FILE "
4949  << rotDragFile << "\n";
4950 
4951  iout << iINFO << "ROTATING DRAG AXIS FILE "
4952  << rotDragAxisFile << "\n";
4953 
4954  iout << iINFO << "ROTATING DRAG PIVOT POINT FILE "
4955  << rotDragPivotFile << "\n";
4956 
4957  iout << iINFO << "ROTATING DRAG GLOBAL ANGULAR VELOCITY (deg/step) "
4958  << rotDragGlobVel << "\n";
4959 
4960  iout << iINFO << "ROTATING DRAG ANGULAR VELOCITY FILE "
4961  << rotDragVelFile << "\n";
4962 
4963  iout << endi;
4964  }
4965 
4966 
4967  // "constant" torque
4968  if (consTorqueOn) {
4969  iout << iINFO << "\"CONSTANT\" TORQUE ACTIVE.\n";
4970 
4971  iout << iINFO << "\"CONSTANT\" TORQUE MAIN PDB FILE "
4972  << consTorqueFile << "\n";
4973 
4974  iout << iINFO << "\"CONSTANT\" TORQUE AXIS FILE "
4975  << consTorqueAxisFile << "\n";
4976 
4977  iout << iINFO << "\"CONSTANT\" TORQUE PIVOT POINT FILE "
4978  << consTorquePivotFile << "\n";
4979 
4980  iout << iINFO << "\"CONSTANT\" TORQUE GLOBAL VALUE (Kcal/(mol*A^2)) "
4981  << consTorqueGlobVal << "\n";
4982 
4983  iout << iINFO << "\"CONSTANT\" TORQUE DACTORS FILE "
4984  << consTorqueValFile << "\n";
4985 
4986  iout << endi;
4987  }
4988 
4989  if (mgridforceOn) {
4990  iout << iINFO << "GRID FORCE ACTIVE\n";
4991  iout << iINFO << " Please include this reference in published work using\n";
4992  iout << iINFO << " the Gridforce module of NAMD: David Wells, Volha Abramkina,\n";
4993  iout << iINFO << " and Aleksei Aksimentiev, J. Chem. Phys. 127:125101-10 (2007).\n";
4994  print_mgrid_params();
4995  }
4996 
4997  //****** BEGIN SMD constraints changes
4998 
4999  if (SMDOn) {
5000  iout << iINFO << "SMD ACTIVE\n";
5001 
5002  iout << iINFO << "SMD VELOCITY "
5003  << SMDVel << " ANGSTROM/TIMESTEP\n";
5004 
5005  iout << iINFO << "SMD DIRECTION "
5006  << SMDDir << "\n";
5007 
5008  iout << iINFO << "SMD K "
5009  << SMDk << "\n";
5010 
5011  iout << iINFO << "SMD K2 "
5012  << SMDk2 << "\n";
5013 
5014  iout << iINFO << "SMD OUTPUT FREQUENCY "
5015  << SMDOutputFreq << " TIMESTEPS\n";
5016 
5017  iout << iINFO << "SMD FILE " << SMDFile << "\n";
5018 
5019  iout << endi;
5020  }
5021 
5022  //****** END SMD constraints changes
5023 
5024  if (TMDOn) {
5025  iout << iINFO << "TMD ACTIVE BETWEEN STEPS " << TMDFirstStep
5026  << " and " << TMDLastStep << "\n";
5027  iout << iINFO << "TMD K " << TMDk << "\n";
5028  iout << iINFO << "TMD FILE " << TMDFile << "\n";
5029  iout << iINFO << "TMD OUTPUT FREQUENCY " << TMDOutputFreq << "\n";
5030  if (TMDInitialRMSD) {
5031  iout << iINFO << "TMD TARGET RMSD AT FIRST STEP " << TMDInitialRMSD << "\n";
5032  } else {
5033  iout << iINFO << "TMD TARGET RMSD AT FIRST STEP COMPUTED FROM INITIAL COORDINATES\n";
5034  }
5035  iout << iINFO << "TMD TARGET RMSD AT FINAL STEP " << TMDFinalRMSD << "\n";
5036  iout << endi;
5037  }
5038 
5039  if (symmetryOn) {
5040  if (symmetryLastStep == -1){
5041  iout << iINFO << "SYMMETRY RESTRAINTS ACTIVE BETWEEN STEPS " << symmetryFirstStep << " and " << "INFINITY" << "\n";
5042  }
5043  else{
5044  iout << iINFO << "SYMMETRY RESTRAINTS ACTIVE BETWEEN STEPS " << symmetryFirstStep << " and " << symmetryLastStep << "\n";
5045  }
5046  // iout << iINFO << "SYMMETRY FILE " << symmetryFile << "\n";
5047 
5048  current = config->find("symmetryFile");
5049  for ( ; current; current = current->next ) {
5050  iout << iINFO << "SYMMETRY FILE " << current->data << "\n";
5051  }
5052 
5053  current = config->find("symmetryMatrixFile");
5054  for ( ; current; current = current->next ) {
5055  iout << iINFO << "SYMMETRY MATRIX FILE " << current->data << "\n";
5056  }
5057  iout << iINFO << "SYMMETRY FORCE CONSTANT " << symmetryk << "\n";
5058  if (symmetryScaleForces){
5059  iout << iINFO << "SYMMETRY SCALE FORCES ON\n";
5060  }
5061  iout << iINFO << "SYMMETRY FIRST FULL STEP " << symmetryFirstFullStep << "\n";
5062  if (symmetryLastFullStep == -1){
5063  iout << iINFO << "SYMMETRY LAST FULL STEP " << "INFINITY" << "\n";
5064  //iout << iINFO << "FULL SYMMETRY FORCE BETWEEN STEPS " << symmetryFirstFullStep << " and " << "INFINITY" << "\n";
5065  }
5066  else {
5067  iout << iINFO << "SYMMETRY LAST FULL STEP " << symmetryLastFullStep << "\n";
5068  // iout << iINFO << "FULL SYMMETRY FORCE BETWEEN STEPS " << symmetryFirstFullStep << " and " << symmetryLastFullStep << "\n";
5069  }
5070 
5071  iout << endi;
5072  }
5073 //Modifications for alchemical fep
5074 // Alchemical FEP status
5075 
5076 // current = config->find("alchOutFile");
5077  if (alchFepOn)
5078  {
5079  iout << iINFO << "ALCHEMICAL FEP ON\n";
5080  iout << iINFO << "FEP CURRENT LAMBDA VALUE "
5081  << alchLambda << "\n";
5082  iout << iINFO << "FEP COMPARISON LAMBDA VALUE "
5083  << alchLambda2 << "\n";
5084  if (alchLambdaIDWS >= 0.) {
5085  iout << iINFO << "FEP ALTERNATE COMPARISON LAMBDA VALUE "
5086  << alchLambdaIDWS << "\n";
5087  }
5088  if (alchLambdaFreq > 0) {
5089  iout << iINFO << "FEP CURRENT LAMBDA VALUE SET TO INCREASE IN EVERY "
5090  << alchLambdaFreq << " STEPS\n";
5091  }
5092  if (!alchDecouple) {
5093  iout << iINFO << "FEP INTRA-ALCHEMICAL NON-BONDED INTERACTIONS WILL BE "
5094  << "DECOUPLED\n";
5095  }else{
5096  iout << iINFO << "FEP INTRA-ALCHEMICAL NON-BONDED INTERACTIONS WILL BE "
5097  << "RETAINED\n";
5098  }
5099  if (alchBondDecouple) {
5100  iout << iINFO << "FEP INTRA-ALCHEMICAL BONDED INTERACTIONS WILL BE "
5101  << "DECOUPLED\n";
5102  }else{
5103  iout << iINFO << "FEP INTRA-ALCHEMICAL BONDED INTERACTIONS WILL BE "
5104  << "RETAINED\n";
5105  }
5106  if (alchWCAOn) {
5107  iout << iINFO << "FEP WEEKS-CHANDLER-ANDERSEN (WCA) VDW DECOUPLING "
5108  << "ACTIVE\n";
5109  } else {
5110  iout << iINFO << "FEP VDW SHIFTING COEFFICIENT "
5111  << alchVdwShiftCoeff << "\n";
5112  }
5113  iout << iINFO << "FEP ELEC. ACTIVE FOR ANNIHILATED "
5114  << "PARTICLES BETWEEN LAMBDA = 0 AND LAMBDA = "
5115  << (1 - alchElecLambdaStart) << "\n";
5116  iout << iINFO << "FEP ELEC. ACTIVE FOR EXNIHILATED "
5117  << "PARTICLES BETWEEN LAMBDA = "
5118  << alchElecLambdaStart << " AND LAMBDA = 1\n";
5119  if (alchWCAOn) {
5120  iout << iINFO << "FEP VDW-REPU. ACTIVE FOR ANNIHILATED PARTICLES "
5121  << "BETWEEN LAMBDA = " << (1 - alchRepLambdaEnd) << " AND LAMBDA "
5122  << "= 1\n";
5123  iout << iINFO << "FEP VDW-REPU. ACTIVE FOR EXNIHILATED PARTICLES "
5124  << "BETWEEN LAMBDA = 0 AND LAMBDA " << alchRepLambdaEnd << "\n";
5125  iout << iINFO << "FEP VDW-ATTR. ACTIVE FOR ANNIHILATED PARTICLES "
5126  << "BETWEEN LAMBDA = " << (1 - alchVdwLambdaEnd) << " AND LAMBDA = "
5127  << (1 - alchRepLambdaEnd) << "\n";
5128  iout << iINFO << "FEP VDW-ATTR. ACTIVE FOR EXNIHILATED PARTICLES "
5129  << "BETWEEN LAMBDA = " << alchRepLambdaEnd << " AND LAMBDA = "
5130  << alchVdwLambdaEnd << "\n";
5131  } else {
5132  iout << iINFO << "FEP VDW ACTIVE FOR ANNIHILATED "
5133  << "PARTICLES BETWEEN LAMBDA = "
5134  << (1 - alchVdwLambdaEnd) << " AND LAMBDA = 1\n";
5135  iout << iINFO << "FEP VDW ACTIVE FOR EXNIHILATED "
5136  << "PARTICLES BETWEEN LAMBDA = 0 AND LAMBDA = "
5137  << alchVdwLambdaEnd << "\n";
5138  }
5139  iout << iINFO << "FEP BOND ACTIVE FOR ANNIHILATED "
5140  << "PARTICLES BETWEEN LAMBDA = "
5141  << (1 - alchBondLambdaEnd) << " AND LAMBDA = 1\n";
5142  iout << iINFO << "FEP BOND ACTIVE FOR EXNIHILATED "
5143  << "PARTICLES BETWEEN LAMBDA = 0 AND LAMBDA = "
5144  << alchBondLambdaEnd << "\n";
5145  }
5146 //fepe
5147 
5148  if (alchThermIntOn)
5149  {
5150  iout << iINFO << "THERMODYNAMIC INTEGRATION (TI) ON\n";
5151  iout << iINFO << "TI LAMBDA VALUE "
5152  << alchLambda << "\n";
5153  if (alchLambdaFreq > 0) {
5154  iout << iINFO << "TI COMPARISON LAMBDA VALUE "
5155  << alchLambda2 << "\n";
5156  iout << iINFO << "TI CURRENT LAMBDA VALUE SET TO INCREASE IN EVERY "
5157  << alchLambdaFreq << " STEPS\n";
5158  }
5159  if (!alchDecouple) {
5160  iout << iINFO << "TI INTRA-ALCHEMICAL NON-BONDED INTERACTIONS WILL BE "
5161  << "DECOUPLED\n";
5162  }else{
5163  iout << iINFO << "TI INTRA-ALCHEMICAL NON-BONDED INTERACTIONS WILL BE "
5164  << "RETAINED\n";
5165  }
5166  if (alchBondDecouple) {
5167  iout << iINFO << "TI INTRA-ALCHEMICAL BONDED INTERACTIONS WILL BE "
5168  << "DECOUPLED\n";
5169  }else{
5170  iout << iINFO << "TI INTRA-ALCHEMICAL BONDED INTERACTIONS WILL BE "
5171  << "RETAINED\n";
5172  }
5173  iout << iINFO << "TI VDW SHIFTING COEFFICIENT "
5174  << alchVdwShiftCoeff << "\n";
5175  iout << iINFO << "TI ELEC. ACTIVE FOR ANNIHILATED "
5176  << "PARTICLES BETWEEN LAMBDA = 0 AND LAMBDA = "
5177  << (1 - alchElecLambdaStart) << "\n";
5178  iout << iINFO << "TI ELEC. ACTIVE FOR EXNIHILATED "
5179  << "PARTICLES BETWEEN LAMBDA = "
5180  << alchElecLambdaStart << " AND LAMBDA = 1\n";
5181  iout << iINFO << "TI VDW ACTIVE FOR ANNIHILATED "
5182  << "PARTICLES BETWEEN LAMBDA = "
5183  << (1 - alchVdwLambdaEnd) << " AND LAMBDA = 1\n";
5184  iout << iINFO << "TI VDW ACTIVE FOR EXNIHILATED "