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