Update default behavior and documentation for Drude
[namd.git] / src / SimParameters.C
1 /**
2 ***  Copyright (c) 1995, 1996, 1997, 1998, 1999, 2000 by
3 ***  The Board of Trustees of the University of Illinois.
4 ***  All rights reserved.
5 **/
6
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
14 /** \file SimParameters.C
15  * SimParameters is just a glorified structure to hold the global
16  * static simulation parameters such as timestep size, cutoff, etc. that
17  * are read in from the configuration file.
18  **/
19
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
80 #include "ComputeNonbondedMICKernel.h"
81
82 //#ifdef NAMD_CUDA
83 //bool one_cuda_device_per_node();
84 //#endif
85 #include "DeviceCUDA.h"
86 #ifdef NAMD_CUDA
87 #ifdef WIN32
88 #define __thread __declspec(thread)
89 #endif
90 extern __thread DeviceCUDA *deviceCUDA;
91 #endif
92
93 //#define DEBUGM
94 #include "Debug.h"
95
96 #define XXXBIGREAL 1.0e32
97
98
99 char* SimParameters::getfromparseopts(const char* name, char *outbuf) {
100   if ( parseopts ) return parseopts->getfromptr(name,outbuf);
101   else return 0;
102 }
103
104 int SimParameters::istrueinparseopts(const char* name) {
105   if ( parseopts ) return parseopts->istruefromptr(name);
106   else return -1;
107 }
108
109 int SimParameters::issetinparseopts(const char* name) {
110   if ( parseopts ) return parseopts->issetfromptr(name);
111   else return -1;
112 }
113
114
115 /************************************************************************/
116 /*                  */
117 /*      FUNCTION initialize_config_data      */
118 /*                  */
119 /*  This function is used by the master process to populate the     */
120 /*   simulation parameters from a ConfigList object that is passed to   */
121 /*   it.  Each parameter is checked to make sure that it has a value    */
122 /*   that makes sense, and that it doesn't conflict with any other      */
123 /*   values that have been given.          */
124 /*                  */
125 /************************************************************************/
126
127 void SimParameters::initialize_config_data(ConfigList *config, char *&cwd)
128
129 {
130
131    parseopts = new ParseOptions;   //  Object to check consistency of config file
132    ParseOptions &opts = *parseopts;
133
134    config_parser(opts);
135
136    ///////////////////////////////// check the internal consistency
137    if (!opts.check_consistency()) 
138    {
139       NAMD_die("Internal error in configuration file parser");
140    }
141
142    // Now, feed the object with the actual configuration options through the
143    // ParseOptions file and make sure everything is OK
144    if (!opts.set(*config)) 
145    {
146       NAMD_die("ERROR(S) IN THE CONFIGURATION FILE");
147    }
148
149    //// now do more logic stuff that can't be done by the ParseOptions object
150
151    check_config(opts,config,cwd);
152
153    print_config(opts,config,cwd);
154
155 }
156
157 /************************************************************************/
158 /*                                                                      */
159 /*      FUNCTION scriptSet                                              */
160 /*                                                                      */
161 /************************************************************************/
162
163 int atobool(const char *s) {
164   return ( (! strncasecmp(s,"yes",8)) ||
165            (! strncasecmp(s,"on",8)) || 
166            (! strncasecmp(s,"true",8)) );
167 };
168                          
169 void SimParameters::scriptSet(const char *param, const char *value) {
170
171   if ( CkMyRank() ) return;
172
173 #define MAX_SCRIPT_PARAM_SIZE 128
174 #define SCRIPT_PARSE_BOOL(NAME,VAR) { if ( ! strncasecmp(param,(NAME),MAX_SCRIPT_PARAM_SIZE) ) { (VAR) = atobool(value); return; } }
175 #define SCRIPT_PARSE_INT(NAME,VAR) { if ( ! strncasecmp(param,(NAME),MAX_SCRIPT_PARAM_SIZE) ) { (VAR) = atoi(value); return; } }
176 #define SCRIPT_PARSE_FLOAT(NAME,VAR) { if ( ! strncasecmp(param,(NAME),MAX_SCRIPT_PARAM_SIZE) ) { (VAR) = atof(value); return; } }
177 #define SCRIPT_PARSE_MOD_FLOAT(NAME,VAR,MOD) { if ( ! strncasecmp(param,(NAME),MAX_SCRIPT_PARAM_SIZE) ) { (VAR) = atof(value) MOD; return; } }
178 #define SCRIPT_PARSE_VECTOR(NAME,VAR) { if ( ! strncasecmp(param,(NAME),MAX_SCRIPT_PARAM_SIZE) ) { (VAR).set(value); return; } }
179 #define SCRIPT_PARSE_STRING(NAME,VAR) { if ( ! strncasecmp(param,(NAME),MAX_SCRIPT_PARAM_SIZE) ) { strcpy(VAR,value); return; } }
180
181   SCRIPT_PARSE_FLOAT("scriptArg1",scriptArg1)
182   SCRIPT_PARSE_FLOAT("scriptArg2",scriptArg2)
183   SCRIPT_PARSE_FLOAT("scriptArg3",scriptArg3)
184   SCRIPT_PARSE_FLOAT("scriptArg4",scriptArg4)
185   SCRIPT_PARSE_FLOAT("scriptArg5",scriptArg5)
186   SCRIPT_PARSE_INT("scriptIntArg1",scriptIntArg1)
187   SCRIPT_PARSE_INT("scriptIntArg2",scriptIntArg2)
188   SCRIPT_PARSE_STRING("scriptStringArg1",scriptStringArg1)
189   SCRIPT_PARSE_STRING("scriptStringArg2",scriptStringArg2)
190   SCRIPT_PARSE_INT("numsteps",N)
191   if ( ! strncasecmp(param,"firsttimestep",MAX_SCRIPT_PARAM_SIZE) ) {
192     N = firstTimestep = atoi(value); return;
193   }
194   SCRIPT_PARSE_FLOAT("reassignTemp",reassignTemp)
195   SCRIPT_PARSE_FLOAT("rescaleTemp",rescaleTemp)
196   SCRIPT_PARSE_BOOL("velocityQuenching",minimizeOn)
197   SCRIPT_PARSE_BOOL("maximumMove",maximumMove)
198   // SCRIPT_PARSE_BOOL("Langevin",langevinOn)
199   if ( ! strncasecmp(param,"Langevin",MAX_SCRIPT_PARAM_SIZE) ) {
200     langevinOn = atobool(value);
201     if ( langevinOn && ! langevinOnAtStartup ) {
202       NAMD_die("Langevin must be enabled at startup to disable and re-enable in script.");
203     }
204     return;
205   }
206   SCRIPT_PARSE_FLOAT("langevinTemp",langevinTemp)
207   SCRIPT_PARSE_BOOL("langevinBAOAB",langevin_useBAOAB) // [!!] Use the BAOAB integrator or not
208   SCRIPT_PARSE_FLOAT("loweAndersenTemp",loweAndersenTemp) // BEGIN LA, END LA
209   SCRIPT_PARSE_BOOL("stochRescale",stochRescaleOn)
210   SCRIPT_PARSE_FLOAT("stochRescaleTemp",stochRescaleTemp)
211   SCRIPT_PARSE_FLOAT("initialTemp",initialTemp)
212   SCRIPT_PARSE_BOOL("useGroupPressure",useGroupPressure)
213   SCRIPT_PARSE_BOOL("useFlexibleCell",useFlexibleCell)
214   SCRIPT_PARSE_BOOL("useConstantArea",useConstantArea)
215   SCRIPT_PARSE_BOOL("fixCellDims",fixCellDims)
216   SCRIPT_PARSE_BOOL("fixCellDimX",fixCellDimX)
217   SCRIPT_PARSE_BOOL("fixCellDimY",fixCellDimY)
218   SCRIPT_PARSE_BOOL("fixCellDimZ",fixCellDimZ)
219   SCRIPT_PARSE_BOOL("useConstantRatio",useConstantRatio)
220   SCRIPT_PARSE_BOOL("LangevinPiston",langevinPistonOn)
221   SCRIPT_PARSE_MOD_FLOAT("LangevinPistonTarget",
222                         langevinPistonTarget,/PRESSUREFACTOR)
223   SCRIPT_PARSE_FLOAT("LangevinPistonPeriod",langevinPistonPeriod)
224   SCRIPT_PARSE_FLOAT("LangevinPistonDecay",langevinPistonDecay)
225   SCRIPT_PARSE_FLOAT("LangevinPistonTemp",langevinPistonTemp)
226   SCRIPT_PARSE_MOD_FLOAT("SurfaceTensionTarget",
227                         surfaceTensionTarget,*(100.0/PRESSUREFACTOR))
228   SCRIPT_PARSE_BOOL("BerendsenPressure",berendsenPressureOn)
229   SCRIPT_PARSE_MOD_FLOAT("BerendsenPressureTarget",
230                         berendsenPressureTarget,/PRESSUREFACTOR)
231   SCRIPT_PARSE_MOD_FLOAT("BerendsenPressureCompressibility",
232                         berendsenPressureCompressibility,*PRESSUREFACTOR)
233   SCRIPT_PARSE_FLOAT("BerendsenPressureRelaxationTime",
234                                 berendsenPressureRelaxationTime)
235   SCRIPT_PARSE_FLOAT("constraintScaling",constraintScaling)
236   SCRIPT_PARSE_FLOAT("consForceScaling",consForceScaling)
237   SCRIPT_PARSE_BOOL("drudeHardWall",drudeHardWallOn)
238   SCRIPT_PARSE_FLOAT("drudeBondConst",drudeBondConst)
239   SCRIPT_PARSE_FLOAT("drudeBondLen",drudeBondLen)
240   SCRIPT_PARSE_STRING("outputname",outputFilename)
241   SCRIPT_PARSE_INT("outputEnergies",outputEnergies)
242   SCRIPT_PARSE_STRING("restartname",restartFilename)
243   SCRIPT_PARSE_INT("DCDfreq",dcdFrequency)
244   if ( ! strncasecmp(param,"DCDfile",MAX_SCRIPT_PARAM_SIZE) ) { 
245     close_dcdfile();  // *** implemented in Output.C ***
246     strcpy(dcdFilename,value);
247     return;
248   }
249   if ( ! strncasecmp(param,"velDCDfile",MAX_SCRIPT_PARAM_SIZE) ) { 
250     close_veldcdfile();  // *** implemented in Output.C ***
251     strcpy(velDcdFilename,value);
252     return;
253   }
254   SCRIPT_PARSE_STRING("tclBCArgs",tclBCArgs)
255   SCRIPT_PARSE_VECTOR("eField",eField)
256   SCRIPT_PARSE_FLOAT("eFieldFreq",eFieldFreq)
257   SCRIPT_PARSE_FLOAT("eFieldPhase",eFieldPhase) 
258   SCRIPT_PARSE_FLOAT("accelMDE",accelMDE) 
259   SCRIPT_PARSE_FLOAT("accelMDalpha",accelMDalpha) 
260   SCRIPT_PARSE_FLOAT("accelMDTE",accelMDTE) 
261   SCRIPT_PARSE_FLOAT("accelMDTalpha",accelMDTalpha) 
262   SCRIPT_PARSE_FLOAT("accelMDGSigma0P",accelMDGSigma0P) 
263   SCRIPT_PARSE_FLOAT("accelMDGSigma0D",accelMDGSigma0D) 
264   SCRIPT_PARSE_STRING("accelMDGRestartFile",accelMDGRestartFile)
265   SCRIPT_PARSE_VECTOR("stirAxis",stirAxis)
266   SCRIPT_PARSE_VECTOR("stirPivot",stirPivot)
267   if ( ! strncasecmp(param,"mgridforcescale",MAX_SCRIPT_PARAM_SIZE) ) {
268     NAMD_die("Can't yet modify mgridforcescale in a script");
269     return;
270   }
271   if ( ! strncasecmp(param,"mgridforcevoff",MAX_SCRIPT_PARAM_SIZE) ) {
272     NAMD_die("Can't yet modify mgridforcevoff in a script");
273     return;
274   }
275    
276   if ( ! strncasecmp(param,"fixedatoms",MAX_SCRIPT_PARAM_SIZE) ) {
277     if ( ! fixedAtomsOn )
278       NAMD_die("FixedAtoms may not be enabled in a script.");
279     if ( ! fixedAtomsForces )
280       NAMD_die("To use fixedAtoms in script first use fixedAtomsForces yes.");
281     fixedAtomsOn = atobool(value);
282     return;
283   }
284
285 //fepb
286   if ( ! strncasecmp(param,"alch",MAX_SCRIPT_PARAM_SIZE) ) {
287     alchOn = atobool(value);
288     if ( alchOn && ! alchOnAtStartup ) {
289        NAMD_die("Alchemy must be enabled at startup to disable and re-enable in script.");
290     }
291     alchFepOn = alchOn && alchFepOnAtStartup;
292     alchThermIntOn = alchOn && alchThermIntOnAtStartup;
293     ComputeNonbondedUtil::select();
294     if ( PMEOn ) ComputePmeUtil::select();
295     return;
296   }
297   SCRIPT_PARSE_INT("alchEquilSteps",alchEquilSteps)
298
299   if ( ! strncasecmp(param,"alchRepLambda",MAX_SCRIPT_PARAM_SIZE) ) {
300     alchRepLambda = atof(value);
301     alchFepWCARepuOn = true;
302     alchFepWCADispOn = false;
303     alchFepElecOn    = false;
304     ComputeNonbondedUtil::select();
305     return;
306   }
307
308   if ( ! strncasecmp(param,"alchDispLambda",MAX_SCRIPT_PARAM_SIZE) ) {
309     alchDispLambda = atof(value);
310     alchFepWCARepuOn = false;
311     alchFepWCADispOn = true;
312     alchFepElecOn    = false;
313     ComputeNonbondedUtil::select();
314     return;
315   }
316
317   if ( ! strncasecmp(param,"alchElecLambda",MAX_SCRIPT_PARAM_SIZE) ) {
318     alchElecLambda = atof(value);
319     alchFepWCARepuOn = false;
320     alchFepWCADispOn = false;
321     alchFepElecOn    = true;
322     ComputeNonbondedUtil::select();
323     return;
324   }
325
326
327   if ( ! strncasecmp(param,"alchFepWCArcut1",MAX_SCRIPT_PARAM_SIZE) ) {
328     alchFepWCArcut1 = atof(value);
329     ComputeNonbondedUtil::select();
330     return;
331   }
332
333   if ( ! strncasecmp(param,"alchFepWCArcut2",MAX_SCRIPT_PARAM_SIZE) ) {
334     alchFepWCArcut2 = atof(value);
335     ComputeNonbondedUtil::select();
336     return;
337   }
338
339   if ( ! strncasecmp(param,"alchFepWCArcut3",MAX_SCRIPT_PARAM_SIZE) ) {
340     alchFepWCArcut3 = atof(value);
341     ComputeNonbondedUtil::select();
342     return;
343   }
344
345   if ( ! strncasecmp(param,"alchLambda",MAX_SCRIPT_PARAM_SIZE) ) {
346     alchLambda = atof(value);
347     if ( alchLambda < 0.0 || 1.0 < alchLambda ) {
348       NAMD_die("Alchemical lambda values should be in the range [0.0, 1.0]\n");
349     }
350     ComputeNonbondedUtil::select();
351     return;
352   }
353
354   if ( ! strncasecmp(param,"alchLambda2",MAX_SCRIPT_PARAM_SIZE) ) {
355     alchLambda2 = atof(value);
356     if ( alchLambda2 < 0.0 || 1.0 < alchLambda2 ) {
357       NAMD_die("Alchemical lambda values should be in the range [0.0, 1.0]\n");
358     }
359     ComputeNonbondedUtil::select();
360     return;
361   }
362
363   if ( ! strncasecmp(param,"alchLambdaIDWS",MAX_SCRIPT_PARAM_SIZE) ) {
364     alchLambdaIDWS = atof(value);
365     setupIDWS();
366     ComputeNonbondedUtil::select();
367     return;
368   }
369
370   if ( ! strncasecmp(param,"alchLambdaFreq",MAX_SCRIPT_PARAM_SIZE) ) {
371     alchLambdaFreq = atoi(value);
372     if ( alchLambdaIDWS >= 0 ) {
373       NAMD_die("alchLambdaIDWS and alchLambdaFreq are not compatible.\n");
374     }
375     ComputeNonbondedUtil::select();
376     return;
377   }
378 //fepe
379
380   if ( ! strncasecmp(param,"nonbondedScaling",MAX_SCRIPT_PARAM_SIZE) ) {
381     nonbondedScaling = atof(value);
382     ComputeNonbondedUtil::select();
383     return;
384   }
385
386   if ( ! strncasecmp(param,"commOnly",MAX_SCRIPT_PARAM_SIZE) ) {
387     commOnly = atobool(value);
388     ComputeNonbondedUtil::select();
389     return;
390   }
391
392   // REST2 - We don't have to make any changes to ComputeNonbondedUtil
393   // other than recalculating the LJTable.  Skip doing the other stuff.
394   if ( ! strncasecmp(param,"soluteScalingFactor",MAX_SCRIPT_PARAM_SIZE)) {
395     soluteScalingFactor = atof(value);
396     if (soluteScalingFactor < 0.0) {
397       NAMD_die("Solute scaling factor should be non-negative\n");
398     }
399     soluteScalingFactorCharge = soluteScalingFactor;
400     soluteScalingFactorVdw = soluteScalingFactor;
401     // update LJTable for CPU
402     ComputeNonbondedUtil::select();
403 #ifdef NAMD_CUDA
404     // update LJTable for GPU, needs CPU update first
405     ComputeCUDAMgr::getComputeCUDAMgr()->update();
406 #endif
407     return;
408   }
409   if ( ! strncasecmp(param,"soluteScalingFactorVdw",MAX_SCRIPT_PARAM_SIZE)) {
410     soluteScalingFactorVdw = atof(value);
411     if (soluteScalingFactorVdw < 0.0) {
412       NAMD_die("Solute scaling factor for van der Waals "
413           "should be non-negative\n");
414     }
415     // update LJTable for CPU
416     ComputeNonbondedUtil::select();
417 #ifdef NAMD_CUDA
418     // update LJTable for GPU, needs CPU update first
419     ComputeCUDAMgr::getComputeCUDAMgr()->update();
420 #endif
421     return;
422   }
423   if ( ! strncasecmp(param,"soluteScalingFactorCharge",MAX_SCRIPT_PARAM_SIZE)) {
424     soluteScalingFactorCharge = atof(value);
425     if (soluteScalingFactorCharge < 0.0) {
426       NAMD_die("Solute scaling factor for electrostatics "
427           "should be non-negative\n");
428     }
429     return;
430   }
431
432   char *error = new char[2 * MAX_SCRIPT_PARAM_SIZE + 100];
433   sprintf(error,"Setting parameter %s from script failed!\n",param);
434   NAMD_die(error);
435
436 }
437
438 void SimParameters::nonbonded_select() {
439     ComputeNonbondedUtil::select();
440 }
441
442 void SimParameters::pme_select() {
443   ComputePmeUtil::select();
444 }
445
446 /************************************************************************/
447 /*                                                                      */
448 /*      FUNCTION config_parser                                          */
449 /*                                                                      */
450 /************************************************************************/
451                          
452 void SimParameters::config_parser(ParseOptions &opts) {
453
454    //  Set all variable to fallback default values.  This is not really
455    //  necessary, as we give default values when we set up the ParseOptions
456    //  object, but it helps the debuggers figure out we've initialized the
457    //  variables.
458    HydrogenBonds = FALSE;
459    useAntecedent = TRUE;
460    aaAngleExp = 2;
461    haAngleExp = 4;
462    distAttExp = 4;
463    distRepExp = 6;
464    dhaCutoffAngle = 100.0;
465    dhaOnAngle = 60.0;
466    dhaOffAngle = 80.0;
467    daCutoffDist = 7.5;
468    daOnDist = 5.5;
469    daOffDist = 6.5;
470
471    config_parser_basic(opts);
472    config_parser_fileio(opts);
473    config_parser_fullelect(opts);
474    config_parser_methods(opts);
475    config_parser_constraints(opts);
476    #ifdef OPENATOM_VERSION
477    config_parser_openatom(opts);
478    #endif // OPENATOM_VERSION
479
480    config_parser_gridforce(opts);
481    config_parser_mgridforce(opts);
482    config_parser_movdrag(opts);
483    config_parser_rotdrag(opts);
484    config_parser_constorque(opts);
485    config_parser_boundary(opts);
486    config_parser_misc(opts);
487
488 }
489
490 void SimParameters::config_parser_basic(ParseOptions &opts) {
491    
492    //  So first we set up the ParseOptions objects so that it will check
493    //  all of the logical rules that the configuration file must follow.
494
495    opts.optional("main", "obsolete", "used to flag obsolete options",
496     PARSE_STRING);
497
498    ////// basic options
499    opts.require("main", "timestep", "size of the timestep, in fs",
500     &dt, 1.0);
501    opts.range("timestep", NOT_NEGATIVE);
502    opts.units("timestep", N_FSEC);
503
504    opts.optional("main", "numsteps", "number of timesteps to perform",
505     &N,0);
506    opts.range("numsteps", NOT_NEGATIVE);
507
508    opts.optional("main", "stepspercycle",
509       "Number of steps between atom migrations", 
510       &stepsPerCycle, 20);
511    opts.range("stepspercycle", POSITIVE);
512
513    opts.require("main", "cutoff", "local electrostatic and Vdw distance", 
514       &cutoff);
515    opts.range("cutoff", POSITIVE);
516    opts.units("cutoff", N_ANGSTROM);
517    
518    opts.optional("main", "nonbondedScaling", "nonbonded scaling factor",
519      &nonbondedScaling, 1.0);
520    opts.range("nonbondedScaling", NOT_NEGATIVE);
521
522    opts.optional("main", "limitDist", "limit nonbonded below this distance",
523      &limitDist, 0.0);
524    opts.range("limitDist", NOT_NEGATIVE);
525
526    opts.require("main", "exclude", "Electrostatic and VDW exclusion policy",
527     PARSE_STRING);
528
529    opts.optional("exclude", "1-4scaling", "1-4 electrostatic scaling factor",
530      &scale14, 1.0);
531    opts.range("1-4scaling", POSITIVE);
532
533    opts.optionalB("main", "switching",
534      "Should a smoothing function be used?", &switchingActive, TRUE);
535    
536    opts.optionalB("switching", "vdwForceSwitching",
537      "Use force switching for vdw?", &vdwForceSwitching, FALSE);
538    
539    opts.optional("switching", "switchdist",
540      "Distance for switching function activation",
541      &switchingDist);
542    opts.range("switchdist", POSITIVE);
543    opts.units("switchdist", N_ANGSTROM);
544    
545    opts.optionalB("main", "martiniSwitching",
546      "Use Martini residue-based coarse-grain switching?", &martiniSwitching, FALSE);
547    opts.optionalB("main", "martiniDielAllow",
548      "Allow use of dielectric != 15.0 when using Martini", &martiniDielAllow, FALSE);
549
550    opts.optional("main", "pairlistdist",  "Pairlist inclusion distance",
551      &pairlistDist);
552    opts.range("pairlistdist", POSITIVE);
553    opts.units("pairlistdist", N_ANGSTROM);
554
555    opts.optional("main", "pairlistMinProcs",  "Min procs for pairlists",
556      &pairlistMinProcs,1);
557    opts.range("pairlistMinProcs", POSITIVE);
558
559    opts.optional("main", "pairlistsPerCycle",  "regenerate x times per cycle",
560      &pairlistsPerCycle,2);
561    opts.range("pairlistsPerCycle", POSITIVE);
562
563    opts.optional("main", "outputPairlists", "how often to print warnings",
564      &outputPairlists, 0);
565    opts.range("outputPairlists", NOT_NEGATIVE);
566
567    opts.optional("main", "pairlistShrink",  "tol *= (1 - x) on regeneration",
568      &pairlistShrink,0.01);
569    opts.range("pairlistShrink", NOT_NEGATIVE);
570
571    opts.optional("main", "pairlistGrow",  "tol *= (1 + x) on trigger",
572      &pairlistGrow, 0.01);
573    opts.range("pairlistGrow", NOT_NEGATIVE);
574
575    opts.optional("main", "pairlistTrigger",  "trigger is atom > (1 - x) * tol",
576      &pairlistTrigger, 0.3);
577    opts.range("pairlistTrigger", NOT_NEGATIVE);
578
579    opts.optional("main", "temperature", "initial temperature",
580      &initialTemp);
581    opts.range("temperature", NOT_NEGATIVE);
582    opts.units("temperature", N_KELVIN);
583
584    opts.optionalB("main", "COMmotion", "allow initial center of mass movement",
585       &comMove, FALSE);
586
587    opts.optionalB("main", "zeroMomentum", "constrain center of mass",
588       &zeroMomentum, FALSE);
589    opts.optionalB("zeroMomentum", "zeroMomentumAlt", "constrain center of mass",
590       &zeroMomentumAlt, FALSE);
591
592    opts.optionalB("main", "wrapWater", "wrap waters around periodic boundaries on output",
593       &wrapWater, FALSE);
594    opts.optionalB("main", "wrapAll", "wrap all clusters around periodic boundaries on output",
595       &wrapAll, FALSE);
596    opts.optionalB("main", "wrapNearest", "wrap to nearest image to cell origin",
597       &wrapNearest, FALSE);
598
599    opts.optional("main", "dielectric", "dielectric constant",
600      &dielectric, 1.0);
601    opts.range("dielectric", POSITIVE); // Hmmm, dielectric < 1 ...
602
603    opts.optional("main", "margin", "Patch width margin", &margin, XXXBIGREAL);
604    opts.range("margin", NOT_NEGATIVE);
605    opts.units("margin", N_ANGSTROM);
606
607    opts.optional("main", "seed", "Initial random number seed", &randomSeed);
608    opts.range("seed", POSITIVE);
609
610    opts.optional("main", "outputEnergies", "How often to print energies in timesteps",
611      &outputEnergies, 1);
612    opts.range("outputEnergies", POSITIVE);
613      
614    opts.optional("main", "outputMomenta", "How often to print linear and angular momenta in timesteps",
615      &outputMomenta, 0);
616    opts.range("outputMomenta", NOT_NEGATIVE);
617      
618    opts.optional("main", "outputTiming", "How often to print timing data in timesteps",
619      &outputTiming);
620    opts.range("outputTiming", NOT_NEGATIVE);
621      
622    opts.optional("main", "outputCudaTiming", "How often to print CUDA timing data in timesteps",
623      &outputCudaTiming, 0);
624    opts.range("outputCudaTiming", NOT_NEGATIVE);
625      
626    opts.optional("main", "outputPressure", "How often to print pressure data in timesteps",
627      &outputPressure, 0);
628    opts.range("outputPressure", NOT_NEGATIVE);
629      
630    opts.optionalB("main", "mergeCrossterms", "merge crossterm energy with dihedral when printing?",
631       &mergeCrossterms, TRUE);
632
633    opts.optional("main", "MTSAlgorithm", "Multiple timestep algorithm",
634     PARSE_STRING);
635
636    opts.optional("main", "longSplitting", "Long range force splitting option",
637     PARSE_STRING);
638
639    opts.optionalB("main", "ignoreMass", "Do not use masses to find hydrogen atoms",
640     &ignoreMass, FALSE);
641
642    opts.optional("main", "splitPatch", "Atom into patch splitting option",
643     PARSE_STRING);
644    opts.optional("main", "hgroupCutoff", "Hydrogen margin", &hgroupCutoff, 2.5);
645
646    opts.optional("main", "extendedSystem",
647     "Initial configuration of extended system variables and periodic cell",
648     PARSE_STRING);
649
650    opts.optional("main", "cellBasisVector1", "Basis vector for periodic cell",
651     &cellBasisVector1);
652    opts.optional("main", "cellBasisVector2", "Basis vector for periodic cell",
653     &cellBasisVector2);
654    opts.optional("main", "cellBasisVector3", "Basis vector for periodic cell",
655     &cellBasisVector3);
656    opts.optional("main", "cellOrigin", "Fixed center of periodic cell",
657     &cellOrigin);
658
659    opts.optionalB("main", "molly", "Rigid bonds to hydrogen",&mollyOn,FALSE);
660    opts.optional("main", "mollyTolerance", "Error tolerance for MOLLY",
661                  &mollyTol, 0.00001);
662    opts.optional("main", "mollyIterations", 
663                  "Max number of iterations for MOLLY", &mollyIter, 100);
664
665    opts.optional("main", "rigidBonds", "Rigid bonds to hydrogen",PARSE_STRING);
666    opts.optional("main", "rigidTolerance", 
667                  "Error tolerance for rigid bonds to hydrogen",
668                  &rigidTol, 1.0e-8);
669    opts.optional("main", "rigidIterations", 
670                  "Max number of SHAKE iterations for rigid bonds to hydrogen",
671                  &rigidIter, 100);
672    opts.optionalB("main", "rigidDieOnError", 
673                  "Die if rigidTolerance is not achieved after rigidIterations",
674                  &rigidDie, TRUE);
675    opts.optionalB("main", "useSettle",
676                   "Use the SETTLE algorithm for rigid waters",
677                  &useSettle, TRUE);
678
679    opts.optional("main", "nonbondedFreq", "Nonbonded evaluation frequency",
680     &nonbondedFrequency, 1);
681    opts.range("nonbondedFreq", POSITIVE);
682
683    opts.optionalB("main", "outputPatchDetails", "print number of atoms in each patch",
684       &outputPatchDetails, FALSE);
685    opts.optionalB("main", "staticAtomAssignment", "never migrate atoms",
686       &staticAtomAssignment, FALSE);
687    opts.optionalB("main", "replicaUniformPatchGrids", "same patch grid size on all replicas",
688       &replicaUniformPatchGrids, FALSE);
689 #ifndef MEM_OPT_VERSION
690    // in standard (non-mem-opt) version, enable lone pairs by default
691    // for compatibility with recent force fields
692    opts.optionalB("main", "lonePairs", "Enable lone pairs", &lonepairs, TRUE);
693 #else
694    // in mem-opt version, disable lone pairs by default
695    // because they are not supported
696    opts.optionalB("main", "lonePairs", "Enable lone pairs", &lonepairs, FALSE);
697 #endif
698    opts.optional("main", "waterModel", "Water model to use", PARSE_STRING);
699    opts.optionalB("main", "LJcorrection", "Apply analytical tail corrections for energy and virial", &LJcorrection, FALSE);
700 }
701
702 void SimParameters::config_parser_fileio(ParseOptions &opts) {
703    
704    /////////////// file I/O
705
706    opts.optional("main", "cwd", "current working directory", PARSE_STRING);
707
708 // In order to include AMBER options, "coordinates", "structure"
709 // and "parameters" are now optional, not required. The presence
710 // of them will be checked later in check_config()
711
712 //   opts.require("main", "coordinates", "initial PDB coordinate file",
713 //    PARSE_STRING);
714    opts.optional("main", "coordinates", "initial PDB coordinate file",
715     PARSE_STRING);
716
717    opts.optional("main", "velocities",
718      "initial velocities, given as a PDB file", PARSE_STRING);
719    opts.optional("main", "binvelocities",
720      "initial velocities, given as a binary restart", PARSE_STRING);
721    opts.optional("main", "bincoordinates",
722      "initial coordinates in a binary restart file", PARSE_STRING);
723 #ifdef MEM_OPT_VERSION
724    opts.optional("main", "binrefcoords",
725      "reference coordinates in a binary restart file", PARSE_STRING);
726 #endif
727
728 //   opts.require("main", "structure", "initial PSF structure file",
729 //    PARSE_STRING);
730    opts.optional("main", "structure", "initial PSF structure file",
731     PARSE_STRING);   
732
733 //   opts.require("main", "parameters",
734 //"CHARMm 19 or CHARMm 22 compatable force field file (multiple "
735 //"inputs allowed)", PARSE_MULTIPLES);
736    opts.optional("main", "parameters",
737 "CHARMm 19 or CHARMm 22 compatable force field file (multiple "
738 "inputs allowed)", PARSE_MULTIPLES);
739
740
741    //****** BEGIN CHARMM/XPLOR type changes
742    //// enable XPLOR as well as CHARMM input files for parameters
743    opts.optionalB("parameters", "paraTypeXplor", "Parameter file in Xplor format?", &paraTypeXplorOn, FALSE);
744    opts.optionalB("parameters", "paraTypeCharmm", "Parameter file in Charmm format?", &paraTypeCharmmOn, FALSE); 
745    //****** END CHARMM/XPLOR type changes
746
747    // Ported by JLai -- JE - Go parameters
748    opts.optionalB("main", "GromacsPair", "Separately calculate pair interactions", &goGroPair, FALSE);
749    opts.optionalB("main", "GoForcesOn", "Go forces will be calculated", &goForcesOn, FALSE);
750    opts.require("GoForcesOn", "GoParameters", "Go parameter file", goParameters);
751    opts.require("GoForcesOn", "GoCoordinates", "target coordinates for Go forces", goCoordinates);
752    // Added by JLai -- Go-method parameter -- switches between using the matrix and sparse matrix representations [6.3.11] 
753    //   opts.optional("GoForcesOn", "GoMethod", "Which type of matrix should be used to store Go contacts?", &goMethod);
754    // Added by JLai -- Go-method parameter [8.2.11]
755    //opts.optional("GoForcesOn", "GoMethod", "Which type of matrix should be used to store Go contacts?", PARSE_STRING); 
756    opts.require("GoForcesOn", "GoMethod", "Which type of matrix should be used to store Go contacts?", PARSE_STRING); 
757   // End of Port -- JL
758    
759    opts.require("main", "outputname",
760     "prefix for the final PDB position and velocity filenames", 
761     outputFilename);
762
763    opts.optional("main", "auxFile", "Filename for data stream output",
764      auxFilename);
765
766    opts.optional("main", "numinputprocs", "Number of pes to use for parallel input", 
767                  &numinputprocs, 0);
768    opts.range("numinputprocs", NOT_NEGATIVE);
769
770    opts.optional("main", "numoutputprocs", "Number of pes to use for parallel output", 
771                  &numoutputprocs, 0);
772    opts.range("numoutputprocs", NOT_NEGATIVE);
773    opts.optional("main", "numoutputwriters", "Number of output processors that simultaneously write to an output file", 
774                  &numoutputwrts, 1);
775    opts.range("numoutputwriters", NOT_NEGATIVE);
776
777    opts.optional("main", "DCDfreq", "Frequency of DCD trajectory output, in "
778     "timesteps", &dcdFrequency, 0);
779    opts.range("DCDfreq", NOT_NEGATIVE);
780    opts.optional("DCDfreq", "DCDfile", "DCD trajectory output file name",
781      dcdFilename);
782    opts.optionalB("DCDfreq", "DCDunitcell", "Store unit cell in dcd timesteps?",
783        &dcdUnitCell);
784
785    opts.optional("main", "velDCDfreq", "Frequency of velocity "
786     "DCD output, in timesteps", &velDcdFrequency, 0);
787    opts.range("velDCDfreq", NOT_NEGATIVE);
788    opts.optional("velDCDfreq", "velDCDfile", "velocity DCD output file name",
789      velDcdFilename);
790    
791    opts.optional("main", "forceDCDfreq", "Frequency of force"
792     "DCD output, in timesteps", &forceDcdFrequency, 0);
793    opts.range("forceDCDfreq", NOT_NEGATIVE);
794    opts.optional("forceDCDfreq", "forceDCDfile", "force DCD output file name",
795      forceDcdFilename);
796    
797    opts.optional("main", "XSTfreq", "Frequency of XST trajectory output, in "
798     "timesteps", &xstFrequency, 0);
799    opts.range("XSTfreq", NOT_NEGATIVE);
800    opts.optional("XSTfreq", "XSTfile", "Extended sytem trajectory output "
801     "file name", xstFilename);
802
803    opts.optional("main", "restartfreq", "Frequency of restart file "
804     "generation", &restartFrequency, 0);
805    opts.range("restartfreq", NOT_NEGATIVE);
806    opts.optional("restartfreq", "restartname", "Prefix for the position and "
807      "velocity PDB files used for restarting", restartFilename);
808    opts.optionalB("restartfreq", "restartsave", "Save restart files with "
809      "unique filenames rather than overwriting", &restartSave, FALSE);
810    opts.optionalB("restartfreq", "restartsavedcd", "Save DCD files with "
811      "unique filenames at each restart", &restartSaveDcd, FALSE);
812
813    opts.optionalB("restartfreq", "binaryrestart", "Specify use of binary restart files ", 
814        &binaryRestart, TRUE);
815
816    opts.optionalB("outputname", "binaryoutput", "Specify use of binary output files ", 
817        &binaryOutput, TRUE);
818
819    opts.optionalB("main", "amber", "Is it AMBER force field?",
820        &amberOn, FALSE);
821    opts.optionalB("amber", "readexclusions", "Read exclusions from parm file?",
822        &readExclusions, TRUE);
823    opts.require("amber", "scnb", "1-4 VDW interactions are divided by scnb",
824        &vdwscale14, 2.0);
825    opts.require("amber", "parmfile", "AMBER parm file", PARSE_STRING);
826    opts.optional("amber", "ambercoor", "AMBER coordinate file", PARSE_STRING);
827
828    /* GROMACS options */
829    opts.optionalB("main", "gromacs", "Use GROMACS-like force field?",
830        &gromacsOn, FALSE);
831    opts.require("gromacs", "grotopfile", "GROMACS topology file",
832                 PARSE_STRING);
833    opts.optional("gromacs", "grocoorfile","GROMACS coordinate file",
834                  PARSE_STRING);
835
836   // OPLS options
837    opts.optionalB("main", "vdwGeometricSigma",
838        "Use geometric mean to combine L-J sigmas, as for OPLS",
839        &vdwGeometricSigma, FALSE);
840
841    // load/store computeMap
842    opts.optional("main", "computeMapFile", "Filename for computeMap",
843      computeMapFilename);
844    opts.optionalB("main", "storeComputeMap", "store computeMap?",
845        &storeComputeMap, FALSE);
846    opts.optionalB("main", "loadComputeMap", "load computeMap?",
847        &loadComputeMap, FALSE);
848 }
849
850
851 void SimParameters::config_parser_fullelect(ParseOptions &opts) {
852    
853    /////////// FMA options
854 #ifdef DPMTA
855    DebugM(1,"DPMTA setup start\n");
856    //  PMTA is included, so really get these values
857    opts.optionalB("main", "FMA", "Should FMA be used?", &FMAOn, FALSE);
858    opts.optional("FMA", "FMALevels", "Tree levels to use in FMA", &FMALevels,
859      5);
860    opts.range("FMALevels", POSITIVE);
861    opts.optional("FMA", "FMAMp", "Number of FMA multipoles", &FMAMp, 8);
862    opts.range("FMAMp", POSITIVE);
863    opts.optionalB("FMA", "FMAFFT", "Use FFT enhancement in FMA?", &FMAFFTOn, TRUE);
864    opts.optional("FMAFFT", "FMAFFTBlock", "FFT blocking factor",
865     &FMAFFTBlock, 4);
866    opts.range("FMAFFTBlock", POSITIVE);
867    DebugM(1,"DPMTA setup end\n");
868 #else
869    //  PMTA is NOT included.  So just set all the values to 0.
870    FMAOn = FALSE;
871    FMALevels = 0;
872    FMAMp = 0;
873    FMAFFTOn = FALSE;
874    FMAFFTBlock = 0;
875 #endif
876
877    opts.optional("main", "fullElectFrequency",
878       "Number of steps between full electrostatic executions", 
879       &fullElectFrequency);
880    opts.range("fullElectFrequency", POSITIVE);
881
882    //  USE OF THIS PARAMETER DISCOURAGED
883    opts.optional("main", "fmaFrequency",
884       "Number of steps between full electrostatic executions", 
885       &fmaFrequency);
886    opts.range("fmaFrequency", POSITIVE);
887
888    opts.optional("main", "fmaTheta",
889       "FMA theta parameter value", 
890       &fmaTheta,0.715);
891    opts.range("fmaTheta", POSITIVE);
892
893    opts.optionalB("main", "FullDirect", "Should direct calculations of full electrostatics be performed?",
894       &fullDirectOn, FALSE);
895
896
897    ///////////  Multilevel Summation Method
898
899    opts.optionalB("main", "MSM",
900        "Use multilevel summation method for electrostatics?",
901        &MSMOn, FALSE);
902    opts.optional("MSM", "MSMQuality", "MSM quality",
903        &MSMQuality, 0);
904    opts.optional("MSM", "MSMApprox", "MSM approximation",
905        &MSMApprox, 0);
906    opts.optional("MSM", "MSMSplit", "MSM splitting",
907        &MSMSplit, 0);
908    opts.optional("MSM", "MSMLevels", "MSM maximum number of levels",
909        &MSMLevels, 0);  // set to 0 adapts to as many as needed
910    opts.optional("MSM", "MSMGridSpacing", "MSM grid spacing (Angstroms)",
911        &MSMGridSpacing, 2.5);
912    opts.optional("MSM", "MSMPadding", "MSM padding (Angstroms)",
913        &MSMPadding, 2.5);
914    opts.optional("MSM", "MSMxmin", "MSM x minimum (Angstroms)", &MSMxmin, 0);
915    opts.optional("MSM", "MSMxmax", "MSM x maximum (Angstroms)", &MSMxmax, 0);
916    opts.optional("MSM", "MSMymin", "MSM y minimum (Angstroms)", &MSMymin, 0);
917    opts.optional("MSM", "MSMymax", "MSM y maximum (Angstroms)", &MSMymax, 0);
918    opts.optional("MSM", "MSMzmin", "MSM z minimum (Angstroms)", &MSMzmin, 0);
919    opts.optional("MSM", "MSMzmax", "MSM z maximum (Angstroms)", &MSMzmax, 0);
920    opts.optional("MSM", "MSMBlockSizeX",
921        "MSM grid block size along X direction (for decomposing parallel work)",
922        &MSMBlockSizeX, 8);
923    opts.optional("MSM", "MSMBlockSizeY",
924        "MSM grid block size along Y direction (for decomposing parallel work)",
925        &MSMBlockSizeY, 8);
926    opts.optional("MSM", "MSMBlockSizeZ",
927        "MSM grid block size along Z direction (for decomposing parallel work)",
928        &MSMBlockSizeZ, 8);
929
930    opts.optionalB("MSM", "MsmSerial",
931        "Use MSM serial version for long-range calculation?",
932        &MsmSerialOn, FALSE);
933
934
935    ///////////  Fast Multipole Method
936
937    opts.optionalB("main", "FMM",
938        "Use fast multipole method for electrostatics?",
939        &FMMOn, FALSE);
940    opts.optional("FMM", "FMMLevels", "FMM number of levels",
941        &FMMLevels, 0);
942    opts.optional("FMM", "FMMPadding", "FMM padding margin (Angstroms)",
943        &FMMPadding, 0);
944
945    opts.optionalB("main", "useCUDA2", "Use new CUDA code", &useCUDA2, TRUE);
946
947    ///////////  Particle Mesh Ewald
948
949    opts.optionalB("main", "PME", "Use particle mesh Ewald for electrostatics?",
950         &PMEOn, FALSE);
951    opts.optional("PME", "PMETolerance", "PME direct space tolerance",
952         &PMETolerance, 1.e-6);
953    opts.optional("PME", "PMEInterpOrder", "PME interpolation order",
954         &PMEInterpOrder, 4);  // cubic interpolation is default
955    opts.optional("PME", "PMEGridSizeX", "PME grid in x dimension",
956         &PMEGridSizeX, 0);
957    opts.optional("PME", "PMEGridSizeY", "PME grid in y dimension",
958         &PMEGridSizeY, 0);
959    opts.optional("PME", "PMEGridSizeZ", "PME grid in z dimension",
960         &PMEGridSizeZ, 0);
961    opts.optional("PME", "PMEGridSpacing", "Maximum PME grid spacing (Angstroms)",
962         &PMEGridSpacing, 0.);
963    opts.range("PMEGridSpacing", NOT_NEGATIVE);
964    opts.optional("PME", "PMEProcessors",
965         "PME FFT and reciprocal sum processor count", &PMEProcessors, 0);
966    opts.optional("PME", "PMEMinSlices",
967         "minimum thickness of PME reciprocal sum slab", &PMEMinSlices, 2);
968    opts.range("PMEMinSlices", NOT_NEGATIVE);
969    opts.optional("PME", "PMEPencils",
970         "PME FFT and reciprocal sum pencil grid size", &PMEPencils, -1);
971    opts.optional("PME", "PMEPencilsX",
972         "PME FFT and reciprocal sum pencil grid size X", &PMEPencilsX, 0);
973    opts.optional("PME", "PMEPencilsY",
974         "PME FFT and reciprocal sum pencil grid size Y", &PMEPencilsY, 0);
975    opts.optional("PME", "PMEPencilsZ",
976         "PME FFT and reciprocal sum pencil grid size Z", &PMEPencilsZ, 0);
977    opts.range("PMEPencilsX", NOT_NEGATIVE);
978    opts.range("PMEPencilsY", NOT_NEGATIVE);
979    opts.range("PMEPencilsZ", NOT_NEGATIVE);
980    opts.optional("PME", "PMEPencilsYLayout",
981         "PME FFT and reciprocal sum Y pencil layout strategy", &PMEPencilsYLayout, 0);
982    opts.optional("PME", "PMEPencilsXLayout",
983         "PME FFT and reciprocal sum X pencil layout strategy", &PMEPencilsXLayout, 1);
984    opts.range("PMEPencilsYLayout", NOT_NEGATIVE);
985    opts.range("PMEPencilsXLayout", NOT_NEGATIVE);
986    opts.optional("PME", "PMESendOrder",
987         "PME message ordering control", &PMESendOrder, 0);
988    opts.range("PMESendOrder", NOT_NEGATIVE);
989    opts.optional("PME", "PMEMinPoints",
990         "minimum points per PME reciprocal sum pencil", &PMEMinPoints, 10000);
991    opts.range("PMEMinPoints", NOT_NEGATIVE);
992    opts.optionalB("main", "PMEBarrier", "Use barrier in PME?",
993         &PMEBarrier, FALSE);
994    opts.optionalB("main", "PMEOffload", "Offload PME to accelerator?",
995         &PMEOffload);
996
997    opts.optionalB("PME", "usePMECUDA", "Use the PME CUDA version", &usePMECUDA, CmiNumPhysicalNodes() < 5);
998
999 #ifdef DPME
1000    opts.optionalB("PME", "useDPME", "Use old DPME code?", &useDPME, FALSE);
1001 #else
1002    useDPME = 0;
1003 #endif
1004    opts.optionalB("main", "FFTWPatient", "Use intensive plan creation to optimize FFTW?",
1005 #ifdef WIN32
1006         &FFTWPatient, TRUE);
1007 #else
1008         &FFTWPatient, FALSE);
1009 #endif
1010
1011    opts.optionalB("main", "FFTWEstimate", "Use estimates to optimize FFTW?",
1012 #ifdef WIN32
1013         &FFTWEstimate, TRUE);
1014 #else
1015         &FFTWEstimate, FALSE);
1016 #endif
1017    opts.optionalB("main", "FFTWUseWisdom", "Read/save wisdom file for FFTW?",
1018 #ifdef WIN32
1019         &FFTWUseWisdom, FALSE);
1020 #else
1021         &FFTWUseWisdom, TRUE);
1022 #endif
1023    opts.optional("FFTWUseWisdom", "FFTWWisdomFile", "File for FFTW wisdom",
1024         FFTWWisdomFile);
1025
1026 }
1027
1028 void SimParameters::config_parser_methods(ParseOptions &opts) {
1029    
1030    /////////// Special Dynamics Methods
1031    opts.optionalB("main", "minimization", "Should minimization be performed?",
1032       &minimizeCGOn, FALSE);
1033    opts.optionalB("main", "minVerbose", "Print extra minimization diagnostics?",
1034       &minVerbose, FALSE);
1035    opts.optional("main", "minTinyStep", "very first minimization steps",
1036       &minTinyStep, 1.0e-6);
1037    opts.range("minTinyStep", POSITIVE);
1038    opts.optional("main", "minBabyStep", "initial minimization steps",
1039       &minBabyStep, 1.0e-2);
1040    opts.range("minBabyStep", POSITIVE);
1041    opts.optional("main", "minLineGoal", "line minimization gradient reduction",
1042       &minLineGoal, 1.0e-3);
1043    opts.range("minLineGoal", POSITIVE);
1044
1045    opts.optionalB("main", "velocityQuenching",
1046       "Should old-style minimization be performed?", &minimizeOn, FALSE);
1047
1048    opts.optional("main", "maximumMove", "Maximum atom movement per step", &maximumMove, 0.0);
1049    opts.range("maximumMove", NOT_NEGATIVE);
1050    opts.units("maximumMove", N_ANGSTROM);
1051
1052    opts.optionalB("main", "Langevin", "Should Langevin dynamics be performed?",
1053       &langevinOn, FALSE);
1054    opts.require("Langevin", "langevinTemp", "Temperature for heat bath in Langevin "
1055      "dynamics", &langevinTemp);
1056    opts.range("langevinTemp", NOT_NEGATIVE);
1057    opts.units("langevinTemp", N_KELVIN);
1058    opts.optional("Langevin", "langevinDamping", "Damping coefficient (1/ps)",
1059       &langevinDamping);
1060    opts.range("langevinDamping", POSITIVE);
1061    opts.optionalB("Langevin", "langevinHydrogen", "Should Langevin dynamics be applied to hydrogen atoms?",
1062       &langevinHydrogen);
1063    opts.optional("Langevin", "langevinFile", "PDB file with temperature "
1064      "coupling terms (B(i)) (default is the PDB input file)",
1065      PARSE_STRING);
1066    opts.optional("Langevin", "langevinCol", "Column in the langevinFile "
1067      "containing the temperature coupling term B(i);\n"
1068      "default is 'O'", PARSE_STRING);
1069
1070    // use BAOAB integration instead of BBK
1071    opts.optionalB("Langevin", "langevinBAOAB",
1072        "Should Langevin dynamics be performed using BAOAB integration?",
1073        &langevin_useBAOAB, FALSE);
1074
1075 // BEGIN LA
1076    opts.optionalB("main", "LoweAndersen", "Should Lowe-Andersen dynamics be performed?",
1077                   &loweAndersenOn, FALSE);
1078    opts.require("LoweAndersen", "loweAndersenTemp", "Temperature for heat bath in Lowe-Andersen "
1079                 "dynamics", &loweAndersenTemp);
1080    opts.range("loweAndersenTemp", NOT_NEGATIVE);
1081    opts.units("loweAndersenTemp", N_KELVIN);
1082    opts.optional("LoweAndersen", "loweAndersenRate", "Collision rate (1/ps)",
1083                  &loweAndersenRate, 50);
1084    opts.range("loweAndersenRate", POSITIVE);
1085    opts.optional("LoweAndersen", "loweAndersenCutoff", "Cutoff radius",
1086                  &loweAndersenCutoff, 2.7);
1087    opts.range("loweAndersenCutoff", POSITIVE);
1088    opts.units("loweAndersenCutoff", N_ANGSTROM);
1089 // END LA
1090
1091 //fepb
1092    opts.optionalB("main", "alch", "Is achemical simulation being performed?",
1093      &alchOn, FALSE);
1094    opts.require("alch", "alchLambda", "Alchemical coupling parameter value",
1095      &alchLambda);
1096    opts.range("alchLambda", NOT_NEGATIVE);
1097
1098    opts.optional("alch", "alchFile", "PDB file with perturbation flags "
1099      "default is the input PDB file", PARSE_STRING);
1100    opts.optional("alch", "alchCol", "Column in the alchFile with the "
1101      "perturbation flag", PARSE_STRING);
1102
1103    opts.optional("alch", "alchOutFreq", "Frequency of alchemical energy"
1104      "output in timesteps", &alchOutFreq, 5);
1105    opts.range("alchoutfreq", NOT_NEGATIVE);
1106    opts.optional("alch", "alchOutFile", "Alchemical energy output filename",
1107      alchOutFile);
1108
1109    // soft-core parameters
1110    opts.optional("alch", "alchVdwShiftCoeff", "Coeff used for generating"
1111      "the altered alchemical vDW interactions", &alchVdwShiftCoeff, 5.);
1112    opts.range("alchVdwShiftCoeff", NOT_NEGATIVE);
1113
1114    // scheduling options for different interaction types
1115    opts.optional("alch", "alchElecLambdaStart", "Lambda at which electrostatic"
1116       "scaling of exnihilated particles begins", &alchElecLambdaStart, 0.5);
1117    opts.range("alchElecLambdaStart", NOT_NEGATIVE);
1118
1119    opts.optional("alch", "alchVdwLambdaEnd", "Lambda at which vdW"
1120       "scaling of exnihilated particles begins", &alchVdwLambdaEnd, 1.0);
1121    opts.range("alchVdwLambdaEnd", NOT_NEGATIVE);
1122
1123    opts.optional("alch", "alchBondLambdaEnd", "Lambda at which bonded"
1124       "scaling of exnihilated particles begins", &alchBondLambdaEnd, 0.0);
1125    opts.range("alchBondLambdaEnd", NOT_NEGATIVE);
1126    
1127    opts.optionalB("alch", "alchDecouple", "Enable alchemical decoupling?",
1128      &alchDecouple, FALSE);
1129    opts.optionalB("alch", "alchBondDecouple", "Enable decoupling of purely "
1130      "alchemical bonds?", &alchBondDecouple, FALSE);
1131
1132    // parameters for alchemical analysis options
1133    opts.optional("alch", "alchType", "Which alchemical method to use?",
1134      PARSE_STRING);
1135    opts.optional("alch", "alchLambda2", "Alchemical coupling comparison value",
1136      &alchLambda2, -1);
1137    opts.optional("alch", "alchLambdaIDWS", "Alchemical coupling comparison value for interleaved double-wide sampling",
1138      &alchLambdaIDWS, -1);
1139    opts.optional("alch", "alchLambdaFreq",
1140      "Frequency of increasing coupling parameter value", &alchLambdaFreq, 0);
1141    opts.range("alchLambdaFreq", NOT_NEGATIVE);
1142    opts.optional("alch", "alchSwitchType", "Switching type flag",
1143      PARSE_STRING);
1144    opts.optional("alch", "alchEquilSteps", "Equilibration steps, before "
1145      "data collection in the alchemical window", &alchEquilSteps, 0);
1146    opts.range("alchEquilSteps", NOT_NEGATIVE);
1147
1148    // WCA decomposition options
1149    opts.optionalB("alch", "alchFepWCARepuOn",
1150      "WCA decomposition repu interaction in use?", &alchFepWCARepuOn, FALSE);
1151    opts.optionalB("alch", "alchFepWCADispOn",
1152      "WCA decomposition disp interaction in use?", &alchFepWCADispOn, FALSE);
1153    opts.optionalB("alch", "alchEnsembleAvg", "Ensemble Average in use?",
1154      &alchEnsembleAvg, TRUE);
1155    opts.optionalB("alch", "alchFepWhamOn",
1156      "Energy output for Wham postprocessing in use?", &alchFepWhamOn, FALSE);
1157    opts.optional("alch", "alchFepWCArcut1",
1158      "WCA repulsion Coeff1 used for generating the altered alchemical vDW "
1159      "interactions", &alchFepWCArcut1, 0.0);
1160    opts.range("alchFepWCArcut1", NOT_NEGATIVE);
1161    opts.optional("alch", "alchFepWCArcut2", "WCA repulsion Coeff2 used for "
1162      "generating the altered alchemical vDW interactions", &alchFepWCArcut2,
1163      1.0);
1164    opts.range("alchFepWCArcut2", NOT_NEGATIVE);
1165    opts.optional("alch", "alchFepWCArcut3",
1166      "WCA repulsion Coeff3 used for generating the altered alchemical vDW "
1167      "interactions", &alchFepWCArcut3, 1.0);
1168    opts.range("alchFepWCArcut3", NOT_NEGATIVE);
1169    // These default to invalid lambda values.
1170    opts.optional("alch", "alchRepLambda", "Lambda of WCA repulsion"
1171      "Coupling parameter value for WCA repulsion", &alchRepLambda, -1.0);
1172    opts.optional("alch", "alchDispLambda", "Lambda of WCA dispersion"
1173      "Coupling parameter value for WCA dispersion", &alchDispLambda, -1.0);
1174    opts.optional("alch", "alchElecLambda", "Lambda of electrostatic "
1175      "perturbation Coupling parameter value for electrostatic perturbation",
1176      &alchElecLambda, -1.0);
1177 //fepe
1178
1179    opts.optionalB("main", "les", "Is locally enhanced sampling enabled?",
1180      &lesOn, FALSE);
1181    opts.require("les", "lesFactor", "Local enhancement factor", &lesFactor);
1182    opts.optional("les", "lesFile", "PDB file with enhancement flags "
1183      "default is the input PDB file", PARSE_STRING); 
1184    opts.optional("les", "lesCol", "Column in the lesFile with the "
1185      "enhancement flag", PARSE_STRING);
1186    opts.optionalB("les", "lesReduceTemp", "Reduce enhanced atom temperature?",
1187      &lesReduceTemp, FALSE);
1188    opts.optionalB("les", "lesReduceMass", "Reduce enhanced atom mass?",
1189      &lesReduceMass, FALSE);
1190
1191    // REST2 (replica exchange solute tempering) parameters
1192    opts.optionalB("main", "soluteScaling",
1193        "Is replica exchange solute tempering enabled?",
1194        &soluteScalingOn, FALSE);
1195    opts.require("soluteScaling", "soluteScalingFactor",
1196        "Solute scaling factor",
1197        &soluteScalingFactor);
1198    opts.range("soluteScalingFactor", NOT_NEGATIVE);
1199    opts.optional("soluteScaling", "soluteScalingFactorCharge",
1200        "Solute scaling factor for electrostatic interactions",
1201        &soluteScalingFactorCharge);
1202    opts.range("soluteScalingFactorCharge", NOT_NEGATIVE);
1203    opts.optional("soluteScaling", "soluteScalingFactorVdw",
1204        "Solute scaling factor for van der Waals interactions",
1205        &soluteScalingFactorVdw);
1206    opts.range("soluteScalingFactorVdw", NOT_NEGATIVE);
1207    opts.optional("soluteScaling", "soluteScalingFile",
1208        "PDB file with scaling flags; if undefined, defaults to main PDB file",
1209        PARSE_STRING);
1210    opts.optional("soluteScaling", "soluteScalingCol",
1211        "Column in the soluteScalingFile providing the scaling flag",
1212        PARSE_STRING);
1213    opts.optionalB("main", "soluteScalingAll",
1214        "Apply scaling also to bond and angle interactions?",
1215        &soluteScalingAll, FALSE);
1216
1217    // Drude oscillators
1218    opts.optionalB("main", "drude", "Perform integration of Drude oscillators?",
1219        &drudeOn, FALSE);
1220    opts.require("drude", "drudeTemp", "Temperature for freezing "
1221        "Drude oscillators", &drudeTemp);
1222    opts.range("drudeTemp", NOT_NEGATIVE);
1223    opts.units("drudeTemp", N_KELVIN);
1224    opts.optional("drude", "drudeDamping", "Damping coefficient (1/ps) for "
1225        "Drude oscillators", &drudeDamping);
1226    opts.range("drudeDamping", POSITIVE);
1227    opts.optional("drude", "drudeNbtholeCut", "Nonbonded Thole interactions "
1228        "interaction radius", &drudeNbtholeCut, 5.0);
1229    opts.range("drudeNbtholeCut", POSITIVE);
1230    opts.optionalB("drude", "drudeHardWall", "Apply maximum Drude bond length "
1231        "restriction?", &drudeHardWallOn, TRUE);
1232    opts.optional("drude", "drudeBondLen", "Drude oscillator bond length "
1233        "beyond which to apply restraint", &drudeBondLen, 0.25);
1234    opts.range("drudeBondLen", POSITIVE);
1235    opts.optional("drude", "drudeBondConst", "Drude oscillator restraining "
1236        "force constant", &drudeBondConst, 40000.0);
1237    opts.range("drudeBondConst", POSITIVE);
1238
1239    // Pair interaction calculations
1240     opts.optionalB("main", "pairInteraction", 
1241         "Are pair interactions calculated?", &pairInteractionOn, FALSE);
1242     opts.optional("pairInteraction", "pairInteractionFile", 
1243         "PDB files with interaction flags " "default is the input PDB file", 
1244         PARSE_STRING);
1245     opts.optional("pairInteraction", "pairInteractionCol", 
1246         "Column in the pairInteractionFile with the interaction flags",
1247         PARSE_STRING);
1248     opts.require("pairInteraction", "pairInteractionGroup1",
1249         "Flag for interaction group 1", &pairInteractionGroup1);
1250     opts.optional("pairInteraction", "pairInteractionGroup2",
1251         "Flag for interaction group 2", &pairInteractionGroup2, -1);
1252     opts.optionalB("pairInteraction", "pairInteractionSelf",
1253         "Compute only within-group interactions?", &pairInteractionSelf, 
1254         FALSE);
1255    // Options for CG simulations
1256    opts.optionalB("main", "cosAngles", "Are some angles cosine-based?", &cosAngles, FALSE);
1257
1258
1259    //  Dihedral angle dynamics
1260    opts.optionalB("main", "globalTest", "Should global integration (for development) be used?",
1261     &globalOn, FALSE);
1262    opts.optionalB("main", "dihedral", "Should dihedral angle dynamics be performed?",
1263     &dihedralOn, FALSE);
1264    COLDOn = FALSE;
1265    opts.optionalB("dihedral", "COLD", "Should overdamped Langevin dynamics be performed?",
1266     &COLDOn, FALSE);
1267    opts.require("COLD", "COLDTemp", "Temperature for heat bath in COLD",
1268     &COLDTemp);
1269    opts.range("COLDTemp", NOT_NEGATIVE);
1270    opts.units("COLDTemp", N_KELVIN);
1271    opts.require("COLD", "COLDRate", "Damping rate for COLD",
1272     &COLDRate, 3000.0);
1273    opts.range("COLDRate", NOT_NEGATIVE);
1274
1275    //  Get the parameters for temperature coupling
1276    opts.optionalB("main", "tcouple", 
1277       "Should temperature coupling be performed?",
1278       &tCoupleOn, FALSE);
1279    opts.require("tcouple", "tCoupleTemp", 
1280     "Temperature for temperature coupling", &tCoupleTemp);
1281    opts.range("tCoupleTemp", NOT_NEGATIVE);
1282    opts.units("tCoupleTemp", N_KELVIN);
1283    opts.optional("tCouple", "tCoupleFile", "PDB file with temperature "
1284      "coupling terms (B(i)) (default is the PDB input file)",
1285      PARSE_STRING);
1286    opts.optional("tCouple", "tCoupleCol", "Column in the tCoupleFile "
1287      "containing the temperature coupling term B(i);\n"
1288      "default is 'O'", PARSE_STRING);
1289
1290    opts.optionalB("main", "stochRescale",
1291       "Should stochastic velocity rescaling be performed?",
1292       &stochRescaleOn, FALSE);
1293    opts.require("stochRescale", "stochRescaleTemp",
1294       "Temperature for stochastic velocity rescaling",
1295        &stochRescaleTemp);
1296    opts.range("stochRescaleTemp", NOT_NEGATIVE);
1297    opts.units("stochRescaleTemp", N_KELVIN);
1298    opts.require("stochRescale", "stochRescalePeriod",
1299       "Time scale for stochastic velocity rescaling (ps)",
1300        &stochRescalePeriod);
1301    opts.range("stochRescalePeriod", POSITIVE);
1302    opts.optional("stochRescale", "stochRescaleFreq",
1303        "Number of steps between stochastic rescalings",
1304         &stochRescaleFreq);
1305    opts.range("stochRescaleFreq", POSITIVE);
1306    opts.optionalB("stochRescale", "stochRescaleHeat",
1307        "Should heat transfer and work be computed?", &stochRescaleHeat, FALSE);
1308
1309    opts.optional("main", "rescaleFreq", "Number of steps between "
1310     "velocity rescaling", &rescaleFreq);
1311    opts.range("rescaleFreq", POSITIVE);
1312    opts.optional("main", "rescaleTemp", "Target temperature for velocity rescaling",
1313     &rescaleTemp);
1314    opts.range("rescaleTemp", NOT_NEGATIVE);
1315    opts.units("rescaleTemp", N_KELVIN);
1316
1317    opts.optional("main", "reassignFreq", "Number of steps between "
1318     "velocity reassignment", &reassignFreq);
1319    opts.range("reassignFreq", POSITIVE);
1320    opts.optional("main", "reassignTemp", "Target temperature for velocity reassignment",
1321     &reassignTemp);
1322    opts.range("reassignTemp", NOT_NEGATIVE);
1323    opts.units("reassignTemp", N_KELVIN);
1324    opts.optional("main", "reassignIncr", "Temperature increment for velocity reassignment",
1325     &reassignIncr);
1326    opts.units("reassignIncr", N_KELVIN);
1327    opts.optional("main", "reassignHold", "Final holding temperature for velocity reassignment",
1328     &reassignHold);
1329    opts.range("reassignHold", NOT_NEGATIVE);
1330    opts.units("reassignHold", N_KELVIN);
1331
1332    ////  Group rather than atomic pressure
1333    opts.optionalB("main", "useGroupPressure", 
1334       "Use group rather than atomic quantities for pressure control?",
1335       &useGroupPressure, FALSE);
1336
1337    ////  Anisotropic cell fluctuations
1338    opts.optionalB("main", "useFlexibleCell",
1339       "Use anisotropic cell fluctuation for pressure control?",
1340       &useFlexibleCell, FALSE);
1341
1342    // Fix specific cell dimensions
1343    opts.optionalB("main", "fixCellDims",
1344       "Fix some cell dimensions?",
1345       &fixCellDims, FALSE);
1346
1347    opts.optionalB("fixCellDims", "fixCellDimX",
1348       "Fix the X dimension?",
1349       &fixCellDimX, FALSE);
1350    opts.optionalB("fixCellDims", "fixCellDimY",
1351       "Fix the Y dimension?",
1352       &fixCellDimY, FALSE);
1353    opts.optionalB("fixCellDims", "fixCellDimZ",
1354       "Fix the Z dimension?",
1355       &fixCellDimZ, FALSE);
1356
1357    ////  Constant dimension ratio in X-Y plane
1358    opts.optionalB("main", "useConstantRatio",
1359       "Use constant X-Y ratio for pressure control?",
1360       &useConstantRatio, FALSE);
1361
1362    ////  Constant area and normal pressure conditions
1363    opts.optionalB("main", "useConstantArea",
1364       "Use constant area for pressure control?",
1365       &useConstantArea, FALSE);
1366
1367    //// Exclude atoms from pressure
1368    opts.optionalB("main", "excludeFromPressure",
1369         "Should some atoms be excluded from pressure rescaling?",
1370         &excludeFromPressure, FALSE);
1371    opts.optional("excludeFromPressure", "excludeFromPressureFile",
1372         "PDB file for atoms to be excluded from pressure",
1373         PARSE_STRING);
1374    opts.optional("excludeFromPressure", "excludeFromPressureCol", 
1375         "Column in the excludeFromPressureFile"
1376         "containing the flags (nonzero means excluded);\n"
1377         "default is 'O'", PARSE_STRING);
1378
1379    ////  Berendsen pressure bath coupling
1380    opts.optionalB("main", "BerendsenPressure", 
1381       "Should Berendsen pressure bath coupling be performed?",
1382       &berendsenPressureOn, FALSE);
1383    opts.require("BerendsenPressure", "BerendsenPressureTarget",
1384     "Target pressure for pressure coupling",
1385     &berendsenPressureTarget);
1386    // opts.units("BerendsenPressureTarget",);
1387    opts.require("BerendsenPressure", "BerendsenPressureCompressibility",
1388     "Isothermal compressibility for pressure coupling",
1389     &berendsenPressureCompressibility);
1390    // opts.units("BerendsenPressureCompressibility",);
1391    opts.require("BerendsenPressure", "BerendsenPressureRelaxationTime",
1392     "Relaxation time for pressure coupling",
1393     &berendsenPressureRelaxationTime);
1394    opts.range("BerendsenPressureRelaxationTime", POSITIVE);
1395    opts.units("BerendsenPressureRelaxationTime", N_FSEC);
1396    opts.optional("BerendsenPressure", "BerendsenPressureFreq",
1397     "Number of steps between volume rescaling",
1398     &berendsenPressureFreq, 1);
1399    opts.range("BerendsenPressureFreq", POSITIVE);
1400
1401    ////  Langevin Piston pressure control
1402    opts.optionalB("main", "LangevinPiston",
1403       "Should Langevin piston pressure control be used?",
1404       &langevinPistonOn, FALSE);
1405    opts.optionalB("LangevinPiston", "LangevinPistonBarrier",
1406       "Should Langevin piston barrier be used?",
1407       &langevinPistonBarrier, TRUE);
1408    opts.require("LangevinPiston", "LangevinPistonTarget",
1409       "Target pressure for pressure control",
1410       &langevinPistonTarget);
1411    opts.require("LangevinPiston", "LangevinPistonPeriod",
1412       "Oscillation period for pressure control",
1413       &langevinPistonPeriod);
1414    opts.range("LangevinPistonPeriod", POSITIVE);
1415    opts.units("LangevinPistonPeriod", N_FSEC);
1416    opts.require("LangevinPiston", "LangevinPistonDecay",
1417       "Decay time for pressure control",
1418       &langevinPistonDecay);
1419    opts.range("LangevinPistonDecay", POSITIVE);
1420    opts.units("LangevinPistonDecay", N_FSEC);
1421    opts.require("LangevinPiston", "LangevinPistonTemp",
1422       "Temperature for pressure control piston",
1423       &langevinPistonTemp);
1424    opts.range("LangevinPistonTemp", POSITIVE);
1425    opts.units("LangevinPistonTemp", N_KELVIN);
1426    opts.optional("LangevinPiston", "StrainRate",
1427       "Initial strain rate for pressure control (x y z)",
1428       &strainRate);
1429
1430    // Multigrator temperature and/or pressure control
1431    opts.optionalB("main", "Multigrator", 
1432       "Should multigrator temperature and/or pressure control be used?",
1433       &multigratorOn, FALSE);
1434    opts.require("Multigrator", "MultigratorPressureTarget",
1435     "Target pressure for pressure coupling",
1436     &multigratorPressureTarget);
1437    opts.require("Multigrator", "MultigratorTemperatureTarget",
1438     "Target temperature for temperature coupling",
1439     &multigratorTemperatureTarget);
1440    opts.require("Multigrator", "MultigratorPressureFreq",
1441     "Number of steps between pressure control moves",
1442     &multigratorPressureFreq);
1443    opts.range("MultigratorPressureFreq", POSITIVE);
1444    opts.optional("Multigrator", "MultigratorPressureRelaxationTime",
1445     "Relaxation time for pressure coupling is fs",
1446     &multigratorPressureRelaxationTime, 30000);
1447    opts.range("MultigratorPressureRelaxationTime", POSITIVE);
1448    opts.units("MultigratorPressureRelaxationTime", N_FSEC);
1449    opts.optional("Multigrator", "MultigratorTemperatureRelaxationTime",
1450     "Relaxation time for temperature coupling is fs",
1451     &multigratorTemperatureRelaxationTime, 1000);
1452    opts.range("MultigratorTemperatureRelaxationTime", POSITIVE);
1453    opts.units("MultigratorTemperatureRelaxationTime", N_FSEC);
1454    opts.require("Multigrator", "MultigratorTemperatureFreq",
1455     "Number of steps between temperature control moves",
1456     &multigratorTemperatureFreq);
1457    opts.range("MultigratorTemperatureFreq", POSITIVE);
1458    opts.optional("Multigrator", "MultigratorNoseHooverChainLength",
1459     "Nose-Hoover chain length",
1460     &multigratorNoseHooverChainLength, 4);
1461    opts.range("MultigratorNoseHooverChainLength", POSITIVE);
1462
1463    //// Surface tension 
1464    opts.optional("main", "SurfaceTensionTarget",
1465       "Surface tension in the x-y plane",
1466       &surfaceTensionTarget, 0);
1467
1468    //// Pressure Profile calculations
1469    opts.optionalB("main", "pressureprofile", "Compute pressure profile?",
1470      &pressureProfileOn, FALSE);
1471    opts.require("pressureprofile", "pressureprofileslabs", 
1472      "Number of pressure profile slabs", &pressureProfileSlabs, 10);
1473    opts.optional("pressureprofile", "pressureprofilefreq",
1474      "How often to store profile data", &pressureProfileFreq, 1);
1475    opts.optional("pressureprofile", "pressureProfileAtomTypes", 
1476      "Number of pressure profile atom types", &pressureProfileAtomTypes, 1);
1477    opts.range("pressureProfileAtomTypes", POSITIVE);
1478    opts.optional("pressureProfile", "pressureProfileAtomTypesFile", 
1479         "PDB files with pressure profile atom types" "default is the input PDB file", 
1480         PARSE_STRING);
1481    opts.optional("pressureProfile", "pressureProfileAtomTypesCol", 
1482         "Column in the pressureProfileAtomTypesFile with the atom types ",
1483         PARSE_STRING);
1484    opts.optionalB("pressureProfile", "pressureProfileEwald", 
1485        "Compute Ewald contribution to pressure profile",
1486        &pressureProfileEwaldOn, FALSE);
1487    opts.optional("pressureProfile", "pressureProfileEwaldX",
1488        "Ewald grid size X", &pressureProfileEwaldX, 10);
1489    opts.range("pressureProfileEwaldX", POSITIVE);
1490    opts.optional("pressureProfile", "pressureProfileEwaldY",
1491        "Ewald grid size Y", &pressureProfileEwaldY, 10);
1492    opts.range("pressureProfileEwaldY", POSITIVE);
1493    opts.optional("pressureProfile", "pressureProfileEwaldZ",
1494        "Ewald grid size Z", &pressureProfileEwaldZ, 10);
1495    opts.range("pressureProfileEwaldZ", POSITIVE);
1496
1497    /// accelerated MD parameters
1498    opts.optionalB("main", "accelMD", "Perform acclerated MD?", &accelMDOn, FALSE);
1499    opts.optional("accelMD", "accelMDFirstStep", "First accelMD step", &accelMDFirstStep, 0);
1500    opts.range("accelMDFirstStep", NOT_NEGATIVE);
1501    opts.optional("accelMD", "accelMDLastStep", "Last accelMD step", &accelMDLastStep, 0);
1502    opts.range("accelMDLastStep", NOT_NEGATIVE);
1503    opts.optional("accelMD", "accelMDOutFreq", "Frequency of accelMD output", &accelMDOutFreq, 1);
1504    opts.range("accelMDOutFreq", POSITIVE);
1505    opts.optionalB("accelMD", "accelMDdihe", "Apply boost to dihedral potential", &accelMDdihe, TRUE);
1506    opts.optionalB("accelMD", "accelMDDebugOn", "Debugging accelMD", &accelMDDebugOn, FALSE);
1507    opts.optional("accelMD", "accelMDE","E for AMD", &accelMDE);
1508    opts.units("accelMDE", N_KCAL);
1509    opts.optional("accelMD", "accelMDalpha","alpha for AMD", &accelMDalpha);
1510    opts.units("accelMDalpha", N_KCAL);
1511    opts.range("accelMDalpha", POSITIVE);
1512    opts.optionalB("accelMD", "accelMDdual", "Apply dual boost", &accelMDdual, FALSE);
1513    opts.optional("accelMDdual", "accelMDTE","E for total potential under accelMDdual mode", &accelMDTE);
1514    opts.units("accelMDTE", N_KCAL);
1515    opts.optional("accelMDdual", "accelMDTalpha","alpha for total potential under accelMDdual mode", &accelMDTalpha);
1516    opts.units("accelMDTalpha", N_KCAL);
1517    opts.range("accelMDTalpha", POSITIVE);
1518    // GaMD parameters
1519    opts.optionalB("accelMD", "accelMDG", "Perform Gaussian accelMD calculation?", &accelMDG, FALSE);
1520    opts.optional("accelMDG", "accelMDGiE", "Flag to set the mode iE in Gaussian accelMD", &accelMDGiE, 1);
1521    opts.optional("accelMDG", "accelMDGcMDSteps", "No. of cMD steps", &accelMDGcMDSteps, 1000000);
1522    opts.range("accelMDGcMDSteps", NOT_NEGATIVE);
1523    opts.optional("accelMDG", "accelMDGEquiSteps", "No. of equilibration steps after adding boost potential", &accelMDGEquiSteps, 1000000);
1524    opts.range("accelMDGEquiSteps", NOT_NEGATIVE);
1525    opts.require("accelMDG", "accelMDGcMDPrepSteps", "No. of preparation cMD steps", &accelMDGcMDPrepSteps, 200000);
1526    opts.range("accelMDGcMDPrepSteps", NOT_NEGATIVE);
1527    opts.require("accelMDG", "accelMDGEquiPrepSteps", "No. of preparation equilibration steps", &accelMDGEquiPrepSteps, 200000);
1528    opts.range("accelMDGEquiPrepSteps", NOT_NEGATIVE);
1529    opts.optional("accelMDG", "accelMDGSigma0P", "Upper limit of std of total potential", &accelMDGSigma0P, 6.0);
1530    opts.units("accelMDGSigma0P", N_KCAL);
1531    opts.range("accelMDGSigma0P", NOT_NEGATIVE);
1532    opts.optional("accelMDG", "accelMDGSigma0D", "Upper limit of std of dihedral potential", &accelMDGSigma0D, 6.0);
1533    opts.units("accelMDGSigma0D", N_KCAL);
1534    opts.range("accelMDGSigma0D", NOT_NEGATIVE);
1535    opts.optionalB("accelMDG", "accelMDGRestart", "Flag to set use restart file in Gaussian accelMD", &accelMDGRestart, FALSE);
1536    opts.require("accelMDGRestart", "accelMDGRestartFile", "Restart file name for Gaussian accelMD", accelMDGRestartFile);
1537    opts.optionalB("accelMDG", "accelMDGresetVaftercmd", "Flag to reset potential after accelMDGcMDSteps steps", 
1538            &accelMDGresetVaftercmd, FALSE);
1539
1540    // Adaptive Temperature Sampling (adaptTemp) parameters
1541    opts.optionalB("main", "adaptTempMD", "Perform adaptive temperature sampling", &adaptTempOn, FALSE);
1542    opts.optional("adaptTempMD", "adaptTempFirstStep", "First adaptTemp step", &adaptTempFirstStep, 0);
1543    opts.range("adaptTempFirstStep", NOT_NEGATIVE);
1544    opts.optional("adaptTempMD", "adaptTempLastStep", "Last adaptTemp step", &adaptTempLastStep, 0);
1545    opts.range("adaptTempLastStep", NOT_NEGATIVE);
1546    opts.optional("adaptTempMD", "adaptTempOutFreq", "Frequency of adaptTemp output", &adaptTempOutFreq, 10);
1547    opts.range("adaptTempOutFreq", POSITIVE);
1548    opts.optional("adaptTempMD", "adaptTempFreq", "Frequency of writing average energies to adaptTempOutFile", &adaptTempFreq, 10);
1549    opts.range("adaptTempFreq", POSITIVE);
1550    opts.optionalB("adaptTempMD", "adaptTempDebug", "Print debug output for adaptTemp", &adaptTempDebug, FALSE);
1551    opts.optional("adaptTempMD", "adaptTempTmin","Minimun temperature for adaptTemp", &adaptTempTmin);
1552    opts.units("adaptTempTmin", N_KELVIN);
1553    opts.range("adaptTempTmin", POSITIVE);
1554    opts.optional("adaptTempMD", "adaptTempTmax","Maximum temperature for adaptTemp", &adaptTempTmax);
1555    opts.units("adaptTempTmax", N_KELVIN);
1556    opts.range("adaptTempTmax", POSITIVE);
1557    opts.optional("adaptTempMD", "adaptTempBins","Number of bins to store average energies", &adaptTempBins,0);
1558    opts.range("adaptTempBins", NOT_NEGATIVE);
1559    opts.optional("adaptTempMD", "adaptTempDt", "Integration timestep for Temp. updates", &adaptTempDt, 0.0001);
1560    opts.units("adaptTempDt", N_FSEC);
1561    opts.range("adaptTempDt", NOT_NEGATIVE);
1562    opts.optional("adaptTempMD", "adaptTempAutoDt", "Average temperature update in percent of temperature range", &adaptTempAutoDt, 0.0);
1563    opts.range("adaptTempAutoDt", NOT_NEGATIVE);
1564    opts.optional("adaptTempMD", "adaptTempCgamma", "Adaptive bin averaging constant", &adaptTempCgamma, 0.1);
1565    opts.range("adaptTempCgamma", NOT_NEGATIVE);
1566    opts.optionalB("adaptTempMD","adaptTempLangevin","Send adaptTemp temperature to langevin thermostat",&adaptTempLangevin,TRUE);
1567    opts.optionalB("adaptTempMD","adaptTempRescaling","Send adaptTemp temperature to velocity rescaling thermostat", &adaptTempRescale,TRUE);
1568    opts.optional("adaptTempMD", "adaptTempInFile", "File containing restart information for adaptTemp", adaptTempInFile);
1569    opts.optional("adaptTempMD", "adaptTempRestartFile", "File for writing adaptTemp restart information", adaptTempRestartFile);
1570    opts.require("adaptTempRestartFile","adaptTempRestartFreq", "Frequency of writing restart file", &adaptTempRestartFreq,0);
1571    opts.range("adaptTempRestartFreq",NOT_NEGATIVE);
1572    opts.optionalB("adaptTempMD", "adaptTempRandom", "Randomly assign a temperature if we step out of range", &adaptTempRandom, FALSE);
1573 }
1574
1575 void SimParameters::config_parser_constraints(ParseOptions &opts) {
1576    
1577    ////  Fixed Atoms
1578    opts.optionalB("main", "fixedatoms", "Are there fixed atoms?",
1579     &fixedAtomsOn, FALSE);
1580    opts.optionalB("fixedatoms", "fixedAtomsForces",
1581      "Calculate forces between fixed atoms?  (Required to unfix during run.)",
1582      &fixedAtomsForces, FALSE);
1583    opts.optional("fixedatoms", "fixedAtomsFile", "PDB file with flags for "
1584      "fixed atoms (default is the PDB input file)",
1585      PARSE_STRING);
1586    opts.optional("fixedatoms", "fixedAtomsCol", "Column in the fixedAtomsFile "
1587      "containing the flags (nonzero means fixed);\n"
1588      "default is 'O'", PARSE_STRING);
1589    opts.optional("fixedatoms", "fixedAtomListFile", "the text input file for fixed atoms "
1590                  "used for parallel input IO", PARSE_STRING);
1591    opts.optionalB("fixedatoms", "fixedAtomsForceOutput",
1592      "Do we write out forces acting on fixed atoms?",
1593      &fixedAtomsForceOutput, FALSE);
1594
1595    ////  Harmonic Constraints
1596    opts.optionalB("main", "constraints", "Are harmonic constraints active?",
1597      &constraintsOn, FALSE);
1598    opts.require("constraints", "consexp", "Exponent for harmonic potential",
1599     &constraintExp, 2);
1600    opts.range("consexp", POSITIVE);
1601 #ifndef MEM_OPT_VERSION
1602    opts.require("constraints", "consref", "PDB file containing reference "
1603     "positions",
1604     PARSE_STRING);
1605    opts.require("constraints", "conskfile", "PDB file containing force "
1606     "constaints in one of the columns", PARSE_STRING);
1607    opts.require("constraints", "conskcol", "Column of conskfile to use "
1608     "for the force constants", PARSE_STRING);
1609 #else
1610    opts.require("constraints", "consAtomListFile", "the text input file for constrained atoms "
1611                  "used for parallel input IO", PARSE_STRING);
1612 #endif
1613    opts.require("constraints", "constraintScaling", "constraint scaling factor",
1614      &constraintScaling, 1.0);
1615    opts.range("constraintScaling", NOT_NEGATIVE);
1616
1617
1618
1619    //****** BEGIN selective restraints (X,Y,Z) changes
1620
1621    //// selective restraints (X,Y,Z) 
1622    opts.optionalB("constraints", "selectConstraints", 
1623    "Restrain only selected Cartesian components of the coordinates?",
1624      &selectConstraintsOn, FALSE);
1625    opts.optionalB("selectConstraints", "selectConstrX",  
1626    "Restrain X components of coordinates ", &constrXOn, FALSE);
1627    opts.optionalB("selectConstraints", "selectConstrY",  
1628    "Restrain Y components of coordinates ", &constrYOn, FALSE);
1629    opts.optionalB("selectConstraints", "selectConstrZ",  
1630    "Restrain Z components of coordinates ", &constrZOn, FALSE);
1631    //****** END selective restraints (X,Y,Z) changes
1632
1633    // spherical constraints
1634    opts.optionalB("constraints", "sphericalConstraints", 
1635    "Restrain only radial spherical component of the coordinates?",
1636      &sphericalConstraintsOn, FALSE);
1637    opts.optional("sphericalConstraints", "sphericalConstrCenter",
1638    "Center of spherical constraints", &sphericalConstrCenter);
1639  
1640    //****** BEGIN moving constraints changes 
1641
1642    //// Moving Harmonic Constraints
1643    opts.optionalB("constraints", "movingConstraints",
1644       "Are some of the constraints moving?", 
1645       &movingConstraintsOn, FALSE);
1646    opts.require("movingConstraints", "movingConsVel",
1647     "Velocity of the movement, A/timestep", &movingConsVel);
1648    //****** END moving constraints changes 
1649
1650    // BEGIN rotating constraints changes
1651    opts.optionalB("constraints", "rotConstraints",
1652       "Are the constraints rotating?", 
1653       &rotConstraintsOn, FALSE);
1654    opts.require("rotConstraints", "rotConsAxis",
1655     "Axis of rotation", &rotConsAxis);
1656    opts.require("rotConstraints", "rotConsPivot",
1657     "Pivot point of rotation", 
1658     &rotConsPivot);
1659    opts.require("rotConstraints", "rotConsVel",
1660     "Velocity of rotation, deg/timestep", &rotConsVel);
1661
1662    // END rotating constraints changes
1663
1664    // external command forces
1665    opts.optionalB("main", "extForces", "External command forces?",
1666       &extForcesOn, FALSE);
1667    opts.require("extForces", "extForcesCommand",
1668       "External forces command", extForcesCommand);
1669    opts.require("extForces", "extCoordFilename",
1670       "External forces coordinate filename", extCoordFilename);
1671    opts.require("extForces", "extForceFilename",
1672       "External forces force filename", extForceFilename);
1673
1674    
1675   // QM/MM forces
1676    opts.optionalB("main", "QMForces", "Apply QM forces?",
1677       &qmForcesOn, FALSE);
1678    opts.require("QMForces", "QMSoftware",
1679       "software whose format will be used for input/output", qmSoftware);
1680    opts.require("QMForces", "QMExecPath",
1681       "path to executable", qmExecPath);
1682    opts.optional("QMForces", "QMChargeMode",
1683       "type of QM atom charges gathered from the QM software", qmChrgModeS);
1684    opts.require("QMForces", "QMColumn",
1685       "column defining QM and MM regions", qmColumn);
1686    opts.require("QMForces", "QMBaseDir",
1687       "base path and name for QM input and output (preferably in memory)", qmBaseDir);
1688    opts.optional("QMForces", "QMConfigLine",
1689       "Configuration line for QM (multiple inputs allowed)", PARSE_MULTIPLES);
1690    opts.optional("QMForces", "QMParamPDB",
1691       "PDB with QM parameters", qmParamPDB);
1692    opts.optional("QMForces", "QMPrepProc",
1693       "initial preparation executable", qmPrepProc);
1694    opts.optional("QMForces", "QMSecProc",
1695       "secondary executable", qmSecProc);
1696    opts.optional("QMForces", "QMCharge",
1697       "charge of the QM group", PARSE_MULTIPLES);
1698    opts.optionalB("QMForces", "QMChargeFromPSF",
1699       "gets charge of the QM group form PSF values", &qmChrgFromPSF, FALSE);
1700    opts.optional("QMForces", "QMMult",
1701       "multiplicity of the QM group", PARSE_MULTIPLES);
1702    opts.optional("QMForces", "QMLinkElement",
1703       "element of link atom", PARSE_MULTIPLES);
1704    opts.optionalB("QMForces", "QMReplaceAll",
1705       "replace all NAMD forces with QM forces", &qmReplaceAll, FALSE);
1706    opts.optional("QMForces", "QMPCStride",
1707       "frequency of selection of point charges", &qmPCSelFreq, 1);
1708    opts.range("QMPCStride", POSITIVE);
1709    opts.optionalB("QMForces", "QMNoPntChrg",
1710       "no point charges will be passed to the QM system(s)", &qmNoPC, FALSE);
1711    opts.optionalB("QMForces", "QMElecEmbed",
1712       "activates electrostatic embedding", &qmElecEmbed, TRUE);
1713    opts.optionalB("QMForces", "QMVdWParams",
1714       "use special VdW parameters for QM atoms", &qmVDW, FALSE);
1715    opts.optional("QMForces", "QMBondColumn",
1716       "column defining QM-MM bomnds", qmBondColumn);
1717    opts.optionalB("QMForces", "QMBondDist",
1718       "values in QMBondColumn defines the distance of new link atom", &qmBondDist, FALSE);
1719    opts.optional("QMForces", "QMBondValueType",
1720       "type of value in bond column: len or ratio", qmBondValueTypeS);
1721    opts.optional("QMForces", "QMBondScheme",
1722       "type of treatment given to QM-MM bonds.", qmBondSchemeS);
1723    opts.optional("QMForces", "QMenergyStride",
1724       "frequency of QM specific energy output (every x steps)", &qmEnergyOutFreq, 1);
1725    opts.optional("QMForces", "QMOutStride",
1726       "frequency of QM specific charge output (every x steps)", &qmOutFreq, 0);
1727    opts.range("QMOutStride", NOT_NEGATIVE);
1728    opts.optional("QMForces", "QMPositionOutStride",
1729       "frequency of QM specific position output (every x steps)", &qmPosOutFreq, 0);
1730    opts.range("QMPositionOutStride", NOT_NEGATIVE);
1731    opts.optional("QMForces", "QMSimsPerNode",
1732       "QM executions per node", &qmSimsPerNode, 1);
1733    opts.range("QMSimsPerNode", POSITIVE);
1734    opts.optionalB("QMForces", "QMSwitching",
1735       "apply switching to point charges.", &qmPCSwitchOn, FALSE);
1736    opts.optional("QMForces", "QMSwitchingType",
1737       "How are charges scaled down to be presented to QM groups.", qmPCSwitchTypeS);
1738    opts.optional("QMForces", "QMPointChargeScheme",
1739       "type of treatment given to the total sum of point charges.", qmPCSchemeS);
1740    opts.optionalB("QMForces", "QMCustomPCSelection",
1741       "custom and fixed selection of point charges per QM group.", &qmCustomPCSel, FALSE);
1742    opts.optional("QMForces", "QMCustomPCFile",
1743       "file with a selection of point charges for a single QM group", PARSE_MULTIPLES);
1744    opts.optionalB("QMForces", "QMLiveSolventSel",
1745       "Continuously update the selection of solvent molecules in QM groups", &qmLSSOn, FALSE);
1746    opts.optional("QMForces", "QMLSSFreq",
1747       "frequency of QM water selection update", &qmLSSFreq, 100);
1748    opts.range("QMLSSFreq", POSITIVE);
1749    opts.optional("QMForces", "QMLSSResname",
1750       "residue name for the solvent molecules (TIP3).", qmLSSResname);
1751    opts.optional("QMForces", "QMLSSMode",
1752       "mode of selection of point solvent molecules", qmLSSModeS);
1753    opts.optional("QMForces", "QMLSSRef",
1754       "for COM mode, defines reference for COM distance calculation", PARSE_MULTIPLES);
1755    opts.optionalB("QMForces", "QMCSMD",
1756       "Do we use Conditional SMD option?", &qmCSMD, FALSE);
1757    opts.optional("QMForces", "QMCSMDFile",
1758                 "File for Conditional SMD information",qmCSMDFile);
1759    
1760    //print which bad contacts are being moved downhill
1761    opts.optionalB("main", "printBadContacts", "Print atoms with huge forces?",
1762       &printBadContacts, FALSE);
1763
1764    /* GBIS generalized born implicit solvent*/
1765
1766    opts.optionalB("main", "GBIS", "Use GB implicit solvent?",
1767       &GBISOn, FALSE);
1768    opts.optionalB("main", "GBISSer", "Use GB implicit solvent?",
1769       &GBISserOn, FALSE);
1770
1771    opts.optional("GBIS", "solventDielectric",
1772       "Solvent Dielectric", &solvent_dielectric, 78.5);
1773    opts.optional("GBIS", "intrinsicRadiusOffset",
1774       "Coulomb Radius Offset", &coulomb_radius_offset, 0.09);
1775    opts.optional("GBIS", "ionConcentration",
1776       "Ion Concentration", &ion_concentration, 0.2); //0.2 mol/L
1777    opts.optional("GBIS", "GBISDelta",
1778       "delta from GBOBC", &gbis_delta, 1.0); //0.8 or 1.0
1779    opts.optional("GBIS", "GBISBeta",
1780       "beta from GBOBC", &gbis_beta, 0.8);   //0.0 or 0.8
1781    opts.optional("GBIS", "GBISGamma",
1782       "gamma from GBOBC", &gbis_gamma, 4.85);//2.290912 or 4.85
1783    opts.optional("GBIS", "alphaCutoff",
1784       "cutoff for calculating effective born radius", &alpha_cutoff, 15);
1785    opts.optional("GBIS", "alphaMax",
1786       "maximum allowable born radius", &alpha_max, 30);
1787    opts.optional("GBIS", "fsMax",
1788       "maximum screened intrinsic radius", &fsMax, 1.728);
1789
1790    opts.optionalB("main", "SASA", "Use Linear Combination of Pairwise Overlaps (LCPO) for calculating SASA",
1791       &LCPOOn, FALSE);
1792    opts.optional("SASA", "surfaceTension",
1793       "Surfce Tension for SASA (kcal/mol/Ang^2)", &surface_tension, 0.005);
1794
1795    //****** BEGIN SMD constraints changes 
1796
1797    // SMD constraints
1798    opts.optionalB("main", "SMD",
1799       "Do we use SMD option?", 
1800       &SMDOn, FALSE);
1801    opts.require("SMD", "SMDVel",
1802                 "Velocity of the movement, A/timestep", &SMDVel);
1803    opts.range("SMDVel", NOT_NEGATIVE);
1804    opts.require("SMD", "SMDDir",
1805                 "Direction of movement", &SMDDir);
1806    opts.require("SMD", "SMDk",
1807                 "Elastic constant for SMD", &SMDk);
1808    opts.optional("SMD", "SMDk2",
1809                 "Transverse elastic constant for SMD", &SMDk2, 0);
1810    opts.range("SMDk", NOT_NEGATIVE);
1811    opts.range("SMDk2", NOT_NEGATIVE);
1812    opts.require("SMD", "SMDFile",
1813                 "File for SMD information",
1814                  SMDFile);
1815    opts.optional("SMD", "SMDOutputFreq",
1816                  "Frequency of output",
1817                  &SMDOutputFreq, 1);
1818    opts.range("SMDOutputFreq", POSITIVE);
1819    
1820    //****** END SMD constraints changes 
1821
1822    //****** BEGIN tabulated energies section
1823    opts.optionalB("main", "tabulatedEnergies", "Do we get energies from a table?", &tabulatedEnergies, FALSE);
1824 //   opts.require("tabulatedEnergies", "tableNumTypes","Number of types for energy tabulation", &tableNumTypes);
1825    opts.require("tabulatedEnergies", "tabulatedEnergiesFile", "File containing energy table", tabulatedEnergiesFile);
1826    opts.require("tabulatedEnergies", "tableInterpType", "Cubic or linear interpolation", tableInterpType);
1827
1828    // TMD parameters
1829    opts.optionalB("main", "TMD", "Perform Targeted MD?", &TMDOn, FALSE);
1830    opts.optional("TMD", "TMDk", "Elastic constant for TMD", &TMDk, 0); 
1831    opts.range("TMDk", NOT_NEGATIVE);
1832    opts.require("TMD", "TMDFile", "File for TMD information", TMDFile);
1833    opts.optionalB("TMD", "TMDDiffRMSD", "Restrain Difference between the RMSD from two structures", &TMDDiffRMSD, FALSE);
1834    opts.require("TMDDiffRMSD", "TMDFile2",  "Second file for TMD information", TMDFile2); 
1835     
1836    opts.optional("TMD", "TMDOutputFreq", "Frequency of TMD output", 
1837        &TMDOutputFreq, 1);
1838    opts.range("TMDOutputFreq", POSITIVE);
1839    opts.require("TMD", "TMDLastStep", "Last TMD timestep", &TMDLastStep);
1840    opts.range("TMDLastStep", POSITIVE);
1841    opts.optional("TMD", "TMDFirstStep", "First TMD step (default 0)", &TMDFirstStep, 0);
1842    opts.optional("TMD", "TMDInitialRMSD", "Target RMSD at first TMD step (default -1 to use initial coordinates)", &TMDInitialRMSD);
1843    TMDInitialRMSD = -1;
1844    opts.optional("TMD", "TMDFinalRMSD", "Target RMSD at last TMD step (default 0 )", &TMDFinalRMSD, 0);
1845    opts.range("TMDInitialRMSD", NOT_NEGATIVE);
1846    // End of TMD parameters
1847
1848    // Symmetry restraint parameters
1849    opts.optionalB("main", "symmetryRestraints", "Enable symmetry restraints?", &symmetryOn, FALSE); 
1850    opts.optional("symmetryRestraints", "symmetryk", "Elastic constant for symmetry restraints", &symmetryk, 0);
1851    opts.range("symmetryk", NOT_NEGATIVE);
1852    opts.optional("symmetryRestraints", "symmetrykfile", "PDB file specifying force contants on a per-atom basis", PARSE_MULTIPLES);
1853    opts.optionalB("symmetryRestraints", "symmetryScaleForces", "Scale applied forces over time?", &symmetryScaleForces, FALSE);
1854    opts.require("symmetryRestraints", "symmetryFile", "File for symmetry information", PARSE_MULTIPLES);
1855    opts.optional("symmetryRestraints", "symmetryMatrixFile", "File(s) for transfromation matrices", PARSE_MULTIPLES);
1856    opts.optional("symmetryRestraints", "symmetryLastStep", "Last symmetry timestep", &symmetryLastStep, -1);
1857    opts.optional("symmetryRestraints", "symmetryFirstStep", "First symmetry step (default 0)", &symmetryFirstStep, 0);
1858    opts.optional("symmetryRestraints", "symmetryLastFullStep", "Last full force symmetry timestep (default symmetryLastStep)", &symmetryLastFullStep, symmetryLastStep);
1859    opts.optional("symmetryRestraints", "symmetryFirstFullStep", "First full force symmetry step (default symmetryFirstStep)", &symmetryFirstFullStep, symmetryFirstStep);
1860   //End of symmetry restraint parameters.
1861
1862    ////  Global Forces / Tcl
1863    opts.optionalB("main", "tclForces", "Are Tcl global forces active?",
1864      &tclForcesOn, FALSE);
1865    opts.require("tclForces", "tclForcesScript",
1866      "Tcl script for global forces", PARSE_MULTIPLES);
1867
1868    ////  Boundary Forces / Tcl
1869    opts.optionalB("main", "tclBC", "Are Tcl boundary forces active?",
1870      &tclBCOn, FALSE);
1871    opts.require("tclBC", "tclBCScript",
1872      "Tcl script defining calcforces for boundary forces", PARSE_STRING);
1873    tclBCScript = 0;
1874    opts.optional("tclBC", "tclBCArgs", "Extra args for calcforces command",
1875      tclBCArgs);
1876    tclBCArgs[0] = 0;
1877
1878    ////  Global Forces / Misc
1879    opts.optionalB("main", "miscForces", "Are misc global forces active?",
1880      &miscForcesOn, FALSE);
1881    opts.optional("miscForces", "miscForcesScript",
1882      "script for misc forces", PARSE_MULTIPLES);
1883
1884    ////  Free Energy Perturbation
1885    opts.optionalB("main", "freeEnergy", "Perform free energy perturbation?",
1886      &freeEnergyOn, FALSE);
1887    opts.require("freeEnergy", "freeEnergyConfig",
1888      "Configuration file for free energy perturbation", PARSE_MULTIPLES);
1889
1890    ////  Constant Force
1891    opts.optionalB("main", "constantforce", "Apply constant force?",
1892      &consForceOn, FALSE);
1893    opts.optional("constantforce", "consForceFile",
1894        "Configuration file for constant forces", PARSE_STRING);
1895    opts.require("constantforce", "consForceScaling",
1896        "Scaling factor for constant forces", &consForceScaling, 1.0);
1897  
1898     //// Collective variables
1899     opts.optionalB("main", "colvars", "Is the colvars module enabled?",
1900       &colvarsOn, FALSE);
1901     opts.optional("colvars", "colvarsConfig",
1902       "configuration for the collective variables", PARSE_STRING);
1903     opts.optional("colvars", "colvarsInput",
1904       "input restart file for the collective variables", PARSE_STRING);
1905
1906 }
1907
1908 #ifdef OPENATOM_VERSION
1909 void SimParameters::config_parser_openatom(ParseOptions &opts) {
1910   opts.optionalB("main", "openatom", "OpenAtom active?", &openatomOn, FALSE);
1911   opts.require("openatom", "openatomDriverFile", "What config file specifies openatom input parameters", PARSE_STRING);
1912   opts.require("openatom", "openatomPhysicsFile", "What structure file specifies openatom input system", PARSE_STRING);
1913   opts.require("openatom", "openatomPdbFile", "NAMD input file defining QM and MM regions", PARSE_STRING);
1914    opts.optional("openatom", "openatomCol", "Column in the openatomPdb with the QM/MM flag", PARSE_STRING);
1915 }
1916 #endif // OPENATOM_VERSION
1917
1918 /* BEGIN gf */
1919 void SimParameters::config_parser_mgridforce(ParseOptions &opts) {
1920     //// Gridforce
1921     opts.optionalB("main", "mgridforce", "Is Multiple gridforce active?", 
1922                    &mgridforceOn, FALSE);
1923     opts.optional("mgridforce", "mgridforcevolts", "Is Gridforce using Volts/eV as units?",
1924                   PARSE_MULTIPLES);
1925     opts.require("mgridforce", "mgridforcescale", "Scale factor by which to multiply "
1926                  "grid forces", PARSE_MULTIPLES);
1927     opts.require("mgridforce", "mgridforcefile", "PDB file containing force "
1928                  "multipliers in one of the columns", PARSE_MULTIPLES);
1929     opts.require("mgridforce", "mgridforcecol", "Column of gridforcefile to "
1930                  "use for force multiplier", PARSE_MULTIPLES);
1931     opts.optional("mgridforce", "mgridforcechargecol", "Column of gridforcefile to "
1932                   "use for charge", PARSE_MULTIPLES);
1933     opts.require("mgridforce", "mgridforcepotfile", "Gridforce potential file",
1934                  PARSE_MULTIPLES);
1935     opts.optional("mgridforce", "mgridforcecont1", "Use continuous grid "
1936                    "in K1 direction?", PARSE_MULTIPLES);
1937     opts.optional("mgridforce", "mgridforcecont2", "Use continuous grid "
1938                    "in K2 direction?", PARSE_MULTIPLES);
1939     opts.optional("mgridforce", "mgridforcecont3", "Use continuous grid "
1940                    "in K3 direction?", PARSE_MULTIPLES);
1941     opts.optional("mgridforce", "mgridforcevoff", "Gridforce potential offsets",
1942                   PARSE_MULTIPLES);
1943     opts.optional("mgridforce", "mgridforcelite", "Use Gridforce Lite?",
1944                   PARSE_MULTIPLES);
1945     opts.optional("mgridforce", "mgridforcechecksize", "Check if grid exceeds PBC cell dimensions?", PARSE_MULTIPLES);
1946 }
1947
1948 void SimParameters::config_parser_gridforce(ParseOptions &opts) {
1949     //// Gridforce
1950     opts.optionalB("main", "gridforce", "Is Gridforce active?", 
1951                    &gridforceOn, FALSE);
1952     opts.optionalB("gridforce", "gridforcevolts", "Is Gridforce using Volts/eV as units?",
1953                    &gridforceVolts, FALSE);
1954     opts.require("gridforce", "gridforcescale", "Scale factor by which to multiply "
1955                  "grid forces", &gridforceScale);
1956     opts.require("gridforce", "gridforcefile", "PDB file containing force "
1957                  "multipliers in one of the columns", PARSE_STRING);
1958     opts.require("gridforce", "gridforcecol", "Column of gridforcefile to "
1959                  "use for force multiplier", PARSE_STRING);
1960     opts.optional("gridforce", "gridforcechargecol", "Column of gridforcefile to "
1961                   "use for charge", PARSE_STRING);
1962     opts.require("gridforce", "gridforcepotfile", "Gridforce potential file",
1963                  PARSE_STRING);
1964     opts.optionalB("gridforce", "gridforcecont1", "Use continuous grid "
1965                    "in A1 direction?", &gridforceContA1, FALSE);
1966     opts.optionalB("gridforce", "gridforcecont2", "Use continuous grid "
1967                    "in A2 direction?", &gridforceContA2, FALSE);
1968     opts.optionalB("gridforce", "gridforcecont3", "Use continuous grid "
1969                    "in A3 direction?", &gridforceContA3, FALSE);
1970     opts.optional("gridforce", "gridforcevoff", "Gridforce potential offsets",
1971                   &gridforceVOffset);
1972     opts.optionalB("gridforce", "gridforcelite", "Use Gridforce Lite?",
1973                    &gridforceLite, FALSE);
1974     opts.optionalB("gridforce", "gridforcechecksize", "Check if grid exceeds PBC cell dimensions?",
1975                    &gridforcechecksize, TRUE);
1976 }
1977 /* END gf */
1978
1979 void SimParameters::config_parser_movdrag(ParseOptions &opts) {
1980    //// moving drag
1981    opts.optionalB("main", "movDragOn", "Do we apply moving drag?",
1982       &movDragOn, FALSE);
1983    opts.require("movDragOn", "movDragFile",
1984       "Main moving drag PDB file", movDragFile);
1985    opts.require("movDragOn", "movDragCol",
1986       "Main moving drag PDB column", PARSE_STRING);
1987    opts.require("movDragOn", "movDragGlobVel",
1988       "Global moving drag velocity (A/step)", &movDragGlobVel);
1989    opts.require("movDragOn", "movDragVelFile",
1990       "Moving drag linear velocity file", movDragVelFile);
1991 }
1992
1993 void SimParameters::config_parser_rotdrag(ParseOptions &opts) {
1994    //// rotating drag
1995    opts.optionalB("main", "rotDragOn", "Do we apply rotating drag?",
1996       &rotDragOn, FALSE);
1997    opts.require("rotDragOn", "rotDragFile",
1998       "Main rotating drag PDB file", rotDragFile);
1999    opts.require("rotDragOn", "rotDragCol",
2000       "Main rotating drag PDB column", PARSE_STRING);
2001    opts.require("rotDragOn", "rotDragAxisFile",
2002       "Rotating drag axis file", rotDragAxisFile);
2003    opts.require("rotDragOn", "rotDragPivotFile",
2004       "Rotating drag pivot point file", rotDragPivotFile);
2005    opts.require("rotDragOn", "rotDragGlobVel",
2006       "Global rotating drag angular velocity (deg/step)", &rotDragGlobVel);
2007    opts.require("rotDragOn", "rotDragVelFile",
2008       "Rotating drag angular velocity file", rotDragVelFile);
2009    opts.require("rotDragOn", "rotDragVelCol",
2010       "Rotating drag angular velocity column", PARSE_STRING);
2011 }
2012
2013 void SimParameters::config_parser_constorque(ParseOptions &opts) {
2014    //// "constant" torque
2015    opts.optionalB("main", "consTorqueOn", "Do we apply \"constant\" torque?",
2016       &consTorqueOn, FALSE);
2017    opts.require("consTorqueOn", "consTorqueFile",
2018       "Main \"constant\" torque PDB file", consTorqueFile);
2019    opts.require("consTorqueOn", "consTorqueCol",
2020       "Main \"constant\" torque PDB column", PARSE_STRING);
2021    opts.require("consTorqueOn", "consTorqueAxisFile",
2022       "\"Constant\" torque axis file", consTorqueAxisFile);
2023    opts.require("consTorqueOn", "consTorquePivotFile",
2024       "\"Constant\" torque pivot point file", consTorquePivotFile);
2025    opts.require("consTorqueOn", "consTorqueGlobVal",
2026       "Global \"constant\" torque value (Kcal/(mol*A^2))", &consTorqueGlobVal);
2027    opts.require("consTorqueOn", "consTorqueValFile",
2028       "\"constant\" torque factors file", consTorqueValFile);
2029    opts.require("consTorqueOn", "consTorqueValCol",
2030       "\"constant\" torque factors column", PARSE_STRING);
2031 }
2032
2033 void SimParameters::config_parser_boundary(ParseOptions &opts) {
2034     
2035    //// Spherical Boundary Conditions
2036    opts.optionalB("main", "sphericalBC", "Are spherical boundary counditions "
2037       "active?", &sphericalBCOn, FALSE);
2038    opts.require("sphericalBC", "sphericalBCCenter",
2039      "Center of spherical boundaries", &sphericalCenter);
2040    opts.require("sphericalBC", "sphericalBCr1", "Radius for first sphere "
2041      "potential", &sphericalBCr1);
2042    opts.range("sphericalBCr1", POSITIVE);
2043    opts.units("sphericalBCr1", N_ANGSTROM);
2044    opts.require("sphericalBC", "sphericalBCk1", "Force constant for first "
2045     "sphere potential (+ is an inward force, - outward)",
2046     &sphericalBCk1);
2047    opts.units("sphericalBCk1", N_KCAL);
2048    opts.optional("sphericalBC", "sphericalBCexp1", "Exponent for first "
2049     "sphere potential", &sphericalBCexp1, 2);
2050    opts.range("sphericalBCexp1", POSITIVE);
2051    
2052    opts.optional("sphericalBCr1", "sphericalBCr2", "Radius for second sphere "
2053      "potential", &sphericalBCr2);
2054    opts.range("sphericalBCr2", POSITIVE);
2055    opts.units("sphericalBCr2", N_ANGSTROM);
2056    opts.require("sphericalBCr2", "sphericalBCk2", "Force constant for second "
2057     "sphere potential (+ is an inward force, - outward)",
2058     &sphericalBCk2);
2059    opts.units("sphericalBCk2", N_KCAL);
2060    opts.optional("sphericalBCr2", "sphericalBCexp2", "Exponent for second "
2061     "sphere potential", &sphericalBCexp2, 2);
2062    opts.range("sphericalBCexp2", POSITIVE);
2063
2064    /////////////// Cylindrical Boundary Conditions
2065    opts.optionalB("main", "cylindricalBC", "Are cylindrical boundary counditions "
2066                   "active?", &cylindricalBCOn, FALSE);
2067    opts.require("cylindricalBC", "cylindricalBCr1", "Radius for first cylinder "
2068                  "potential", &cylindricalBCr1);
2069    opts.range("cylindricalBCr1", POSITIVE);
2070    opts.units("cylindricalBCr1", N_ANGSTROM);
2071    opts.require("cylindricalBC", "cylindricalBCk1", "Force constant for first "
2072                 "cylinder potential (+ is an inward force, - outward)",
2073                 &cylindricalBCk1);
2074    opts.units("cylindricalBCk1", N_KCAL);
2075    opts.optional("cylindricalBC", "cylindricalBCexp1", "Exponent for first "
2076                 "cylinder potential", &cylindricalBCexp1, 2);
2077    opts.range("cylindricalBCexp1", POSITIVE);
2078
2079
2080 // additions beyond those already found in spherical parameters    JJU
2081    opts.optional("cylindricalBC", "cylindricalBCAxis", "Cylinder axis (defaults to x)",
2082     PARSE_STRING);
2083    opts.require("cylindricalBC", "cylindricalBCCenter",
2084      "Center of cylindrical boundaries", &cylindricalCenter);
2085    opts.require ("cylindricalBC", "cylindricalBCl1", "Length of first cylinder",
2086                  &cylindricalBCl1);
2087    opts.range("cylindricalBCl1", POSITIVE);
2088    opts.units("cylindricalBCl1", N_ANGSTROM);
2089    opts.optional ("cylindricalBCl1", "cylindricalBCl2", "Length of second cylinder",
2090                   &cylindricalBCl2);
2091    opts.range ("cylindricalBCl2", POSITIVE);
2092    opts.units ("cylindricalBCl2", N_ANGSTROM);
2093 // end  additions
2094
2095    opts.optional("cylindricalBCr1", "cylindricalBCr2", "Radius for second cylinder "
2096                  "potential", &cylindricalBCr2);
2097    opts.range("cylindricalBCr2", POSITIVE);
2098    opts.units("cylindricalBCr2", N_ANGSTROM);
2099    opts.require("cylindricalBCr2", "cylindricalBCk2", "Force constant for second "
2100                 "cylinder potential (+ is an inward force, - outward)",
2101                 &cylindricalBCk2);
2102    opts.units("cylindricalBCk2", N_KCAL);
2103    opts.optional("cylindricalBCr2", "cylindricalBCexp2", "Exponent for second "
2104                 "cylinder potential", &cylindricalBCexp2, 2);
2105    opts.range("cylindricalBCexp2", POSITIVE);
2106
2107    ///////////////  Electric field options
2108    opts.optionalB("main", "eFieldOn", "Should an electric field be applied",
2109                  &eFieldOn, FALSE);
2110    opts.optionalB("eFieldOn", "eFieldNormalized", "Is eField vector scaled by cell basis vectors?",
2111                  &eFieldNormalized, FALSE);
2112    opts.require("eFieldOn", "eField", "Electric field vector", &eField);
2113    opts.optional("eFieldOn", "eFieldFreq", "Electric field frequency", &eFieldFreq);
2114    opts.optional("eFieldOn", "eFieldPhase", "Electric field phase", &eFieldPhase);
2115
2116       ///////////////  Stir options
2117    opts.optionalB("main", "stirOn", "Should stirring torque be applied",
2118                  &stirOn, FALSE);
2119    opts.optional("stirOn", "stirFilename", "PDB file with flags for "
2120      "stirred atoms (default is the PDB input file)",
2121                  PARSE_STRING);
2122    opts.optional("stirOn", "stirredAtomsCol", "Column in the stirredAtomsFile "
2123                  "containing the flags (nonzero means fixed);\n"
2124                  "default is 'O'", PARSE_STRING);
2125    opts.require("stirOn", "stirStartingTheta", "Stir starting theta offset", &stirStartingTheta);
2126    opts.require("stirOn", "stirK", "Stir force harmonic spring constant", &stirK);
2127    //should make this optional, compute from firsttimestep * stirVel
2128    opts.require("stirOn", "stirVel", "Stir angular velocity (deg/timestep)", &stirVel);
2129    opts.require("stirOn", "stirAxis", "Stir axis (direction vector)", &stirAxis);
2130    opts.require("stirOn", "stirPivot", "Stir pivot point (coordinate)", &stirPivot);
2131
2132     //////////  Extra bonds options
2133    opts.optionalB("main", "extraBonds",
2134                 "Should extra bonded forces be applied",
2135                  &extraBondsOn, FALSE);
2136    opts.optional("extraBonds", "extraBondsFile",
2137                 "file with list of extra bonds",
2138                  PARSE_MULTIPLES);
2139    opts.optionalB("extraBonds", "extraBondsCosAngles",
2140                 "Should extra angles be cosine-based to match ancient bug",
2141                 &extraBondsCosAngles, TRUE);
2142
2143 }
2144
2145 void SimParameters::config_parser_misc(ParseOptions &opts) {
2146    
2147    ///////////////  Load balance options
2148    opts.optional("main", "ldBalancer", "Load balancer",
2149      loadBalancer);
2150    opts.optional("main", "ldbStrategy", "Load balancing strategy",
2151      loadStrategy);
2152    opts.optional("main", "ldbPeriod", "steps between load balancing", 
2153      &ldbPeriod);
2154    opts.range("ldbPeriod", POSITIVE);
2155    opts.optional("main", "firstLdbStep", "when to start load balancing",
2156      &firstLdbStep);
2157    opts.range("firstLdbStep", POSITIVE);
2158    opts.optional("main", "lastLdbStep", "when to stop load balancing",
2159      &lastLdbStep);
2160    opts.range("lastLdbStep", POSITIVE);
2161    opts.optional("main", "hybridGroupSize", "Hybrid load balancing group size",
2162      &hybridGroupSize);
2163    opts.optional("main", "ldbBackgroundScaling",
2164      "background load scaling", &ldbBackgroundScaling);
2165    opts.range("ldbBackgroundScaling", NOT_NEGATIVE);
2166    opts.optional("main", "ldbPMEBackgroundScaling",
2167      "PME node background load scaling", &ldbPMEBackgroundScaling);
2168    opts.range("ldbPMEBackgroundScaling", NOT_NEGATIVE);
2169    opts.optional("main", "ldbHomeBackgroundScaling",
2170      "home node background load scaling", &ldbHomeBackgroundScaling);
2171    opts.range("ldbHomeBackgroundScaling", NOT_NEGATIVE);
2172    opts.optional("main", "ldbRelativeGrainsize",
2173      "fraction of average load per compute", &ldbRelativeGrainsize, 0.);
2174    opts.range("ldbRelativeGrainsize", NOT_NEGATIVE);
2175    
2176    opts.optional("main", "traceStartStep", "when to start tracing", &traceStartStep);
2177    opts.range("traceStartStep", POSITIVE);
2178    opts.optional("main", "numTraceSteps", "the number of timesteps to be traced", &numTraceSteps);
2179    opts.range("numTraceSteps", POSITIVE);
2180  
2181 #ifdef MEASURE_NAMD_WITH_PAPI
2182    opts.optionalB("main", "papiMeasure", "whether use PAPI to measure performacne", &papiMeasure, FALSE);
2183    opts.optional("main", "papiMeasureStartStep", "when to measure performacne using PAPI", &papiMeasureStartStep);
2184    opts.range("papiMeasureStartStep", POSITIVE);
2185    opts.optional("main", "numPapiMeasureSteps", "the number of timesteps to be measured using PAPI", &numPapiMeasureSteps);
2186    opts.range("numPapiMeasureSteps", POSITIVE);
2187 #endif
2188
2189    opts.optionalB("main", "outputMaps", "whether to dump compute map and patch map for analysis just before load balancing", &outputMaps, FALSE);
2190    opts.optionalB("main", "benchTimestep", "whether to do benchmarking timestep in which case final file output is disabled", &benchTimestep, FALSE);
2191    opts.optional("main", "useCkLoop", "whether to use CkLoop library to parallelize a loop in a function like OpenMP", &useCkLoop,
2192     #if CMK_SMP && USE_CKLOOP
2193      ( CkNumPes() < 2 * CkNumNodes() ? 0 : CKLOOP_CTRL_PME_FORWARDFFT ) );
2194     #else
2195      0);
2196     #endif
2197    opts.range("useCkLoop", NOT_NEGATIVE);
2198
2199    opts.optionalB("main", "simulateInitialMapping", "whether to study the initial mapping scheme", &simulateInitialMapping, FALSE);
2200    opts.optional("main", "simulatedPEs", "the number of PEs to be used for studying initial mapping", &simulatedPEs);
2201    opts.range("simulatedPEs", POSITIVE);
2202    opts.optional("main", "simulatedNodeSize", "the node size to be used for studying initial mapping", &simulatedNodeSize);
2203    opts.range("simulatedNodeSize", POSITIVE);
2204    opts.optionalB("main", "disableTopology", "ignore torus information during patch placement", &disableTopology, FALSE);
2205    opts.optionalB("main", "verboseTopology", "print torus information during patch placement", &verboseTopology, FALSE);
2206
2207    opts.optionalB("main", "ldbUnloadPME", "no load on PME nodes",
2208      &ldbUnloadPME, FALSE);
2209    opts.optionalB("main", "ldbUnloadZero", "no load on pe zero",
2210      &ldbUnloadZero, FALSE);
2211    opts.optionalB("main", "ldbUnloadOne", "no load on pe one",
2212      &ldbUnloadOne, FALSE);
2213    opts.optionalB("main", "ldbUnloadOutputPEs", "no load on output PEs",
2214      &ldbUnloadOutputPEs, FALSE);
2215    opts.optionalB("main", "noPatchesOnZero", "no patches on pe zero",
2216      &noPatchesOnZero, FALSE);
2217    opts.optionalB("main", "noPatchesOnOutputPEs", "no patches on Output PEs",
2218      &noPatchesOnOutputPEs, FALSE);
2219    opts.optionalB("main", "noPatchesOnOne", "no patches on pe one",
2220      &noPatchesOnOne, FALSE);
2221    opts.optionalB("main", "useCompressedPsf", "The structure file psf is in the compressed format",
2222                   &useCompressedPsf, FALSE);
2223    opts.optionalB("main", "genCompressedPsf", "Generate the compressed version of the psf file",
2224                   &genCompressedPsf, FALSE);
2225    opts.optionalB("main", "usePluginIO", "Use the plugin I/O to load the molecule system", 
2226                   &usePluginIO, FALSE);   
2227    opts.optionalB("main", "mallocTest", "test how much memory all PEs can allocate", 
2228                   &mallocTest, FALSE);   
2229    opts.optionalB("main", "printExclusions", "print exclusion lists to stdout", 
2230                   &printExclusions, FALSE);   
2231    opts.optional("main", "proxySendSpanningTree", "using spanning tree to send proxies",
2232                   &proxySendSpanningTree, -1);
2233    opts.optional("main", "proxyRecvSpanningTree", "using spanning tree to receive proxies",
2234                   &proxyRecvSpanningTree, 0);  // default off due to memory leak -1);
2235    opts.optional("main", "proxyTreeBranchFactor", "the branch factor when building a spanning tree",
2236                   &proxyTreeBranchFactor, 0);  // actual default in ProxyMgr.C
2237    opts.optionalB("main", "twoAwayX", "half-size patches in 1st dimension",
2238      &twoAwayX, -1);
2239    opts.optionalB("main", "twoAwayY", "half-size patches in 2nd dimension",
2240      &twoAwayY, -1);
2241    opts.optionalB("main", "twoAwayZ", "half-size patches in 3rd dimension",
2242      &twoAwayZ, -1);
2243    opts.optional("main", "maxPatches", "maximum patch count", &maxPatches, -1);
2244
2245    /////  Restart timestep option
2246    opts.optional("main", "firsttimestep", "Timestep to start simulation at",
2247      &firstTimestep, 0);
2248    opts.range("firsttimestep", NOT_NEGATIVE);
2249  
2250    /////  Test mode options
2251    opts.optionalB("main", "test", "Perform self-tests rather than simulation",
2252                 &testOn, FALSE);
2253    opts.optionalB("main", "commOnly", "Do not evaluate forces or integrate",
2254                 &commOnly, FALSE);
2255
2256    opts.optionalB("main", "statsOn", "counters in machine layer",
2257                 &statsOn, FALSE);
2258    ///////////////  hydrogen bond computation options
2259    opts.optionalB("main", "hbonds", "Use explicit hydrogen bond term",
2260                  &HydrogenBonds, FALSE);
2261    opts.optionalB("hbonds","hbAntecedents","Include Antecedent in hbond term",
2262                  &useAntecedent, TRUE);
2263    opts.optional("hbonds","hbAAexp","Hbond AA-A-H angle cos exponential",
2264                  &aaAngleExp, 2);
2265    opts.optional("hbonds","hbHAexp","Hbond D-H-A angle cos exponential",
2266                  &haAngleExp, 4);
2267    opts.optional("hbonds","hbDistAexp","Hbond A-D dist attractive exponential",
2268                  &distAttExp, 4);
2269    opts.optional("hbonds","hbDistRexp","Hbond A-D dist repulstive exponential",
2270                  &distRepExp, 6);
2271    opts.optional("hbonds","hbCutoffAngle","Hbond D-H-A cutoff angle",
2272                  &dhaCutoffAngle, 100.0);
2273    opts.range("hbCutoffAngle", NOT_NEGATIVE);
2274    opts.optional("hbonds","hbOnAngle","Hbond D-H-A switch function on angle",
2275                  &dhaOnAngle, 60.0);
2276    opts.range("hbOnAngle", NOT_NEGATIVE);
2277    opts.optional("hbonds","hbOffAngle","Hbond D-H-A switch function off angle",
2278                  &dhaOffAngle, 80.0);
2279    opts.range("hbOffAngle", NOT_NEGATIVE);
2280    opts.optional("hbonds","hbCutoffDist","Hbond A-D cutoff distance",
2281                  &daCutoffDist, 7.5);
2282    opts.range("hbCutoffDist", POSITIVE);
2283    opts.units("hbCutoffDist", N_ANGSTROM);
2284    opts.optional("hbonds","hbOnDist","Hbond A-D switch function on distance",
2285                  &daOnDist, 5.5);
2286    opts.range("hbOnDist", POSITIVE);
2287    opts.units("hbOnDist", N_ANGSTROM);
2288    opts.optional("hbonds","hbOffDist","Hbond A-D switch function off distance",
2289                  &daOffDist, 6.5);
2290    opts.range("hbOffDist", POSITIVE);
2291    opts.units("hbOffDist", N_ANGSTROM);
2292
2293    // IMD options
2294    opts.optionalB("main","IMDon","Connect using IMD?",&IMDon, FALSE);
2295    opts.require("IMDon","IMDport", "Port to which to bind", &IMDport);
2296    opts.range("IMDport",POSITIVE);
2297    opts.require("IMDon","IMDfreq", "Frequency at which to report", &IMDfreq);
2298    opts.range("IMDfreq",POSITIVE);
2299    opts.optionalB("IMDon","IMDwait","Pause until IMD connection?",&IMDwait,
2300      FALSE);
2301    opts.optionalB("IMDon","IMDignore","Ignore any user input?",&IMDignore,
2302      FALSE);
2303    opts.optionalB("IMDon","IMDignoreForces","Ignore forces ONLY?",&IMDignoreForces,
2304      FALSE);
2305    // Maximum Partition options
2306    opts.optional("ldBalancer", "maxSelfPart", 
2307      "maximum number of self partitions in one patch", &maxSelfPart, 20);
2308    opts.range("maxSelfPart",POSITIVE);
2309    opts.optional("ldBalancer", "maxPairPart", 
2310      "maximum number of pair partitions in one patch", &maxPairPart, 8);
2311    opts.range("maxPairPart",POSITIVE);
2312    opts.optional("ldBalancer", "numAtomsSelf", 
2313                  "maximum number of atoms in one self compute distribution", 
2314                  &numAtomsSelf, 154);
2315    opts.range("numAtomsSelf",NOT_NEGATIVE);
2316
2317    opts.optional("ldBalancer", "numAtomsSelf2", 
2318                  "maximum number of atoms in one self compute distribution", 
2319                  &numAtomsSelf2, 154);
2320    opts.range("numAtomsSelf2",NOT_NEGATIVE);
2321
2322    opts.optional("ldBalancer", "numAtomsPair", 
2323                  "maximum number of atoms in one pair compute distribution", 
2324                  &numAtomsPair, 318);
2325    opts.range("numAtomsPair",NOT_NEGATIVE);
2326    opts.optional("ldBalancer", "numAtomsPair2", 
2327                "maximum number of atoms in one pair compute distribution", 
2328                &numAtomsPair2, 637);
2329    opts.range("numAtomsPair2",NOT_NEGATIVE);
2330    opts.optional("main", "minAtomsPerPatch", 
2331                "minimum average atoms per patch", 
2332                &minAtomsPerPatch, 40);
2333    opts.range("minAtomsPerPatch",NOT_NEGATIVE);
2334
2335    // Maximum exclusion flags per atom
2336    opts.optional("main", "maxExclusionFlags", 
2337      "maximum number of exclusion flags per atom", &maxExclusionFlags, 256);
2338    opts.range("maxExclusionFlags",POSITIVE);
2339
2340    // Bonded interactions on GPU
2341    opts.optional("main", "bondedCUDA", "Bitmask for calculating bonded interactions on GPU", &bondedCUDA, 255);
2342
2343    // Automatically disable individual CUDA kernels that are
2344    // incompatible with simulation options.
2345    // Set FALSE to manually control kernel use for development.
2346    opts.optionalB("main", "useCUDAdisable", "Disable kernels to maintain feature compatibility with CUDA", &useCUDAdisable, TRUE);
2347
2348    // MIC specific parameters
2349    opts.optional("main", "mic_unloadMICPEs", "Indicates whether or not the load balancer should unload PEs driving Xeon Phi cards", &mic_unloadMICPEs, 1);
2350    opts.optional("main", "mic_singleKernel", "Set to non-zero to have all MIC work to be placed in a single kernel", &mic_singleKernel, 1);
2351    opts.optional("main", "mic_deviceThreshold", "Threshold to use for directing computes to Xeon Phi devices", &mic_deviceThreshold, -1);
2352    opts.optional("main", "mic_hostSplit", "DMK - reserved", &mic_hostSplit, -1);
2353    opts.optional("main", "mic_numParts_self_p1", "MIC-Specific NumParts SELF Parameter 1", &mic_numParts_self_p1, -1);
2354    opts.optional("main", "mic_numParts_pair_p1", "MIC-Specific NumParts PAIR Parameter 1", &mic_numParts_pair_p1, -1);
2355    opts.optional("main", "mic_numParts_pair_p2", "MIC-Specific NumParts PAIR Parameter 2", &mic_numParts_pair_p2, -1);
2356    opts.range("mic_unloadMICPEs", NOT_NEGATIVE);
2357    opts.range("mic_singleKernel", NOT_NEGATIVE);
2358 }
2359
2360 void SimParameters::readExtendedSystem(const char *filename, Lattice *latptr) {
2361
2362      if ( ! latptr ) {
2363        iout << iINFO << "EXTENDED SYSTEM FILE   " << filename << "\n" << endi;
2364      }
2365
2366      ifstream xscFile(filename);
2367      if ( ! xscFile ) NAMD_die("Unable to open extended system file.\n");
2368
2369      char labels[1024];
2370      do {
2371        if ( ! xscFile ) NAMD_die("Error reading extended system file.\n");
2372        xscFile.getline(labels,1023);
2373      } while ( strncmp(labels,"#$LABELS ",9) );
2374
2375      int a_x, a_y, a_z, b_x, b_y, b_z, c_x, c_y, c_z;
2376      a_x = a_y = a_z = b_x = b_y = b_z = c_x = c_y = c_z = -1;
2377      int o_x, o_y, o_z, s_u, s_v, s_w, s_x, s_y, s_z;
2378      o_x = o_y = o_z = s_u = s_v = s_w = s_x = s_y = s_z = -1;
2379
2380      int pos = 0;
2381      char *l_i = labels + 8;
2382      while ( *l_i ) {
2383        if ( *l_i == ' ' ) { ++l_i; continue; }
2384        char *l_i2;
2385        for ( l_i2 = l_i; *l_i2 && *l_i2 != ' '; ++l_i2 );
2386        if ( (l_i2 - l_i) == 3 && (l_i[1] == '_') ) {
2387          if (l_i[0] == 'a' && l_i[2] == 'x') a_x = pos;
2388          if (l_i[0] == 'a' && l_i[2] == 'y') a_y = pos;
2389          if (l_i[0] == 'a' && l_i[2] == 'z') a_z = pos;
2390          if (l_i[0] == 'b' && l_i[2] == 'x') b_x = pos;
2391          if (l_i[0] == 'b' && l_i[2] == 'y') b_y = pos;
2392          if (l_i[0] == 'b' && l_i[2] == 'z') b_z = pos;
2393          if (l_i[0] == 'c' && l_i[2] == 'x') c_x = pos;
2394          if (l_i[0] == 'c' && l_i[2] == 'y') c_y = pos;
2395          if (l_i[0] == 'c' && l_i[2] == 'z') c_z = pos;
2396          if (l_i[0] == 'o' && l_i[2] == 'x') o_x = pos;
2397          if (l_i[0] == 'o' && l_i[2] == 'y') o_y = pos;
2398          if (l_i[0] == 'o' && l_i[2] == 'z') o_z = pos;
2399          if (l_i[0] == 's' && l_i[2] == 'u') s_u = pos;
2400          if (l_i[0] == 's' && l_i[2] == 'v') s_v = pos;
2401          if (l_i[0] == 's' && l_i[2] == 'w') s_w = pos;
2402          if (l_i[0] == 's' && l_i[2] == 'x') s_x = pos;
2403          if (l_i[0] == 's' && l_i[2] == 'y') s_y = pos;
2404          if (l_i[0] == 's' && l_i[2] == 'z') s_z = pos;
2405        }
2406        ++pos;
2407        l_i = l_i2;
2408      }
2409      int numpos = pos;
2410
2411      for ( pos = 0; pos < numpos; ++pos ) {
2412        double tmp;
2413        xscFile >> tmp;
2414        if ( ! xscFile ) NAMD_die("Error reading extended system file.\n");
2415        if ( pos == a_x ) cellBasisVector1.x = tmp;
2416        if ( pos == a_y ) cellBasisVector1.y = tmp;
2417        if ( pos == a_z ) cellBasisVector1.z = tmp;
2418        if ( pos == b_x ) cellBasisVector2.x = tmp;
2419        if ( pos == b_y ) cellBasisVector2.y = tmp;
2420        if ( pos == b_z ) cellBasisVector2.z = tmp;
2421        if ( pos == c_x ) cellBasisVector3.x = tmp;
2422        if ( pos == c_y ) cellBasisVector3.y = tmp;
2423        if ( pos == c_z ) cellBasisVector3.z = tmp;
2424        if ( pos == o_x ) cellOrigin.x = tmp;
2425        if ( pos == o_y ) cellOrigin.y = tmp;
2426        if ( pos == o_z ) cellOrigin.z = tmp;
2427        if ( pos == s_u ) strainRate2.x = tmp;
2428        if ( pos == s_v ) strainRate2.y = tmp;
2429        if ( pos == s_w ) strainRate2.z = tmp;
2430        if ( pos == s_x ) strainRate.x = tmp;
2431        if ( pos == s_y ) strainRate.y = tmp;
2432        if ( pos == s_z ) strainRate.z = tmp;
2433      }
2434
2435    if ( latptr ) {
2436      Lattice test;
2437      test.set(cellBasisVector1,cellBasisVector2,cellBasisVector3,cellOrigin);
2438     
2439      if ( test.a_p() && ! lattice.a_p() ) {
2440        NAMD_die("cellBasisVector1 added during atom reinitialization");
2441      }
2442      if ( lattice.a_p() && ! test.a_p() ) {
2443        NAMD_die("cellBasisVector1 dropped during atom reinitialization");
2444      }
2445      if ( test.b_p() && ! lattice.b_p() ) {
2446        NAMD_die("cellBasisVector2 added during atom reinitialization");
2447      }
2448      if ( lattice.b_p() && ! test.b_p() ) {
2449        NAMD_die("cellBasisVector2 dropped during atom reinitialization");
2450      }
2451      if ( test.c_p() && ! lattice.c_p() ) {
2452        NAMD_die("cellBasisVector3 added during atom reinitialization");
2453      }
2454      if ( lattice.c_p() && ! test.c_p() ) {
2455        NAMD_die("cellBasisVector3 dropped during atom reinitialization");
2456      }
2457
2458      latptr->set(cellBasisVector1,cellBasisVector2,cellBasisVector3,cellOrigin);
2459    }
2460
2461 }
2462
2463 #ifdef MEM_OPT_VERSION
2464 //This global var is defined in mainfunc.C
2465 extern char *gWorkDir;
2466 #endif
2467
2468 void SimParameters::check_config(ParseOptions &opts, ConfigList *config, char *&cwd) {
2469    
2470    int len;    //  String length
2471    StringList *current; //  Pointer to config option list
2472
2473 #ifdef MEM_OPT_VERSION
2474    char *namdWorkDir = NULL;
2475 #endif
2476
2477   if ( opts.defined("obsolete") ) {
2478     iout << iWARN <<
2479       "\"obsolete\" defined, silently ignoring obsolete options\n" << endi;
2480   }
2481
2482    //  Take care of cwd processing
2483    if (opts.defined("cwd"))
2484    {
2485     //  First allocate and get the cwd value
2486     current = config->find("cwd");
2487
2488     len = strlen(current->data);
2489
2490     if ( CHDIR(current->data) )
2491     {
2492       NAMD_die("chdir() to given cwd failed!");
2493     } else {
2494       iout << iINFO << "Changed directory to " << current->data << "\n" << endi;
2495     }
2496
2497     if (current->data[len-1] != PATHSEP)
2498       len++;
2499
2500     cwd = new char[len+1];
2501
2502     strcpy(cwd, current->data);
2503
2504     if (current->data[strlen(current->data)-1] != PATHSEP)
2505       strcat(cwd, PATHSEPSTR);
2506    }
2507
2508 #ifdef MEM_OPT_VERSION
2509    if(cwd!=NULL)namdWorkDir = cwd;     
2510    else namdWorkDir = gWorkDir;
2511    int dirlen = strlen(namdWorkDir);
2512    //only support the path representation on UNIX-like platforms
2513    char *tmpDir;
2514    if(namdWorkDir[dirlen-1]=='/'){
2515      tmpDir = new char[dirlen+1];
2516      tmpDir[dirlen] = 0;
2517    }else{
2518      tmpDir = new char[dirlen+2];
2519      tmpDir[dirlen]='/';
2520      tmpDir[dirlen+1]=0;
2521    } 
2522    memcpy(tmpDir, namdWorkDir, dirlen);
2523    namdWorkDir = tmpDir;
2524  //finished recording the per atom files, free the space for gWorkDir
2525    delete [] gWorkDir;
2526 #endif
2527
2528
2529    // Don't try to specify coordinates with pluginIO
2530    if ( usePluginIO && opts.defined("coordinates") ) {
2531      NAMD_die("Separate coordinates file not allowed with plugin IO, coordinates will be taken from structure file.");
2532    }
2533
2534    // If it's not AMBER||GROMACS, then "coordinates", "structure"
2535    // and "parameters" must be specified.
2536    if (!amberOn && !gromacsOn) {
2537 #ifndef MEM_OPT_VERSION
2538      if (useCompressedPsf)
2539        NAMD_die("useCompressedPsf requires memory-optimized build!");
2540      if (!usePluginIO && !genCompressedPsf && !opts.defined("coordinates"))
2541        NAMD_die("coordinates not found in the configuration file!");
2542 #else
2543      if(!usePluginIO && !opts.defined("bincoordinates")) {
2544        NAMD_die("bincoordinates not found in the configuration file for the memory optimized version!");
2545      }
2546      if(!usePluginIO && opts.defined("coordinates")) {
2547        NAMD_die("coordinates not allowed in the configuration file for the memory optimized version!");
2548      }
2549 #endif
2550      if (!opts.defined("structure"))
2551        NAMD_die("structure not found in the configuration file!");
2552      if (!opts.defined("parameters"))
2553        NAMD_die("parameters not found in the configuration file!");
2554    }
2555    
2556    // In any case, there should be either "coordinates" or
2557    // "ambercoor", but not both
2558    if (opts.defined("coordinates") && opts.defined("ambercoor"))
2559      NAMD_die("Cannot specify both coordinates and ambercoor!");
2560 #ifndef MEM_OPT_VERSION
2561    if (!genCompressedPsf && !opts.defined("coordinates") && !opts.defined("ambercoor")
2562        && !opts.defined("grocoorfile") && !usePluginIO)
2563      NAMD_die("Coordinate file not found!");
2564 #endif
2565
2566    //  Make sure that both a temperature and a velocity PDB were
2567    //  specified
2568    if (opts.defined("temperature") &&
2569        (opts.defined("velocities") || opts.defined("binvelocities")) ) 
2570    {
2571       NAMD_die("Cannot specify both an initial temperature and a velocity file");
2572    }
2573
2574 #ifdef MEM_OPT_VERSION
2575 //record the absolute file name for binAtomFile, binCoorFile and binVelFile etc.
2576    binAtomFile = NULL;
2577    binCoorFile = NULL;
2578    binVelFile = NULL;   
2579    binRefFile = NULL;
2580
2581    char *curfile = NULL;
2582    dirlen = strlen(namdWorkDir);
2583    current = config->find("structure");;
2584    curfile = current->data;
2585    int filelen = strlen(curfile);
2586    if(*curfile == '/' || *curfile=='~') {
2587      //check whether it is an absolute path
2588      //WARNING: Only works on Unix-like platforms!
2589      //Needs to fix on Windows platform.
2590      //-Chao Mei     
2591      //adding 5 because of ".bin"+"\0"
2592      binAtomFile = new char[filelen+5];
2593      memcpy(binAtomFile, curfile, filelen);
2594      memcpy(binAtomFile+filelen, ".bin", 4);
2595      binAtomFile[filelen+4] = 0;
2596    }else{
2597      binAtomFile = new char[dirlen+filelen+5];
2598      memcpy(binAtomFile, namdWorkDir, dirlen);
2599      memcpy(binAtomFile+dirlen, curfile, filelen);
2600      memcpy(binAtomFile+dirlen+filelen, ".bin", 4);
2601      binAtomFile[dirlen+filelen+4] = 0;
2602    }
2603
2604    current = config->find("bincoordinates");
2605    curfile = current->data;
2606    filelen = strlen(curfile);
2607    if(*curfile == '/' || *curfile=='~') {
2608      binCoorFile = new char[filelen+1];
2609      memcpy(binCoorFile, curfile, filelen);
2610      binCoorFile[filelen] = 0;
2611    }else{
2612      binCoorFile = new char[dirlen+filelen+1];
2613      memcpy(binCoorFile, namdWorkDir, dirlen);
2614      memcpy(binCoorFile+dirlen, curfile, filelen);
2615      binCoorFile[dirlen+filelen] = 0;
2616    }
2617
2618    if(opts.defined("binvelocities")){
2619      current = config->find("binvelocities");
2620      curfile = current->data;
2621      filelen = strlen(curfile);
2622      if(*curfile == '/' || *curfile=='~') {
2623        binVelFile = new char[filelen+1];
2624        memcpy(binVelFile, curfile, filelen);
2625        binVelFile[filelen] = 0;
2626      }else{
2627        binVelFile = new char[dirlen+filelen+1];
2628        memcpy(binVelFile, namdWorkDir, dirlen);
2629        memcpy(binVelFile+dirlen, curfile, filelen);
2630        binVelFile[dirlen+filelen] = 0;
2631      }
2632    }
2633
2634    if(opts.defined("binrefcoords")){
2635      current = config->find("binrefcoords");
2636      curfile = current->data;
2637      filelen = strlen(curfile);
2638      if(*curfile == '/' || *curfile=='~') {
2639        binRefFile = new char[filelen+1];
2640        memcpy(binRefFile, curfile, filelen);
2641        binRefFile[filelen] = 0;
2642      }else{
2643        binRefFile = new char[dirlen+filelen+1];
2644        memcpy(binRefFile, namdWorkDir, dirlen);
2645        memcpy(binRefFile+dirlen, curfile, filelen);
2646        binRefFile[dirlen+filelen] = 0;
2647      }
2648    }
2649
2650    //deal with output file name to make it absolute path for parallel output
2651    if(outputFilename[0] != '/' && outputFilename[0]!='~') {
2652      filelen = strlen(outputFilename);
2653      char *tmpout = new char[filelen];
2654      memcpy(tmpout, outputFilename, filelen);
2655      CmiAssert(filelen+dirlen <= 120); //leave 8 chars for file suffix
2656      memcpy(outputFilename, namdWorkDir, dirlen);
2657      memcpy(outputFilename+dirlen, tmpout, filelen);
2658      outputFilename[filelen+dirlen] = 0;     
2659      delete [] tmpout;
2660    }
2661
2662    if ( dcdFrequency && opts.defined("dcdfile") &&
2663         dcdFilename[0] != '/' && dcdFilename[0]!='~' ) {
2664      filelen = strlen(dcdFilename);
2665      char *tmpout = new char[filelen];
2666      memcpy(tmpout, dcdFilename, filelen);
2667      CmiAssert(filelen+dirlen <= 120); //leave 8 chars for file suffix
2668      memcpy(dcdFilename, namdWorkDir, dirlen);
2669      memcpy(dcdFilename+dirlen, tmpout, filelen);
2670      dcdFilename[filelen+dirlen] = 0;     
2671      delete [] tmpout;
2672    }
2673
2674    if ( velDcdFrequency && opts.defined("veldcdfile") &&
2675         velDcdFilename[0] != '/' && velDcdFilename[0]!='~' ) {
2676      filelen = strlen(velDcdFilename);
2677      char *tmpout = new char[filelen];
2678      memcpy(tmpout, velDcdFilename, filelen);
2679      CmiAssert(filelen+dirlen <= 120); //leave 8 chars for file suffix
2680      memcpy(velDcdFilename, namdWorkDir, dirlen);
2681      memcpy(velDcdFilename+dirlen, tmpout, filelen);
2682      velDcdFilename[filelen+dirlen] = 0;     
2683      delete [] tmpout;
2684    }
2685
2686    if ( forceDcdFrequency && opts.defined("forcedcdfile") &&
2687         forceDcdFilename[0] != '/' && forceDcdFilename[0]!='~' ) {
2688      filelen = strlen(forceDcdFilename);
2689      char *tmpout = new char[filelen];
2690      memcpy(tmpout, forceDcdFilename, filelen);
2691      CmiAssert(filelen+dirlen <= 120); //leave 8 chars for file suffix
2692      memcpy(forceDcdFilename, namdWorkDir, dirlen);
2693      memcpy(forceDcdFilename+dirlen, tmpout, filelen);
2694      forceDcdFilename[filelen+dirlen] = 0;     
2695      delete [] tmpout;
2696    }
2697
2698    if ( restartFrequency && opts.defined("restartname") &&
2699         restartFilename[0] != '/' && restartFilename[0]!='~' ) {
2700      filelen = strlen(restartFilename);
2701      char *tmpout = new char[filelen];
2702      memcpy(tmpout, restartFilename, filelen);
2703      CmiAssert(filelen+dirlen <= 120); //leave 8 chars for file suffix
2704      memcpy(restartFilename, namdWorkDir, dirlen);
2705      memcpy(restartFilename+dirlen, tmpout, filelen);
2706      restartFilename[filelen+dirlen] = 0;     
2707      delete [] tmpout;
2708    }
2709
2710    delete [] namdWorkDir;
2711
2712    if (opts.defined("numinputprocs")) { 
2713      if(numinputprocs > CkNumPes()) {
2714        iout << iWARN << "The number of input processors exceeds the total number of processors. Resetting to half of the number of total processors.\n" << endi;
2715        numinputprocs = (CkNumPes()>>1)+(CkNumPes()&1);
2716      }
2717    }
2718
2719    if (opts.defined("numoutputprocs")) {        
2720      if(numoutputprocs > CkNumPes()) {
2721        iout << iWARN << "The number of output processors exceeds the total number of processors. Resetting to half of the number of total processors.\n" << endi;
2722        numoutputprocs = (CkNumPes()>>1)+(CkNumPes()&1);
2723      }
2724    }
2725
2726 #ifndef OUTPUT_SINGLE_FILE
2727 #error OUTPUT_SINGLE_FILE not defined!
2728 #endif
2729
2730    #if !OUTPUT_SINGLE_FILE
2731    //create directories for multi-file output scheme   
2732    create_output_directories("coor");
2733    create_output_directories("vel");
2734    if(dcdFrequency) {
2735            create_output_directories("dcd");
2736            if(opts.defined("dcdfile")){
2737                    iout << iWARN << "The dcd file output has been changed to directory: " << outputFilename << ".\n" << endi; 
2738            }
2739    }
2740    if (velDcdFrequency) {
2741            create_output_directories("veldcd");
2742            if(opts.defined("veldcdfile")){       
2743                    iout << iWARN << "The veldcd file output has been changed to directory: " << outputFilename << ".\n" << endi;
2744            }
2745    }
2746    if (forceDcdFrequency) {
2747            create_output_directories("forcedcd");
2748            if(opts.defined("forcedcdfile")){       
2749                    iout << iWARN << "The forcedcd file output has been changed to directory: " << outputFilename << ".\n" << endi;
2750            }
2751    }
2752    #endif
2753 #endif
2754
2755    if (! opts.defined("auxFile")) {
2756      strcpy(auxFilename,outputFilename);
2757      strcat(auxFilename,".aux");
2758    }
2759
2760    //  Check for frequencies
2761    if (dcdFrequency) {
2762      if (! opts.defined("dcdfile")) {
2763        strcpy(dcdFilename,outputFilename);
2764        strcat(dcdFilename,".dcd");
2765      }
2766    } else {
2767      dcdFilename[0] = STRINGNULL;
2768    }
2769
2770    if (velDcdFrequency) {
2771      if (! opts.defined("veldcdfile")) {
2772        strcpy(velDcdFilename,outputFilename);
2773        strcat(velDcdFilename,".veldcd");
2774      }
2775    } else {
2776      velDcdFilename[0] = STRINGNULL;
2777    }
2778    
2779    if (forceDcdFrequency) {
2780      if (! opts.defined("forcedcdfile")) {
2781        strcpy(forceDcdFilename,outputFilename);
2782        strcat(forceDcdFilename,".forcedcd");
2783      }
2784    } else {
2785      forceDcdFilename[0] = STRINGNULL;
2786    }
2787    
2788    if (xstFrequency) {
2789      if (! opts.defined("xstfile")) {
2790        strcpy(xstFilename,outputFilename);
2791        strcat(xstFilename,".xst");
2792      }
2793    } else {
2794      xstFilename[0] = STRINGNULL;
2795    }
2796
2797    if (restartFrequency) {
2798      if (! opts.defined("restartname")) {
2799        strcpy(restartFilename,outputFilename);
2800        if ( ! restartSave ) strcat(restartFilename,".restart");
2801      }
2802    } else {
2803      restartFilename[0] = STRINGNULL;
2804      restartSave = FALSE;
2805      binaryRestart = FALSE;
2806    }
2807
2808    if (storeComputeMap || loadComputeMap) {
2809      if (! opts.defined("computeMapFile")) {
2810        strcpy(computeMapFilename,"computeMapFile");
2811        strcat(computeMapFilename,".txt");
2812      }
2813    }
2814
2815
2816    if (!amberOn)
2817    { //****** BEGIN CHARMM/XPLOR type changes
2818      //// set default
2819      if (!paraTypeXplorOn && !paraTypeCharmmOn) 
2820      {
2821        paraTypeXplorOn = TRUE;
2822      }
2823      //// make sure that there is just one type of input parameters specified
2824      if (paraTypeXplorOn && paraTypeCharmmOn) 
2825      {
2826        NAMD_die("Please specify either XPLOR or CHARMM format for parameters!");
2827      }
2828      //****** END CHARMM/XPLOR type changes
2829    }
2830
2831    
2832    //  If minimization isn't on, must have a temp or velocity
2833    if (!(minimizeOn||minimizeCGOn) && !opts.defined("temperature") && 
2834        !opts.defined("velocities") && !opts.defined("binvelocities") ) 
2835    {
2836       NAMD_die("Must have either an initial temperature or a velocity file");
2837    }
2838
2839    if (minimizeOn||minimizeCGOn) { initialTemp = 0.0; }
2840    if (opts.defined("velocities") || opts.defined("binvelocities") )
2841    {
2842   initialTemp = -1.0;
2843    }
2844
2845    ///// periodic cell parameters
2846
2847    if ( opts.defined("extendedSystem") ) readExtendedSystem(config->find("extendedSystem")->data);
2848
2849 #ifdef MEM_OPT_VERSION
2850    if ( LJcorrection ) {
2851       NAMD_die("LJ tail corrections not yet available for memory optimized builds");
2852    }
2853 #endif
2854
2855    if ( LJcorrection && ! cellBasisVector3.length2() ) {
2856      NAMD_die("Can't use LJ tail corrections without periodic boundary conditions!");
2857    }
2858
2859    if ( cellBasisVector3.length2() && ! cellBasisVector2.length2() ) {
2860      NAMD_die("Used cellBasisVector3 without cellBasisVector2!");
2861    }
2862
2863    if ( cellBasisVector2.length2() && ! cellBasisVector1.length2() ) {
2864      NAMD_die("Used cellBasisVector2 without cellBasisVector1!");
2865    }
2866
2867    if ( cellOrigin.length2() && ! cellBasisVector1.length2() ) {
2868      NAMD_die("Used cellOrigin without cellBasisVector1!");
2869    }
2870
2871    lattice.set(cellBasisVector1,cellBasisVector2,cellBasisVector3,cellOrigin);
2872
2873    if (! opts.defined("DCDunitcell")) {
2874       dcdUnitCell = lattice.a_p() && lattice.b_p() && lattice.c_p();
2875    }
2876
2877    char s[129];
2878
2879    ///// cylindricalBC stuff
2880    if ( ! opts.defined("cylindricalBCAxis") )
2881    {
2882       cylindricalBCAxis = 'x';
2883    }
2884    else
2885    {
2886      opts.get("cylindricalBCAxis", s);
2887
2888      if (!strcasecmp(s, "x"))
2889      {
2890       cylindricalBCAxis = 'x';
2891      }
2892      else if (!strcasecmp(s, "y"))
2893      {
2894       cylindricalBCAxis = 'y';
2895      }
2896      else if (!strcasecmp(s, "z"))
2897      {
2898       cylindricalBCAxis = 'z';
2899      }
2900    else
2901      {
2902       char err_msg[128];
2903
2904       sprintf(err_msg, "Illegal value '%s' for 'cylindricalBCAxis' in configuration file", s);
2905       NAMD_die(err_msg);
2906      }
2907    }
2908
2909    if (!opts.defined("splitPatch"))
2910    {
2911      splitPatch = SPLIT_PATCH_HYDROGEN;
2912    }
2913    else
2914    {
2915      opts.get("splitPatch", s);
2916      if (!strcasecmp(s, "position"))
2917        splitPatch = SPLIT_PATCH_POSITION;
2918      else if (!strcasecmp(s,"hydrogen"))
2919        splitPatch = SPLIT_PATCH_HYDROGEN;
2920      else
2921      {
2922        char err_msg[129];
2923        sprintf(err_msg, 
2924           "Illegal value '%s' for 'splitPatch' in configuration file", 
2925        s);
2926        NAMD_die(err_msg);
2927      }
2928    }
2929
2930    ///// exclude stuff
2931    opts.get("exclude", s);
2932
2933    if (!strcasecmp(s, "none"))
2934    {
2935       exclude = NONE;
2936       splitPatch = SPLIT_PATCH_POSITION;
2937    }
2938    else if (!strcasecmp(s, "1-2"))
2939    {
2940       exclude = ONETWO;
2941       splitPatch = SPLIT_PATCH_POSITION;
2942    }
2943    else if (!strcasecmp(s, "1-3"))
2944    {
2945       exclude = ONETHREE;
2946    }
2947    else if (!strcasecmp(s, "1-4"))
2948    {
2949       exclude = ONEFOUR;
2950    }
2951    else if (!strcasecmp(s, "scaled1-4"))
2952    {
2953       exclude = SCALED14;
2954    }
2955    else
2956    {
2957       char err_msg[128];
2958
2959       sprintf(err_msg, "Illegal value '%s' for 'exclude' in configuration file",
2960    s);
2961       NAMD_die(err_msg);
2962    }
2963
2964    if (scale14 != 1.0 && exclude != SCALED14)
2965    {
2966       iout << iWARN << "Exclude is not scaled1-4; 1-4scaling ignored.\n" << endi;
2967    }
2968
2969    // water model stuff
2970    if (!opts.defined("waterModel")) {
2971      watmodel = WAT_TIP3;
2972    } else {
2973      opts.get("waterModel", s);
2974      if (!strncasecmp(s, "tip4", 4)) {
2975        iout << iINFO << "Using TIP4P water model.\n" << endi;
2976        watmodel = WAT_TIP4;
2977      } else if (!strncasecmp(s, "tip3", 4)) {
2978        iout << iINFO << "Using TIP3P water model.\n" << endi;
2979        watmodel = WAT_TIP3;
2980      } else if (!strncasecmp(s, "swm4", 4)) {
2981        iout << iINFO << "Using SWM4-DP water model.\n" << endi;
2982        watmodel = WAT_SWM4;
2983      } else {
2984        char err_msg[128];
2985        sprintf(err_msg,
2986            "Illegal value %s for 'waterModel' in configuration file", s);
2987        NAMD_die(err_msg);
2988      }
2989    }
2990    if (watmodel == WAT_SWM4 && !drudeOn) {
2991      NAMD_die("Must have 'drudeOn' enabled to use SWM4-DP water model.");
2992    }
2993    if (drudeOn && watmodel != WAT_SWM4) {
2994      watmodel = WAT_SWM4;
2995      iout << iWARN
2996        << "Setting water model to 'swm4' (SWM4-DP) for Drude polarization.\n"
2997        << endi;
2998    }
2999
3000    // Drude water model uses "lonepairs"
3001    if (watmodel == WAT_SWM4) {
3002      lonepairs = TRUE;
3003    }
3004
3005    // Added by JLai -- 8.2.11 -- Checks if Go method is defined
3006    if (goForcesOn) {
3007      iout << iINFO << "Go forces are on\n" << endi;
3008    // Added by JLai -- 6.3.11 -- Checks if Go method is defined
3009      int * gomethod = &goMethod;
3010      if (!opts.defined("GoMethod")) {
3011        *gomethod = 0;
3012        //     printf("GO METHOD IS NOT DEFINED SO WE'LL SET IT TO SOME WEIRD VALUE\n");
3013      } else {
3014        opts.get("GoMethod",s);
3015        // printf("GO METHOD IS DEFINED SO WE'LL PRINT IT OUT: %s\n",s);
3016        *gomethod = atoi(s);
3017      }
3018      if (!strcasecmp(s, "matrix")) {
3019        goMethod = 1;
3020        //GoMethod = GO_MATRIX;
3021      } else if (!strcasecmp(s, "faster")) {
3022        goMethod = 2;
3023        //GoMethod = GO_FASTER;
3024      } else if (!strcasecmp(s, "lowmem")) {
3025        goMethod = 3;
3026        //GoMethod = GO_LOWMEM;
3027      }
3028      else {
3029        char err_msg[129];     
3030        sprintf(err_msg,
3031                "Illegal value '%s' for 'GoMethod' in configuration file",
3032                s);
3033        NAMD_die(err_msg);
3034      }
3035    }  // End of NAMD code to check goMethod
3036    // End of Port -- JL
3037
3038    //  Get multiple timestep integration scheme
3039    if (!opts.defined("MTSAlgorithm"))
3040    {
3041   MTSAlgorithm = VERLETI;
3042    }
3043    else
3044    {
3045   opts.get("MTSAlgorithm", s);
3046
3047   if (!strcasecmp(s, "naive"))
3048   {
3049     MTSAlgorithm = NAIVE;
3050   }
3051   else if (!strcasecmp(s, "constant"))
3052   {
3053     MTSAlgorithm = NAIVE;
3054   }
3055   else if (!strcasecmp(s, "impulse"))
3056   {
3057     MTSAlgorithm = VERLETI;
3058   }
3059   else if (!strcasecmp(s, "verleti"))
3060   {
3061     MTSAlgorithm = VERLETI;
3062   }
3063   else
3064   {
3065     char err_msg[129];
3066
3067     sprintf(err_msg, 
3068        "Illegal value '%s' for 'MTSAlgorithm' in configuration file", 
3069        s);
3070     NAMD_die(err_msg);
3071   }
3072    }
3073
3074    //  Get the long range force splitting specification
3075    if (!opts.defined("longSplitting"))
3076    {
3077   longSplitting = C1;
3078    }
3079    else
3080    {
3081   opts.get("longSplitting", s);
3082   if (!strcasecmp(s, "sharp"))
3083     longSplitting = SHARP;
3084   else if (!strcasecmp(s, "xplor"))
3085     longSplitting = XPLOR;
3086   else if (!strcasecmp(s, "c1"))
3087     longSplitting = C1;
3088   else if (!strcasecmp(s, "c2"))
3089     longSplitting = C2;
3090   else
3091   {
3092     char err_msg[129];
3093
3094     sprintf(err_msg, 
3095        "Illegal value '%s' for 'longSplitting' in configuration file", 
3096        s);
3097     NAMD_die(err_msg);
3098   }
3099    }
3100
3101    // take care of rigid bond options
3102    if (!opts.defined("rigidBonds"))
3103    {
3104       rigidBonds = RIGID_NONE;
3105    }
3106    else
3107    {
3108       opts.get("rigidBonds", s); 
3109       if (!strcasecmp(s, "all"))
3110       {
3111           rigidBonds = RIGID_ALL;
3112       }
3113       else if (!strcasecmp(s, "water"))
3114       {
3115            rigidBonds = RIGID_WATER;
3116       } 
3117       else if (!strcasecmp(s, "none"))
3118       {
3119            rigidBonds = RIGID_NONE;
3120       } 
3121       else
3122       {
3123         char err_msg[256];
3124         sprintf(err_msg, 
3125           "Illegal value '%s' for 'rigidBonds' in configuration file", s);
3126         NAMD_die(err_msg);
3127       }
3128    }
3129
3130    // TIP4P and SWM4-DP water models require rigid water
3131    if ((watmodel == WAT_TIP4 || watmodel == WAT_SWM4)
3132        && rigidBonds == RIGID_NONE) {
3133      char err_msg[256];
3134      sprintf(err_msg,
3135          "Water model %s requires rigidBonds set to \"all\" or \"water\"",
3136          (watmodel == WAT_TIP4 ? "TIP4P" : "SWM4-DP"));
3137      NAMD_die(err_msg);
3138    }
3139
3140    //  Take care of switching stuff
3141    if (switchingActive)
3142    {
3143
3144      if (!opts.defined("switchDist")) {
3145        NAMD_die("switchDist must be defined when switching is enabled");
3146      }
3147
3148      if ( (switchingDist>cutoff) || (switchingDist<0) )
3149      {
3150        char err_msg[129];
3151
3152        sprintf(err_msg, 
3153          "switchDist muct be between 0 and cutoff, which is %f", cutoff);
3154        NAMD_die(err_msg);
3155      }
3156
3157    }
3158
3159    if ( martiniSwitching )
3160    {
3161      if ( ! switchingActive ) 
3162      { 
3163        NAMD_die("martiniSwitching requires switching");
3164      }
3165      if ( vdwForceSwitching ) 
3166      { 
3167        NAMD_die("martiniSwitching and vdwForceSwitching are exclusive to one another. Select only one."); 
3168      }
3169      if ( dielectric != 15.0 && ! martiniDielAllow ) 
3170      {
3171        iout << iWARN << "USE DIELECTRIC OF 15.0 WITH MARTINI.\n";
3172        iout << iWARN << "SETTING dielectric 15.0\n";
3173        iout << iWARN << "FOR NON-STANDARD DIELECTRIC WITH MARTINI, SET: martiniDielAllow on\n";
3174        dielectric = 15.0;
3175      }
3176      if ( ! cosAngles )
3177      {
3178        iout << iWARN << "USE COSINE BASED ANGLES WITH MARTINI.\n";
3179        iout << iWARN << "SETTING cosAngles on\n";
3180        cosAngles = TRUE;
3181      }
3182      if ( PMEOn )
3183      {
3184        NAMD_die("Do not use Particle Mesh Ewald with Martini.  Set: PME off");
3185      }
3186      if ( MSMOn )
3187      {
3188        NAMD_die("Do not use Multilevel Summation Method with Martini.  Set: MSM off");
3189      }
3190      if ( FMMOn )
3191      {
3192        NAMD_die("Do not use Fast Multipole Method with Martini.  Set: FMM off");
3193      }
3194
3195    }
3196
3197
3198    if (!opts.defined("pairlistDist"))
3199    {
3200   pairlistDist = cutoff;
3201    }
3202    else if (pairlistDist < cutoff)
3203    {
3204   NAMD_die("pairlistDist must be >= cutoff distance");
3205    }
3206
3207    patchDimension = pairlistDist;
3208
3209    if ( splitPatch == SPLIT_PATCH_HYDROGEN ) {
3210      patchDimension += hgroupCutoff;
3211    }
3212
3213    BigReal defaultMargin = 0.0;
3214    if (berendsenPressureOn || langevinPistonOn) {
3215       defaultMargin = ( useFlexibleCell ? 0.06 : 0.03 ) * patchDimension;
3216    }
3217    if ( margin == XXXBIGREAL ) {
3218      margin = defaultMargin;
3219    }
3220    if ( defaultMargin != 0.0 && margin == 0.0 ) {
3221      margin = defaultMargin;
3222      iout << iWARN << "ALWAYS USE NON-ZERO MARGIN WITH CONSTANT PRESSURE!\n";
3223      iout << iWARN << "CHANGING MARGIN FROM 0 to " << margin << "\n" << endi;
3224    }
3225
3226    patchDimension += margin;
3227
3228     //ensure patch can handle alpha_cutoff for gbis
3229     if (GBISOn) {
3230       //Check compatibility
3231       if (fullDirectOn) {
3232         NAMD_die("GBIS not compatible with FullDirect");
3233       }
3234       if (PMEOn) {
3235         NAMD_die("GBIS not compatible with PME");
3236       }
3237       if (MSMOn) {
3238         NAMD_die("GBIS not compatible with MSM");
3239       }
3240       if (FMMOn) {
3241         NAMD_die("GBIS not compatible with FMM");
3242       }
3243       if (alchOn) {
3244         NAMD_die("GBIS not compatible with Alchemical Transformations");
3245       }
3246       if (lesOn) {
3247         NAMD_die("GBIS not compatible with Locally Enhanced Sampling");
3248       }
3249       if (FMAOn) {
3250         NAMD_die("GBIS not compatible with FMA");
3251       }
3252       if (drudeOn) {
3253         NAMD_die("GBIS not compatible with Drude Polarization");
3254       }
3255
3256       if (alpha_cutoff > patchDimension) {
3257         patchDimension = alpha_cutoff; 
3258       }
3259       //calculate kappa
3260       BigReal tmp = (initialTemp > 0) ? initialTemp : 300;
3261       kappa = 50.29216*sqrt(ion_concentration/solvent_dielectric/tmp);
3262       /*magic number = 1/sqrt(eps0*kB/(2*nA*e^2*1000))*/
3263     } // GBISOn
3264
3265     if (LCPOOn) {
3266 #ifdef MEM_OPT_VERSION
3267       NAMD_die("SASA not yet available for memory optimized builds");
3268 #endif
3269       if ( lattice.volume() > 0 ) {
3270         NAMD_die("SASA does not yet support periodic boundary conditions.");
3271       }
3272       //LCPO requires patches to be at least 16.2Ang in each dimension
3273       // twoAway[XYZ} is ignored for now
3274     }
3275
3276    //  Turn on global integration if not explicitly specified
3277
3278    if ( dihedralOn ) globalOn = TRUE;
3279
3280 #ifdef NAMD_CUDA
3281    if (loweAndersenOn) {
3282        NAMD_die("Lowe-Andersen dynamics not compatible with CUDA at this time");
3283    }
3284 #endif
3285    // END LA
3286
3287    // BEGIN LA
3288    if (loweAndersenOn && (langevinOn || tCoupleOn))
3289    {
3290       NAMD_die("Lowe-Andersen dynamics, Langevin dynamics and temperature coupling are mutually exclusive dynamics modes");
3291    }
3292    // END LA
3293
3294    if (tCoupleOn && opts.defined("rescaleFreq") )
3295    {
3296       NAMD_die("Temperature coupling and temperature rescaling are mutually exclusive");
3297    }
3298
3299    if (globalOn && CkNumPes() > 1)
3300    {
3301       NAMD_die("Global integration does not run in parallel (yet).");
3302    }
3303
3304    if (COLDOn && langevinOn)
3305    {
3306       NAMD_die("COLD and Langevin dynamics are mutually exclusive dynamics modes");
3307    }
3308    if (COLDOn && minimizeOn)
3309    {
3310       NAMD_die("COLD and minimization are mutually exclusive dynamics modes");
3311    }
3312    if (COLDOn && tCoupleOn)
3313    {
3314       NAMD_die("COLD and temperature coupling are mutually exclusive dynamics modes");
3315    }
3316    if (COLDOn && opts.defined("rescaleFreq"))
3317    {
3318       NAMD_die("COLD and velocity rescaling are mutually exclusive dynamics modes");
3319    }
3320
3321    if (splitPatch == SPLIT_PATCH_POSITION && mollyOn )
3322    {
3323       NAMD_die("splitPatch hydrogen is required for MOLLY");
3324    }
3325
3326    if (splitPatch == SPLIT_PATCH_POSITION && rigidBonds != RIGID_NONE)
3327    {
3328       NAMD_die("splitPatch hydrogen is required for rigidBonds");
3329    }
3330
3331    if (accelMDOn) {
3332        if(accelMDG){
3333            char msg[128];
3334            if(accelMDGiE < 1 || accelMDGiE > 2){
3335                sprintf(msg, "accelMDGiE was set to %d but it should be 1 or 2", accelMDGiE);
3336                NAMD_die(msg);
3337            }
3338            if(accelMDGRestart && accelMDGcMDSteps == 0)
3339                accelMDGcMDPrepSteps = 0;
3340            else if(accelMDGcMDSteps - accelMDGcMDPrepSteps < 2)
3341                NAMD_die("'accelMDGcMDSteps' should be larger than 'accelMDGcMDPrepSteps'");
3342
3343            if(accelMDGEquiSteps == 0)
3344                accelMDGEquiPrepSteps = 0;
3345            else if(accelMDGresetVaftercmd){
3346                if(accelMDGEquiPrepSteps <= 0)
3347                    NAMD_die("'accelMDGEquiPrepSteps' should be non-zero");
3348                if(accelMDGEquiSteps - accelMDGEquiPrepSteps < 1)
3349                    NAMD_die("'accelMDGEquiSteps' should be larger than 'accelMDGEquiPrepSteps'");
3350            }
3351
3352            //warn user that accelMD params will be ignored
3353            if(opts.defined("accelMDE"))
3354                iout << iWARN << "accelMDE will be ignored with accelMDG on.\n" << endi;
3355            if(opts.defined("accelMDalpha"))
3356                iout << iWARN << "accelMDalpha will be ignored with accelMDG on.\n" << endi;
3357            if(opts.defined("accelMDTE"))
3358                iout << iWARN << "accelMDTE will be ignored with accelMDG on.\n" << endi;
3359            if(opts.defined("accelMDTalpha"))
3360                iout << iWARN << "accelMDTalpha will be ignored with accelMDG on.\n" << endi;
3361        }
3362        else{
3363            if(!opts.defined("accelMDE") || !opts.defined("accelMDalpha"))
3364                NAMD_die("accelMDE and accelMDalpha are required for accelMD with accelMDG off");
3365
3366            if(accelMDdual && (!opts.defined("accelMDTE") || !opts.defined("accelMDTalpha"))){
3367                NAMD_die("accelMDTE and accelMDTalpha are required for accelMDdual with accelMDG off");
3368            }
3369        }
3370    }
3371
3372    //  Set the default value for the maximum movement parameter
3373    //  for minimization
3374    if (minimizeOn && (maximumMove == 0.0)) 
3375    {
3376       maximumMove = 0.75 * pairlistDist/stepsPerCycle;
3377    }
3378    if (adaptTempOn) {
3379      if (!adaptTempRescale && !adaptTempLangevin) 
3380         NAMD_die("Adaptive tempering needs to be coupled to either the Langevin thermostat or velocity rescaling.");
3381      if (opts.defined("adaptTempInFile") && (opts.defined("adaptTempTmin") ||
3382                                              opts.defined("adaptTempTmax") ||
3383                                              adaptTempBins != 0)) 
3384         NAMD_die("cannot simultaneously specify adaptTempInFile and any of {adaptTempTmin, adaptTempTmax,adaptTempBins} as these are read from the input file");
3385      if (!opts.defined("adaptTempInFile") && !(opts.defined("adaptTempTmin") &&
3386                                              opts.defined("adaptTempTmax") &&
3387                                              adaptTempBins != 0 ))  
3388         NAMD_die("Need to specify either adaptTempInFile or all of {adaptTempTmin, adaptTempTmax,adaptTempBins} if adaptTempMD is on.");
3389    }
3390
3391    langevinOnAtStartup = langevinOn;
3392    if (langevinOn) {
3393      if ( ! opts.defined("langevinDamping") ) langevinDamping = 0.0;
3394      if ( ! opts.defined("langevinHydrogen") ) langevinHydrogen = TRUE;
3395      if ( (opts.defined("langevinDamping") || opts.defined("langevinHydrogen"))
3396        && (opts.defined("langevinFile") || opts.defined("langevinCol")) )
3397        NAMD_die("To specify Langevin dynamics parameters, use either langevinDamping and langevinHydrogen or langevinFile and langevinCol.  Do not combine them.");
3398      if ( opts.defined("langevinHydrogen") && langevinDamping == 0.0 )
3399        NAMD_die("langevinHydrogen requires langevinDamping to be set.");
3400    }
3401    
3402    // BEGIN LA
3403    if (loweAndersenOn) {
3404        if (!opts.defined("loweAndersenRate")) loweAndersenRate = 100;
3405        if (!opts.defined("loweAndersenCutoff")) loweAndersenCutoff = 2.7;
3406    }
3407    // END LA
3408
3409    // BKR - stochastic velocity rescaling
3410    if (stochRescaleOn) {
3411      if (langevinOn || loweAndersenOn || tCoupleOn ||
3412          opts.defined("rescaleFreq") || opts.defined("reassignFreq"))
3413        NAMD_die("Stochastic velocity rescaling is incompatible with other temperature control methods");
3414      // This is largely the same default used in GROMACS.
3415      if (!opts.defined("stochRescaleFreq")) stochRescaleFreq = stepsPerCycle;
3416    }
3417
3418    if (opts.defined("rescaleFreq"))
3419    {
3420   if (!opts.defined("rescaleTemp"))
3421   {
3422     if (opts.defined("temperature"))
3423     {
3424       rescaleTemp = initialTemp;
3425     }
3426     else
3427     {
3428       NAMD_die("Must give a rescale temperature if rescaleFreq is defined");
3429     }
3430   }
3431    }
3432    else
3433    {
3434   rescaleFreq = -1;
3435   rescaleTemp = 0.0;
3436    }
3437
3438    if (opts.defined("rescaleTemp"))
3439    {
3440   if (!opts.defined("rescaleFreq"))
3441   {
3442     NAMD_die("Must give a rescale freqency if rescaleTemp is given");
3443   }
3444    }
3445
3446    if (opts.defined("reassignFreq"))
3447    {
3448   if (!opts.defined("reassignTemp"))
3449   {
3450     if (opts.defined("temperature"))
3451     {
3452       reassignTemp = initialTemp;
3453     }
3454     else
3455     {
3456       NAMD_die("Must give a reassign temperature if reassignFreq is defined");
3457     }
3458   }
3459    }
3460    else
3461    {
3462   reassignFreq = -1;
3463   reassignTemp = 0.0;
3464    }
3465
3466    if (opts.defined("reassignTemp"))
3467    {
3468   if (!opts.defined("reassignFreq"))
3469   {
3470     NAMD_die("Must give a reassignment freqency if reassignTemp is given");
3471   }
3472    }
3473
3474    if (opts.defined("reassignIncr"))
3475    {
3476   if (!opts.defined("reassignFreq"))