Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | File List | Namespace Members | Class Members | File Members

SimParameters.C

Go to the documentation of this file.
00001 
00007 /*****************************************************************************
00008  * $Source: /home/cvs/namd/cvsroot/namd2/src/SimParameters.C,v $
00009  * $Author: dhardy $
00010  * $Date: 2013/05/21 17:27:10 $
00011  * $Revision: 1.1416 $
00012  *****************************************************************************/
00013 
00020 #include "InfoStream.h"
00021 #include "ComputeNonbondedUtil.h"
00022 #include "ConfigList.h"
00023 #include "SimParameters.h"
00024 #include "ParseOptions.h"
00025 #include "structures.h"
00026 #include "Communicate.h"
00027 #include "MStream.h"
00028 #include "Output.h"
00029 #include <stdio.h>
00030 #include <time.h>
00031 #ifdef NAMD_FFTW
00032 #ifdef NAMD_FFTW_3
00033 #include <fftw3.h>
00034 #else
00035 // fftw2 doesn't have these defined
00036 #define fftwf_malloc fftw_malloc
00037 #define fftwf_free fftw_free
00038 #ifdef NAMD_FFTW_NO_TYPE_PREFIX
00039 #include <fftw.h>
00040 #include <rfftw.h>
00041 #else
00042 #include <sfftw.h>
00043 #include <srfftw.h>
00044 #endif
00045 #endif
00046 #endif
00047 #if defined(WIN32) && !defined(__CYGWIN__)
00048 #include <direct.h>
00049 #define CHDIR _chdir
00050 #define MKDIR(X) mkdir(X)
00051 #define PATHSEP '\\'
00052 #define PATHSEPSTR "\\"
00053 #else
00054 #include <unistd.h>
00055 #define CHDIR chdir
00056 #define MKDIR(X) mkdir(X,0777)
00057 #define PATHSEP '/'
00058 #define PATHSEPSTR "/"
00059 #endif
00060 #include <fcntl.h>
00061 #include <sys/stat.h>
00062 #ifdef WIN32
00063 #include <io.h>
00064 #define access(PATH,MODE) _access(PATH,00)
00065 #endif
00066 #include <fstream>
00067 using namespace std;
00068 
00069 #ifdef WIN32
00070 extern "C" {
00071   double erfc(double);
00072 }
00073 #endif
00074 
00075 #include "strlib.h"    //  For strcasecmp and strncasecmp
00076 
00077 //#define DEBUGM
00078 #include "Debug.h"
00079 
00080 #define XXXBIGREAL 1.0e32
00081 
00082 /************************************************************************/
00083 /*                  */
00084 /*      FUNCTION initialize_config_data      */
00085 /*                  */
00086 /*  This function is used by the master process to populate the     */
00087 /*   simulation parameters from a ConfigList object that is passed to   */
00088 /*   it.  Each parameter is checked to make sure that it has a value    */
00089 /*   that makes sense, and that it doesn't conflict with any other      */
00090 /*   values that have been given.          */
00091 /*                  */
00092 /************************************************************************/
00093 
00094 void SimParameters::initialize_config_data(ConfigList *config, char *&cwd)
00095 
00096 {
00097 
00098    ParseOptions opts;   //  Object to check consistency of config file
00099 
00100    config_parser(opts);
00101 
00103    if (!opts.check_consistency()) 
00104    {
00105       NAMD_die("Internal error in configuration file parser");
00106    }
00107 
00108    // Now, feed the object with the actual configuration options through the
00109    // ParseOptions file and make sure everything is OK
00110    if (!opts.set(*config)) 
00111    {
00112       NAMD_die("ERROR(S) IN THE CONFIGURATION FILE");
00113    }
00114 
00116 
00117    check_config(opts,config,cwd);
00118 
00119    print_config(opts,config,cwd);
00120 
00121 }
00122 
00123 /************************************************************************/
00124 /*                                                                      */
00125 /*      FUNCTION scriptSet                                              */
00126 /*                                                                      */
00127 /************************************************************************/
00128 
00129 int atobool(const char *s) {
00130   return ( (! strncasecmp(s,"yes",8)) ||
00131            (! strncasecmp(s,"on",8)) || 
00132            (! strncasecmp(s,"true",8)) );
00133 };
00134                          
00135 void SimParameters::scriptSet(const char *param, const char *value) {
00136 
00137   if ( CkMyRank() ) return;
00138 
00139 #define MAX_SCRIPT_PARAM_SIZE 128
00140 #define SCRIPT_PARSE_BOOL(NAME,VAR) { if ( ! strncasecmp(param,(NAME),MAX_SCRIPT_PARAM_SIZE) ) { (VAR) = atobool(value); return; } }
00141 #define SCRIPT_PARSE_INT(NAME,VAR) { if ( ! strncasecmp(param,(NAME),MAX_SCRIPT_PARAM_SIZE) ) { (VAR) = atoi(value); return; } }
00142 #define SCRIPT_PARSE_FLOAT(NAME,VAR) { if ( ! strncasecmp(param,(NAME),MAX_SCRIPT_PARAM_SIZE) ) { (VAR) = atof(value); return; } }
00143 #define SCRIPT_PARSE_MOD_FLOAT(NAME,VAR,MOD) { if ( ! strncasecmp(param,(NAME),MAX_SCRIPT_PARAM_SIZE) ) { (VAR) = atof(value) MOD; return; } }
00144 #define SCRIPT_PARSE_VECTOR(NAME,VAR) { if ( ! strncasecmp(param,(NAME),MAX_SCRIPT_PARAM_SIZE) ) { (VAR).set(value); return; } }
00145 #define SCRIPT_PARSE_STRING(NAME,VAR) { if ( ! strncasecmp(param,(NAME),MAX_SCRIPT_PARAM_SIZE) ) { strcpy(VAR,value); return; } }
00146 
00147   SCRIPT_PARSE_FLOAT("scriptArg1",scriptArg1)
00148   SCRIPT_PARSE_FLOAT("scriptArg2",scriptArg2)
00149   SCRIPT_PARSE_FLOAT("scriptArg3",scriptArg3)
00150   SCRIPT_PARSE_FLOAT("scriptArg4",scriptArg4)
00151   SCRIPT_PARSE_FLOAT("scriptArg5",scriptArg5)
00152   SCRIPT_PARSE_INT("numsteps",N)
00153   SCRIPT_PARSE_INT("firsttimestep",firstTimestep)
00154   SCRIPT_PARSE_FLOAT("reassignTemp",reassignTemp)
00155   SCRIPT_PARSE_FLOAT("rescaleTemp",rescaleTemp)
00156   SCRIPT_PARSE_BOOL("velocityQuenching",minimizeOn)
00157   SCRIPT_PARSE_BOOL("maximumMove",maximumMove)
00158   // SCRIPT_PARSE_BOOL("Langevin",langevinOn)
00159   SCRIPT_PARSE_FLOAT("langevinTemp",langevinTemp)
00160   SCRIPT_PARSE_BOOL("langevinBAOAB",langevin_useBAOAB) // [!!] Use the BAOAB integrator or not
00161   SCRIPT_PARSE_FLOAT("loweAndersenTemp",loweAndersenTemp) // BEGIN LA, END LA
00162   SCRIPT_PARSE_FLOAT("initialTemp",initialTemp)
00163   SCRIPT_PARSE_BOOL("useGroupPressure",useGroupPressure)
00164   SCRIPT_PARSE_BOOL("useFlexibleCell",useFlexibleCell)
00165   SCRIPT_PARSE_BOOL("useConstantArea",useConstantArea)
00166   SCRIPT_PARSE_BOOL("fixCellDims",fixCellDims)
00167   SCRIPT_PARSE_BOOL("fixCellDimX",fixCellDimX)
00168   SCRIPT_PARSE_BOOL("fixCellDimY",fixCellDimY)
00169   SCRIPT_PARSE_BOOL("fixCellDimZ",fixCellDimZ)
00170   SCRIPT_PARSE_BOOL("useConstantRatio",useConstantRatio)
00171   SCRIPT_PARSE_BOOL("LangevinPiston",langevinPistonOn)
00172   SCRIPT_PARSE_MOD_FLOAT("LangevinPistonTarget",
00173                         langevinPistonTarget,/PRESSUREFACTOR)
00174   SCRIPT_PARSE_FLOAT("LangevinPistonPeriod",langevinPistonPeriod)
00175   SCRIPT_PARSE_FLOAT("LangevinPistonDecay",langevinPistonDecay)
00176   SCRIPT_PARSE_FLOAT("LangevinPistonTemp",langevinPistonTemp)
00177   SCRIPT_PARSE_MOD_FLOAT("SurfaceTensionTarget",
00178                         surfaceTensionTarget,*(100.0/PRESSUREFACTOR))
00179   SCRIPT_PARSE_BOOL("BerendsenPressure",berendsenPressureOn)
00180   SCRIPT_PARSE_MOD_FLOAT("BerendsenPressureTarget",
00181                         berendsenPressureTarget,/PRESSUREFACTOR)
00182   SCRIPT_PARSE_MOD_FLOAT("BerendsenPressureCompressibility",
00183                         berendsenPressureCompressibility,*PRESSUREFACTOR)
00184   SCRIPT_PARSE_FLOAT("BerendsenPressureRelaxationTime",
00185                                 berendsenPressureRelaxationTime)
00186   SCRIPT_PARSE_FLOAT("constraintScaling",constraintScaling)
00187   SCRIPT_PARSE_FLOAT("consForceScaling",consForceScaling)
00188   SCRIPT_PARSE_STRING("outputname",outputFilename)
00189   SCRIPT_PARSE_STRING("tclBCArgs",tclBCArgs)
00190   SCRIPT_PARSE_VECTOR("eField",eField)
00191   SCRIPT_PARSE_FLOAT("eFieldFreq",eFieldFreq)
00192   SCRIPT_PARSE_FLOAT("eFieldPhase",eFieldPhase) 
00193   SCRIPT_PARSE_FLOAT("accelMDE",accelMDE) 
00194   SCRIPT_PARSE_FLOAT("accelMDalpha",accelMDalpha) 
00195   SCRIPT_PARSE_FLOAT("accelMDTE",accelMDTE) 
00196   SCRIPT_PARSE_FLOAT("accelMDTalpha",accelMDTalpha) 
00197   SCRIPT_PARSE_VECTOR("stirAxis",stirAxis)
00198   SCRIPT_PARSE_VECTOR("stirPivot",stirPivot)
00199   if ( ! strncasecmp(param,"mgridforcescale",MAX_SCRIPT_PARAM_SIZE) ) {
00200     NAMD_die("Can't yet modify mgridforcescale in a script");
00201     return;
00202   }
00203   if ( ! strncasecmp(param,"mgridforcevoff",MAX_SCRIPT_PARAM_SIZE) ) {
00204     NAMD_die("Can't yet modify mgridforcevoff in a script");
00205     return;
00206   }
00207    
00208   if ( ! strncasecmp(param,"fixedatoms",MAX_SCRIPT_PARAM_SIZE) ) {
00209     if ( ! fixedAtomsOn )
00210       NAMD_die("FixedAtoms may not be enabled in a script.");
00211     if ( ! fixedAtomsForces )
00212       NAMD_die("To use fixedAtoms in script first use fixedAtomsForces yes.");
00213     fixedAtomsOn = atobool(value);
00214     return;
00215   }
00216 
00217 //Modifications for alchemical fep
00218   SCRIPT_PARSE_INT("alchEquilSteps",alchEquilSteps)
00219 
00220   if ( ! strncasecmp(param,"alchRepLambda",MAX_SCRIPT_PARAM_SIZE) ) {
00221     alchRepLambda = atof(value);
00222     alchFepWCARepuOn = true;
00223     alchFepWCADispOn = false;
00224     alchFepElecOn    = false;
00225     ComputeNonbondedUtil::select();
00226     return;
00227   }
00228 
00229   if ( ! strncasecmp(param,"alchDispLambda",MAX_SCRIPT_PARAM_SIZE) ) {
00230     alchDispLambda = atof(value);
00231     alchFepWCARepuOn = false;
00232     alchFepWCADispOn = true;
00233     alchFepElecOn    = false;
00234     ComputeNonbondedUtil::select();
00235     return;
00236   }
00237 
00238   if ( ! strncasecmp(param,"alchElecLambda",MAX_SCRIPT_PARAM_SIZE) ) {
00239     alchElecLambda = atof(value);
00240     alchFepWCARepuOn = false;
00241     alchFepWCADispOn = false;
00242     alchFepElecOn    = true;
00243     ComputeNonbondedUtil::select();
00244     return;
00245   }
00246 
00247 
00248   if ( ! strncasecmp(param,"alchFepWCArcut1",MAX_SCRIPT_PARAM_SIZE) ) {
00249     alchFepWCArcut1 = atof(value);
00250     ComputeNonbondedUtil::select();
00251     return;
00252   }
00253 
00254   if ( ! strncasecmp(param,"alchFepWCArcut2",MAX_SCRIPT_PARAM_SIZE) ) {
00255     alchFepWCArcut2 = atof(value);
00256     ComputeNonbondedUtil::select();
00257     return;
00258   }
00259 
00260   if ( ! strncasecmp(param,"alchFepWCArcut3",MAX_SCRIPT_PARAM_SIZE) ) {
00261     alchFepWCArcut3 = atof(value);
00262     ComputeNonbondedUtil::select();
00263     return;
00264   }
00265 
00266   if ( ! strncasecmp(param,"alchLambda",MAX_SCRIPT_PARAM_SIZE) ) {
00267     alchLambda = atof(value);
00268     ComputeNonbondedUtil::select();
00269     return;
00270   }
00271 
00272   if ( ! strncasecmp(param,"alchLambda2",MAX_SCRIPT_PARAM_SIZE) ) {
00273     alchLambda2 = atof(value);
00274     ComputeNonbondedUtil::select();
00275     return;
00276   }
00277 //fepe
00278 
00279 // REDUNDANT TI BEGINS
00280 //  if ( ! strncasecmp(param,"tiLambda",MAX_SCRIPT_PARAM_SIZE) ) {
00281 //    alchLambda = atof(value);
00282 //    ComputeNonbondedUtil::select();
00283 //    return;
00284 //  }
00285 // REDUNDANT TI ENDS
00286 
00287   if ( ! strncasecmp(param,"nonbondedScaling",MAX_SCRIPT_PARAM_SIZE) ) {
00288     nonbondedScaling = atof(value);
00289     ComputeNonbondedUtil::select();
00290     return;
00291   }
00292   if ( ! strncasecmp(param,"commOnly",MAX_SCRIPT_PARAM_SIZE) ) {
00293     commOnly = atobool(value);
00294     ComputeNonbondedUtil::select();
00295     return;
00296   }
00297 
00298   char *error = new char[2 * MAX_SCRIPT_PARAM_SIZE + 100];
00299   sprintf(error,"Setting parameter %s from script failed!\n",param);
00300   NAMD_die(error);
00301 
00302 }
00303 
00304 /************************************************************************/
00305 /*                                                                      */
00306 /*      FUNCTION config_parser                                          */
00307 /*                                                                      */
00308 /************************************************************************/
00309                          
00310 void SimParameters::config_parser(ParseOptions &opts) {
00311 
00312    //  Set all variable to fallback default values.  This is not really
00313    //  necessary, as we give default values when we set up the ParseOptions
00314    //  object, but it helps the debuggers figure out we've initialized the
00315    //  variables.
00316    HydrogenBonds = FALSE;
00317    useAntecedent = TRUE;
00318    aaAngleExp = 2;
00319    haAngleExp = 4;
00320    distAttExp = 4;
00321    distRepExp = 6;
00322    dhaCutoffAngle = 100.0;
00323    dhaOnAngle = 60.0;
00324    dhaOffAngle = 80.0;
00325    daCutoffDist = 7.5;
00326    daOnDist = 5.5;
00327    daOffDist = 6.5;
00328 
00329    config_parser_basic(opts);
00330    config_parser_fileio(opts);
00331    config_parser_fullelect(opts);
00332    config_parser_methods(opts);
00333    config_parser_constraints(opts);
00334    #ifdef OPENATOM_VERSION
00335    config_parser_openatom(opts);
00336    #endif // OPENATOM_VERSION
00337 
00338    config_parser_gridforce(opts);
00339    config_parser_mgridforce(opts);
00340    config_parser_movdrag(opts);
00341    config_parser_rotdrag(opts);
00342    config_parser_constorque(opts);
00343    config_parser_boundary(opts);
00344    config_parser_misc(opts);
00345 
00346 }
00347 
00348 void SimParameters::config_parser_basic(ParseOptions &opts) {
00349    
00350    //  So first we set up the ParseOptions objects so that it will check
00351    //  all of the logical rules that the configuration file must follow.
00352 
00353    opts.optional("main", "obsolete", "used to flag obsolete options",
00354     PARSE_STRING);
00355 
00357    opts.require("main", "timestep", "size of the timestep, in fs",
00358     &dt, 1.0);
00359    opts.range("timestep", NOT_NEGATIVE);
00360    opts.units("timestep", N_FSEC);
00361 
00362    opts.optional("main", "numsteps", "number of timesteps to perform",
00363     &N,0);
00364    opts.range("numsteps", NOT_NEGATIVE);
00365 
00366    opts.optional("main", "stepspercycle",
00367       "Number of steps between atom migrations", 
00368       &stepsPerCycle, 20);
00369    opts.range("stepspercycle", POSITIVE);
00370 
00371    opts.require("main", "cutoff", "local electrostatic and Vdw distance", 
00372       &cutoff);
00373    opts.range("cutoff", POSITIVE);
00374    opts.units("cutoff", N_ANGSTROM);
00375    
00376    opts.optional("main", "nonbondedScaling", "nonbonded scaling factor",
00377      &nonbondedScaling, 1.0);
00378    opts.range("nonbondedScaling", NOT_NEGATIVE);
00379 
00380    opts.optional("main", "limitDist", "limit nonbonded below this distance",
00381      &limitDist, 0.0);
00382    opts.range("limitDist", NOT_NEGATIVE);
00383 
00384    opts.require("main", "exclude", "Electrostatic and VDW exclusion policy",
00385     PARSE_STRING);
00386 
00387    opts.optional("exclude", "1-4scaling", "1-4 electrostatic scaling factor",
00388      &scale14, 1.0);
00389    opts.range("1-4scaling", POSITIVE);
00390 
00391    opts.optionalB("main", "switching",
00392      "Should a smoothing function be used?", &switchingActive, TRUE);
00393    
00394    opts.optionalB("switching", "vdwForceSwitching",
00395      "Use force switching for vdw?", &vdwForceSwitching, FALSE);
00396    
00397    opts.optional("switching", "switchdist",
00398      "Distance for switching function activation",
00399      &switchingDist);
00400    opts.range("switchdist", POSITIVE);
00401    opts.units("switchdist", N_ANGSTROM);
00402    
00403    opts.optionalB("main", "martiniSwitching",
00404      "Use Martini residue-based coarse-grain switching?", &martiniSwitching, FALSE);
00405    opts.optionalB("main", "martiniDielAllow",
00406      "Allow use of dielectric != 15.0 when using Martini", &martiniDielAllow, FALSE);
00407 
00408    opts.optional("main", "pairlistdist",  "Pairlist inclusion distance",
00409      &pairlistDist);
00410    opts.range("pairlistdist", POSITIVE);
00411    opts.units("pairlistdist", N_ANGSTROM);
00412 
00413    opts.optional("main", "pairlistMinProcs",  "Min procs for pairlists",
00414      &pairlistMinProcs,1);
00415    opts.range("pairlistMinProcs", POSITIVE);
00416 
00417    opts.optional("main", "pairlistsPerCycle",  "regenerate x times per cycle",
00418      &pairlistsPerCycle,2);
00419    opts.range("pairlistsPerCycle", POSITIVE);
00420 
00421    opts.optional("main", "outputPairlists", "how often to print warnings",
00422      &outputPairlists, 0);
00423    opts.range("outputPairlists", NOT_NEGATIVE);
00424 
00425    opts.optional("main", "pairlistShrink",  "tol *= (1 - x) on regeneration",
00426      &pairlistShrink,0.01);
00427    opts.range("pairlistShrink", NOT_NEGATIVE);
00428 
00429    opts.optional("main", "pairlistGrow",  "tol *= (1 + x) on trigger",
00430      &pairlistGrow, 0.01);
00431    opts.range("pairlistGrow", NOT_NEGATIVE);
00432 
00433    opts.optional("main", "pairlistTrigger",  "trigger is atom > (1 - x) * tol",
00434      &pairlistTrigger, 0.3);
00435    opts.range("pairlistTrigger", NOT_NEGATIVE);
00436 
00437    opts.optional("main", "temperature", "initial temperature",
00438      &initialTemp);
00439    opts.range("temperature", NOT_NEGATIVE);
00440    opts.units("temperature", N_KELVIN);
00441 
00442    opts.optionalB("main", "COMmotion", "allow initial center of mass movement",
00443       &comMove, FALSE);
00444 
00445    opts.optionalB("main", "zeroMomentum", "constrain center of mass",
00446       &zeroMomentum, FALSE);
00447    opts.optionalB("zeroMomentum", "zeroMomentumAlt", "constrain center of mass",
00448       &zeroMomentumAlt, FALSE);
00449 
00450    opts.optionalB("main", "wrapWater", "wrap waters around periodic boundaries on output",
00451       &wrapWater, FALSE);
00452    opts.optionalB("main", "wrapAll", "wrap all clusters around periodic boundaries on output",
00453       &wrapAll, FALSE);
00454    opts.optionalB("main", "wrapNearest", "wrap to nearest image to cell origin",
00455       &wrapNearest, FALSE);
00456 
00457    opts.optional("main", "dielectric", "dielectric constant",
00458      &dielectric, 1.0);
00459    opts.range("dielectric", POSITIVE); // Hmmm, dielectric < 1 ...
00460 
00461    opts.optional("main", "margin", "Patch width margin", &margin, XXXBIGREAL);
00462    opts.range("margin", NOT_NEGATIVE);
00463    opts.units("margin", N_ANGSTROM);
00464 
00465    opts.optional("main", "seed", "Initial random number seed", &randomSeed);
00466    opts.range("seed", POSITIVE);
00467 
00468    opts.optional("main", "outputEnergies", "How often to print energies in timesteps",
00469      &outputEnergies, 1);
00470    opts.range("outputEnergies", POSITIVE);
00471      
00472    opts.optional("main", "outputMomenta", "How often to print linear and angular momenta in timesteps",
00473      &outputMomenta, 0);
00474    opts.range("outputMomenta", NOT_NEGATIVE);
00475      
00476    opts.optional("main", "outputTiming", "How often to print timing data in timesteps",
00477      &outputTiming);
00478    opts.range("outputTiming", NOT_NEGATIVE);
00479      
00480    opts.optional("main", "outputCudaTiming", "How often to print CUDA timing data in timesteps",
00481      &outputCudaTiming, 0);
00482    opts.range("outputCudaTiming", NOT_NEGATIVE);
00483      
00484    opts.optional("main", "outputPressure", "How often to print pressure data in timesteps",
00485      &outputPressure, 0);
00486    opts.range("outputPressure", NOT_NEGATIVE);
00487      
00488    opts.optionalB("main", "mergeCrossterms", "merge crossterm energy with dihedral when printing?",
00489       &mergeCrossterms, TRUE);
00490 
00491    opts.optional("main", "MTSAlgorithm", "Multiple timestep algorithm",
00492     PARSE_STRING);
00493 
00494    opts.optional("main", "longSplitting", "Long range force splitting option",
00495     PARSE_STRING);
00496 
00497    opts.optionalB("main", "ignoreMass", "Do not use masses to find hydrogen atoms",
00498     &ignoreMass, FALSE);
00499 
00500    opts.optional("main", "splitPatch", "Atom into patch splitting option",
00501     PARSE_STRING);
00502    opts.optional("main", "hgroupCutoff", "Hydrogen margin", &hgroupCutoff, 2.5);
00503 
00504    opts.optional("main", "extendedSystem",
00505     "Initial configuration of extended system variables and periodic cell",
00506     PARSE_STRING);
00507 
00508    opts.optional("main", "cellBasisVector1", "Basis vector for periodic cell",
00509     &cellBasisVector1);
00510    opts.optional("main", "cellBasisVector2", "Basis vector for periodic cell",
00511     &cellBasisVector2);
00512    opts.optional("main", "cellBasisVector3", "Basis vector for periodic cell",
00513     &cellBasisVector3);
00514    opts.optional("main", "cellOrigin", "Fixed center of periodic cell",
00515     &cellOrigin);
00516 
00517    opts.optionalB("main", "molly", "Rigid bonds to hydrogen",&mollyOn,FALSE);
00518    opts.optional("main", "mollyTolerance", "Error tolerance for MOLLY",
00519                  &mollyTol, 0.00001);
00520    opts.optional("main", "mollyIterations", 
00521                  "Max number of iterations for MOLLY", &mollyIter, 100);
00522 
00523    opts.optional("main", "rigidBonds", "Rigid bonds to hydrogen",PARSE_STRING);
00524    opts.optional("main", "rigidTolerance", 
00525                  "Error tolerance for rigid bonds to hydrogen",
00526                  &rigidTol, 1.0e-8);
00527    opts.optional("main", "rigidIterations", 
00528                  "Max number of SHAKE iterations for rigid bonds to hydrogen",
00529                  &rigidIter, 100);
00530    opts.optionalB("main", "rigidDieOnError", 
00531                  "Die if rigidTolerance is not achieved after rigidIterations",
00532                  &rigidDie, TRUE);
00533    opts.optionalB("main", "useSettle",
00534                   "Use the SETTLE algorithm for rigid waters",
00535                  &useSettle, TRUE);
00536 
00537    opts.optional("main", "nonbondedFreq", "Nonbonded evaluation frequency",
00538     &nonbondedFrequency, 1);
00539    opts.range("nonbondedFreq", POSITIVE);
00540 
00541    opts.optionalB("main", "outputPatchDetails", "print number of atoms in each patch",
00542       &outputPatchDetails, FALSE);
00543    opts.optionalB("main", "staticAtomAssignment", "never migrate atoms",
00544       &staticAtomAssignment, FALSE);
00545    opts.optionalB("main", "lonePairs", "Enable lone pairs", &lonepairs, FALSE);
00546    opts.optional("main", "waterModel", "Water model to use", PARSE_STRING);
00547    opts.optionalB("main", "LJcorrection", "Apply analytical tail corrections for energy and virial", &LJcorrection, FALSE);
00548 }
00549 
00550 void SimParameters::config_parser_fileio(ParseOptions &opts) {
00551    
00553 
00554    opts.optional("main", "cwd", "current working directory", PARSE_STRING);
00555 
00556 // In order to include AMBER options, "coordinates", "structure"
00557 // and "parameters" are now optional, not required. The presence
00558 // of them will be checked later in check_config()
00559 
00560 //   opts.require("main", "coordinates", "initial PDB coordinate file",
00561 //    PARSE_STRING);
00562    opts.optional("main", "coordinates", "initial PDB coordinate file",
00563     PARSE_STRING);
00564 
00565    opts.optional("main", "velocities",
00566      "initial velocities, given as a PDB file", PARSE_STRING);
00567    opts.optional("main", "binvelocities",
00568      "initial velocities, given as a binary restart", PARSE_STRING);
00569    opts.optional("main", "bincoordinates",
00570      "initial coordinates in a binary restart file", PARSE_STRING);
00571 #ifdef MEM_OPT_VERSION
00572    opts.optional("main", "binrefcoords",
00573      "reference coordinates in a binary restart file", PARSE_STRING);
00574 #endif
00575 
00576 //   opts.require("main", "structure", "initial PSF structure file",
00577 //    PARSE_STRING);
00578    opts.optional("main", "structure", "initial PSF structure file",
00579     PARSE_STRING);   
00580 
00581 //   opts.require("main", "parameters",
00582 //"CHARMm 19 or CHARMm 22 compatable force field file (multiple "
00583 //"inputs allowed)", PARSE_MULTIPLES);
00584    opts.optional("main", "parameters",
00585 "CHARMm 19 or CHARMm 22 compatable force field file (multiple "
00586 "inputs allowed)", PARSE_MULTIPLES);
00587 
00588 
00589    //****** BEGIN CHARMM/XPLOR type changes
00591    opts.optionalB("parameters", "paraTypeXplor", "Parameter file in Xplor format?", &paraTypeXplorOn, FALSE);
00592    opts.optionalB("parameters", "paraTypeCharmm", "Parameter file in Charmm format?", &paraTypeCharmmOn, FALSE); 
00593    //****** END CHARMM/XPLOR type changes
00594 
00595    // Ported by JLai -- JE - Go parameters
00596    opts.optionalB("main", "GromacsPair", "Separately calculate pair interactions", &goGroPair, FALSE);
00597    opts.optionalB("main", "GoForcesOn", "Go forces will be calculated", &goForcesOn, FALSE);
00598    opts.require("GoForcesOn", "GoParameters", "Go parameter file", goParameters);
00599    opts.require("GoForcesOn", "GoCoordinates", "target coordinates for Go forces", goCoordinates);
00600    // Added by JLai -- Go-method parameter -- switches between using the matrix and sparse matrix representations [6.3.11] 
00601    //   opts.optional("GoForcesOn", "GoMethod", "Which type of matrix should be used to store Go contacts?", &goMethod);
00602    // Added by JLai -- Go-method parameter [8.2.11]
00603    //opts.optional("GoForcesOn", "GoMethod", "Which type of matrix should be used to store Go contacts?", PARSE_STRING); 
00604    opts.require("GoForcesOn", "GoMethod", "Which type of matrix should be used to store Go contacts?", PARSE_STRING); 
00605   // End of Port -- JL
00606    
00607    opts.require("main", "outputname",
00608     "prefix for the final PDB position and velocity filenames", 
00609     outputFilename);
00610 
00611    opts.optional("main", "auxFile", "Filename for data stream output",
00612      auxFilename);
00613 
00614    opts.optional("main", "numinputprocs", "Number of pes to use for parallel input", 
00615                  &numinputprocs, 0);
00616    opts.range("numinputprocs", NOT_NEGATIVE);
00617 
00618    opts.optional("main", "numoutputprocs", "Number of pes to use for parallel output", 
00619                  &numoutputprocs, 0);
00620    opts.range("numoutputprocs", NOT_NEGATIVE);
00621    opts.optional("main", "numoutputwriters", "Number of output processors that simultaneously write to an output file", 
00622                  &numoutputwrts, 1);
00623    opts.range("numoutputwriters", NOT_NEGATIVE);
00624 
00625    opts.optional("main", "DCDfreq", "Frequency of DCD trajectory output, in "
00626     "timesteps", &dcdFrequency, 0);
00627    opts.range("DCDfreq", NOT_NEGATIVE);
00628    opts.optional("DCDfreq", "DCDfile", "DCD trajectory output file name",
00629      dcdFilename);
00630    opts.optionalB("DCDfreq", "DCDunitcell", "Store unit cell in dcd timesteps?",
00631        &dcdUnitCell);
00632 
00633    opts.optional("main", "velDCDfreq", "Frequency of velocity "
00634     "DCD output, in timesteps", &velDcdFrequency, 0);
00635    opts.range("velDCDfreq", NOT_NEGATIVE);
00636    opts.optional("velDCDfreq", "velDCDfile", "velocity DCD output file name",
00637      velDcdFilename);
00638    
00639    opts.optional("main", "forceDCDfreq", "Frequency of force"
00640     "DCD output, in timesteps", &forceDcdFrequency, 0);
00641    opts.range("forceDCDfreq", NOT_NEGATIVE);
00642    opts.optional("forceDCDfreq", "forceDCDfile", "force DCD output file name",
00643      forceDcdFilename);
00644    
00645    opts.optional("main", "XSTfreq", "Frequency of XST trajectory output, in "
00646     "timesteps", &xstFrequency, 0);
00647    opts.range("XSTfreq", NOT_NEGATIVE);
00648    opts.optional("XSTfreq", "XSTfile", "Extended sytem trajectory output "
00649     "file name", xstFilename);
00650 
00651    opts.optional("main", "restartfreq", "Frequency of restart file "
00652     "generation", &restartFrequency, 0);
00653    opts.range("restartfreq", NOT_NEGATIVE);
00654    opts.optional("restartfreq", "restartname", "Prefix for the position and "
00655      "velocity PDB files used for restarting", restartFilename);
00656    opts.optionalB("restartfreq", "restartsave", "Save restart files with "
00657      "unique filenames rather than overwriting", &restartSave, FALSE);
00658 
00659    opts.optionalB("restartfreq", "binaryrestart", "Specify use of binary restart files ", 
00660        &binaryRestart, TRUE);
00661 
00662    opts.optionalB("outputname", "binaryoutput", "Specify use of binary output files ", 
00663        &binaryOutput, TRUE);
00664 
00665    opts.optionalB("main", "amber", "Is it AMBER force field?",
00666        &amberOn, FALSE);
00667    opts.optionalB("amber", "readexclusions", "Read exclusions from parm file?",
00668        &readExclusions, TRUE);
00669    opts.require("amber", "scnb", "1-4 VDW interactions are divided by scnb",
00670        &vdwscale14, 2.0);
00671    opts.require("amber", "parmfile", "AMBER parm file", PARSE_STRING);
00672    opts.optional("amber", "ambercoor", "AMBER coordinate file", PARSE_STRING);
00673 
00674 //Modifications for alchemical fep
00675 // begin fep output options    
00676    opts.optional("alch", "alchoutfreq", "Frequency of alchemical energy"
00677      "output in timesteps", &alchOutFreq, 5);
00678    opts.range("alchoutfreq", NOT_NEGATIVE);
00679    opts.optional("alchoutfreq", "alchoutfile", "Alchemical energy output"
00680        "filename", alchOutFile);
00681 // end fep output options
00682 //fepe
00683 
00684 // REDUNDANT TI BEGINS 
00685 //   opts.optional("thermInt", "tioutfreq", "Frequency of TI energy output in "
00686 //     "timesteps", &tiOutFreq, 5);
00687 //   opts.range("tioutfreq", NOT_NEGATIVE);
00688 //   opts.optional("tioutfreq", "tioutfile", "TI energy output filename",
00689 //     tiOutFile);
00690 // REDUNDANT TI ENDS
00691 
00692    /* GROMACS options */
00693    opts.optionalB("main", "gromacs", "Use GROMACS-like force field?",
00694        &gromacsOn, FALSE);
00695    opts.require("gromacs", "grotopfile", "GROMACS topology file",
00696                 PARSE_STRING);
00697    opts.optional("gromacs", "grocoorfile","GROMACS coordinate file",
00698                  PARSE_STRING);
00699 
00700   // OPLS options
00701    opts.optionalB("main", "vdwGeometricSigma",
00702        "Use geometric mean to combine L-J sigmas, as for OPLS",
00703        &vdwGeometricSigma, FALSE);
00704 
00705    // load/store computeMap
00706    opts.optional("main", "computeMapFile", "Filename for computeMap",
00707      computeMapFilename);
00708    opts.optionalB("main", "storeComputeMap", "store computeMap?",
00709        &storeComputeMap, FALSE);
00710    opts.optionalB("main", "loadComputeMap", "load computeMap?",
00711        &loadComputeMap, FALSE);
00712 }
00713 
00714 
00715 void SimParameters::config_parser_fullelect(ParseOptions &opts) {
00716    
00718 #ifdef DPMTA
00719    DebugM(1,"DPMTA setup start\n");
00720    //  PMTA is included, so really get these values
00721    opts.optionalB("main", "FMA", "Should FMA be used?", &FMAOn, FALSE);
00722    opts.optional("FMA", "FMALevels", "Tree levels to use in FMA", &FMALevels,
00723      5);
00724    opts.range("FMALevels", POSITIVE);
00725    opts.optional("FMA", "FMAMp", "Number of FMA multipoles", &FMAMp, 8);
00726    opts.range("FMAMp", POSITIVE);
00727    opts.optionalB("FMA", "FMAFFT", "Use FFT enhancement in FMA?", &FMAFFTOn, TRUE);
00728    opts.optional("FMAFFT", "FMAFFTBlock", "FFT blocking factor",
00729     &FMAFFTBlock, 4);
00730    opts.range("FMAFFTBlock", POSITIVE);
00731    DebugM(1,"DPMTA setup end\n");
00732 #else
00733    //  PMTA is NOT included.  So just set all the values to 0.
00734    FMAOn = FALSE;
00735    FMALevels = 0;
00736    FMAMp = 0;
00737    FMAFFTOn = FALSE;
00738    FMAFFTBlock = 0;
00739 #endif
00740 
00741    opts.optional("main", "fullElectFrequency",
00742       "Number of steps between full electrostatic executions", 
00743       &fullElectFrequency);
00744    opts.range("fullElectFrequency", POSITIVE);
00745 
00746    //  USE OF THIS PARAMETER DISCOURAGED
00747    opts.optional("main", "fmaFrequency",
00748       "Number of steps between full electrostatic executions", 
00749       &fmaFrequency);
00750    opts.range("fmaFrequency", POSITIVE);
00751 
00752    opts.optional("main", "fmaTheta",
00753       "FMA theta parameter value", 
00754       &fmaTheta,0.715);
00755    opts.range("fmaTheta", POSITIVE);
00756 
00757    opts.optionalB("main", "FullDirect", "Should direct calculations of full electrostatics be performed?",
00758       &fullDirectOn, FALSE);
00759 
00760 
00762 
00763    opts.optionalB("main", "MSM",
00764        "Use multilevel summation method for electrostatics?",
00765        &MSMOn, FALSE);
00766    opts.optional("MSM", "MSMQuality", "MSM quality",
00767        &MSMQuality, 0);
00768    opts.optional("MSM", "MSMApprox", "MSM approximation",
00769        &MSMApprox, 0);
00770    opts.optional("MSM", "MSMSplit", "MSM splitting",
00771        &MSMSplit, 0);
00772    opts.optional("MSM", "MSMLevels", "MSM maximum number of levels",
00773        &MSMLevels, 0);  // set to 0 adapts to as many as needed
00774    opts.optional("MSM", "MSMGridSpacing", "MSM grid spacing (Angstroms)",
00775        &MSMGridSpacing, 2.5);
00776    opts.optional("MSM", "MSMPadding", "MSM padding (Angstroms)",
00777        &MSMPadding, 2.5);
00778    opts.optional("MSM", "MSMBlockSizeX",
00779        "MSM grid block size along X direction (for decomposing parallel work)",
00780        &MSMBlockSizeX, 8);
00781    opts.optional("MSM", "MSMBlockSizeY",
00782        "MSM grid block size along Y direction (for decomposing parallel work)",
00783        &MSMBlockSizeY, 8);
00784    opts.optional("MSM", "MSMBlockSizeZ",
00785        "MSM grid block size along Z direction (for decomposing parallel work)",
00786        &MSMBlockSizeZ, 8);
00787 
00788    opts.optionalB("MSM", "MsmSerial",
00789        "Use MSM serial version for long-range calculation?",
00790        &MsmSerialOn, FALSE);
00791 
00792 
00794 
00795    opts.optionalB("main", "FMM",
00796        "Use fast multipole method for electrostatics?",
00797        &FMMOn, FALSE);
00798    opts.optional("FMM", "FMMLevels", "FMM number of levels",
00799        &FMMLevels, 0);
00800    opts.optional("FMM", "FMMPadding", "FMM padding margin (Angstroms)",
00801        &FMMPadding, 0);
00802 
00803 
00805 
00806    opts.optionalB("main", "PME", "Use particle mesh Ewald for electrostatics?",
00807         &PMEOn, FALSE);
00808    opts.optional("PME", "PMETolerance", "PME direct space tolerance",
00809         &PMETolerance, 1.e-6);
00810    opts.optional("PME", "PMEInterpOrder", "PME interpolation order",
00811         &PMEInterpOrder, 4);  // cubic interpolation is default
00812    opts.optional("PME", "PMEGridSizeX", "PME grid in x dimension",
00813         &PMEGridSizeX, 0);
00814    opts.optional("PME", "PMEGridSizeY", "PME grid in y dimension",
00815         &PMEGridSizeY, 0);
00816    opts.optional("PME", "PMEGridSizeZ", "PME grid in z dimension",
00817         &PMEGridSizeZ, 0);
00818    opts.optional("PME", "PMEGridSpacing", "Maximum PME grid spacing (Angstroms)",
00819         &PMEGridSpacing, 0.);
00820    opts.range("PMEGridSpacing", NOT_NEGATIVE);
00821    opts.optional("PME", "PMEProcessors",
00822         "PME FFT and reciprocal sum processor count", &PMEProcessors, 0);
00823    opts.optional("PME", "PMEMinSlices",
00824         "minimum thickness of PME reciprocal sum slab", &PMEMinSlices, 2);
00825    opts.range("PMEMinSlices", NOT_NEGATIVE);
00826    opts.optional("PME", "PMEPencils",
00827         "PME FFT and reciprocal sum pencil grid size", &PMEPencils, -1);
00828    opts.optional("PME", "PMEPencilsX",
00829         "PME FFT and reciprocal sum pencil grid size X", &PMEPencilsX, 0);
00830    opts.optional("PME", "PMEPencilsY",
00831         "PME FFT and reciprocal sum pencil grid size Y", &PMEPencilsY, 0);
00832    opts.optional("PME", "PMEPencilsZ",
00833         "PME FFT and reciprocal sum pencil grid size Z", &PMEPencilsZ, 0);
00834    opts.range("PMEPencilsX", NOT_NEGATIVE);
00835    opts.range("PMEPencilsY", NOT_NEGATIVE);
00836    opts.range("PMEPencilsZ", NOT_NEGATIVE);
00837    opts.optional("PME", "PMEPencilsYLayout",
00838         "PME FFT and reciprocal sum Y pencil layout strategy", &PMEPencilsYLayout, 0);
00839    opts.optional("PME", "PMEPencilsXLayout",
00840         "PME FFT and reciprocal sum X pencil layout strategy", &PMEPencilsXLayout, 0);
00841    opts.range("PMEPencilsYLayout", NOT_NEGATIVE);
00842    opts.range("PMEPencilsXLayout", NOT_NEGATIVE);
00843    opts.optional("PME", "PMESendOrder",
00844         "PME message ordering control", &PMESendOrder, 0);
00845    opts.range("PMESendOrder", NOT_NEGATIVE);
00846    opts.optional("PME", "PMEMinPoints",
00847         "minimum points per PME reciprocal sum pencil", &PMEMinPoints, 10000);
00848    opts.range("PMEMinPoints", NOT_NEGATIVE);
00849    opts.optionalB("main", "PMEBarrier", "Use barrier in PME?",
00850         &PMEBarrier, FALSE);
00851 
00852    opts.optionalB("PME", "useOptPME", "Use the new scalable PME optimization", &useOptPME, FALSE);
00853    opts.optionalB("PME", "useManyToMany", "Use the many-to-many PME optimization", &useManyToMany, FALSE);
00854    if (PMEOn && !useOptPME)
00855      useManyToMany = false;
00856 
00857 #ifdef DPME
00858    opts.optionalB("PME", "useDPME", "Use old DPME code?", &useDPME, FALSE);
00859 #else
00860    useDPME = 0;
00861 #endif
00862    opts.optionalB("main", "FFTWPatient", "Use intensive plan creation to optimize FFTW?",
00863 #ifdef WIN32
00864         &FFTWPatient, TRUE);
00865 #else
00866         &FFTWPatient, FALSE);
00867 #endif
00868 
00869    opts.optionalB("main", "FFTWEstimate", "Use estimates to optimize FFTW?",
00870 #ifdef WIN32
00871         &FFTWEstimate, TRUE);
00872 #else
00873         &FFTWEstimate, FALSE);
00874 #endif
00875    opts.optionalB("main", "FFTWUseWisdom", "Read/save wisdom file for FFTW?",
00876 #ifdef WIN32
00877         &FFTWUseWisdom, FALSE);
00878 #else
00879         &FFTWUseWisdom, TRUE);
00880 #endif
00881    opts.optional("FFTWUseWisdom", "FFTWWisdomFile", "File for FFTW wisdom",
00882         FFTWWisdomFile);
00883 
00884 }
00885 
00886 void SimParameters::config_parser_methods(ParseOptions &opts) {
00887    
00889    opts.optionalB("main", "minimization", "Should minimization be performed?",
00890       &minimizeCGOn, FALSE);
00891    opts.optionalB("main", "minVerbose", "Print extra minimization diagnostics?",
00892       &minVerbose, FALSE);
00893    opts.optional("main", "minTinyStep", "very first minimization steps",
00894       &minTinyStep, 1.0e-6);
00895    opts.range("minTinyStep", POSITIVE);
00896    opts.optional("main", "minBabyStep", "initial minimization steps",
00897       &minBabyStep, 1.0e-2);
00898    opts.range("minBabyStep", POSITIVE);
00899    opts.optional("main", "minLineGoal", "line minimization gradient reduction",
00900       &minLineGoal, 1.0e-3);
00901    opts.range("minLineGoal", POSITIVE);
00902 
00903    opts.optionalB("main", "velocityQuenching",
00904       "Should old-style minimization be performed?", &minimizeOn, FALSE);
00905 
00906    opts.optional("main", "maximumMove", "Maximum atom movement per step", &maximumMove, 0.0);
00907    opts.range("maximumMove", NOT_NEGATIVE);
00908    opts.units("maximumMove", N_ANGSTROM);
00909 
00910    opts.optionalB("main", "Langevin", "Should Langevin dynamics be performed?",
00911       &langevinOn, FALSE);
00912    opts.require("Langevin", "langevinTemp", "Temperature for heat bath in Langevin "
00913      "dynamics", &langevinTemp);
00914    opts.range("langevinTemp", NOT_NEGATIVE);
00915    opts.units("langevinTemp", N_KELVIN);
00916    opts.optional("Langevin", "langevinDamping", "Damping coefficient (1/ps)",
00917       &langevinDamping);
00918    opts.range("langevinDamping", POSITIVE);
00919    opts.optionalB("Langevin", "langevinHydrogen", "Should Langevin dynamics be applied to hydrogen atoms?",
00920       &langevinHydrogen);
00921    opts.optional("Langevin", "langevinFile", "PDB file with temperature "
00922      "coupling terms (B(i)) (default is the PDB input file)",
00923      PARSE_STRING);
00924    opts.optional("Langevin", "langevinCol", "Column in the langevinFile "
00925      "containing the temperature coupling term B(i);\n"
00926      "default is 'O'", PARSE_STRING);
00927 
00928    // use BAOAB integration instead of BBK
00929    opts.optionalB("Langevin", "langevinBAOAB",
00930        "Should Langevin dynamics be performed using BAOAB integration?",
00931        &langevin_useBAOAB, FALSE);
00932 
00933 // BEGIN LA
00934    opts.optionalB("main", "LoweAndersen", "Should Lowe-Andersen dynamics be performed?",
00935                   &loweAndersenOn, FALSE);
00936    opts.require("LoweAndersen", "loweAndersenTemp", "Temperature for heat bath in Lowe-Andersen "
00937                 "dynamics", &loweAndersenTemp);
00938    opts.range("loweAndersenTemp", NOT_NEGATIVE);
00939    opts.units("loweAndersenTemp", N_KELVIN);
00940    opts.optional("LoweAndersen", "loweAndersenRate", "Collision rate (1/ps)",
00941                  &loweAndersenRate, 50);
00942    opts.range("loweAndersenRate", POSITIVE);
00943    opts.optional("LoweAndersen", "loweAndersenCutoff", "Cutoff radius",
00944                  &loweAndersenCutoff, 2.7);
00945    opts.range("loweAndersenCutoff", POSITIVE);
00946    opts.units("loweAndersenCutoff", N_ANGSTROM);
00947 // END LA
00948 
00949 //Modifications for alchemical fep
00950 //  alchemical fep options
00951    opts.optionalB("main", "alch", "Is achemical simulation being performed?",
00952      &alchOn, FALSE);
00953    opts.optional("alch", "alchType", "Which alchemical method to use?", 
00954        PARSE_STRING);
00955    opts.optionalB("alch", "alchFepWCARepuOn", "WCA decomposition repu interaction in use?",
00956      &alchFepWCARepuOn, FALSE);
00957    opts.optionalB("alch", "alchFepWCADispOn", "WCA decomposition disp interaction in use?",
00958      &alchFepWCADispOn, FALSE);
00959    opts.optionalB("alch", "alchEnsembleAvg", "Ensemble Average in use?",
00960      &alchEnsembleAvg, TRUE);
00961    opts.optionalB("alch", "alchFepWhamOn", "Energy output for Wham postprocessing in use?",
00962      &alchFepWhamOn, FALSE);   
00963    opts.optional("alch", "alchFepWCArcut1", "WCA repulsion Coeff1 used for generating"
00964      "the altered alchemical vDW interactions", &alchFepWCArcut1, 0.0);
00965    opts.optional("alch", "alchFepWCArcut2", "WCA repulsion Coeff2 used for generating"
00966      "the altered alchemical vDW interactions", &alchFepWCArcut2, 1.0);
00967    opts.optional("alch", "alchFepWCArcut3", "WCA repulsion Coeff3 used for generating"
00968      "the altered alchemical vDW interactions", &alchFepWCArcut3, 1.0);
00969    opts.range("alchFepWCArcut1", NOT_NEGATIVE); 
00970    opts.range("alchFepWCArcut2", NOT_NEGATIVE);
00971    opts.range("alchFepWCArcut3", NOT_NEGATIVE);
00972 
00973    opts.optional("alch", "alchRepLambda", "Lambda of WCA repulsion"
00974      "Coupling parameter value for WCA repulsion", &alchRepLambda, -1.0);       // an invalid lambda value
00975    opts.optional("alch", "alchDispLambda", "Lambda of WCA dispersion"
00976      "Coupling parameter value for WCA dispersion", &alchDispLambda, -1.0);     // an invalid lambda value
00977    opts.optional("alch", "alchElecLambda", "Lambda of electrostatic perturbation"
00978      "Coupling parameter value for electrostatic perturbation", &alchElecLambda, -1.0); // an invalid lambda value
00979 
00980 //   opts.require("alch", "alchLambda", "Coupling parameter value", 
00981 //       &alchLambda);
00982 //   opts.require("alch", "alchLambda2", "Coupling comparison value",
00983 //       &alchLambda2);
00984    opts.optional("alch", "alchFile", "PDB file with perturbation flags "
00985      "default is the input PDB file", PARSE_STRING); 
00986    opts.optional("alch", "alchCol", "Column in the alchFile with the "
00987      "perturbation flag", PARSE_STRING);
00988    opts.optional("alch", "alchEquilSteps", "Equilibration steps, before "
00989      "data collection in the alchemical window", &alchEquilSteps, 0);
00990    opts.range("alchEquilSteps", NOT_NEGATIVE);
00991    opts.optional("alch", "alchVdwShiftCoeff", "Coeff used for generating"
00992      "the altered alchemical vDW interactions", &alchVdwShiftCoeff, 5.);
00993    opts.range("alchVdwShiftCoeff", NOT_NEGATIVE);
00994    
00995    opts.optional("alch", "alchElecLambdaStart", "Lambda at which electrostatic"
00996       "scaling of exnihilated particles begins", &alchElecLambdaStart, 0.5); 
00997    opts.range("alchElecLambdaStart", NOT_NEGATIVE);
00998    
00999    opts.optional("alch", "alchVdwLambdaEnd", "Lambda at which vdW"
01000       "scaling of exnihilated particles begins", &alchVdwLambdaEnd, 1.0); 
01001    opts.range("alchVdwLambdaEnd", NOT_NEGATIVE);  
01002 // end FEP options
01003 //fepe
01004 
01005 // REDUNDANT TI BEGINS
01006 // Modifications for TI
01007 // lots of duplication of FEP, we'd be better off without all this but 
01008 // keeping it in place for now for compatibility
01009 //   opts.optionalB("main", "thermInt", "Perform thermodynamic integration?",
01010 //     &thermInt, FALSE);
01011 // opts.require("thermInt", "tilambda", "Coupling parameter value", &tiLambda);
01012 // opts.optional("thermInt", "tiFile", "PDB file with perturbation flags "
01013 //     "default is the input PDB file", PARSE_STRING); 
01014 //   opts.optional("thermInt", "tiCol", "Column in the tiFile with the "
01015 //     "perturbation flag", PARSE_STRING);
01016 //   opts.optional("thermInt", "tiEquilSteps", "Equilibration steps, before "
01017 //     "data collection at each tiLambda value", &tiEquilSteps, 0);
01018 //   opts.range("tiEquilSteps", NOT_NEGATIVE);
01019 //   opts.optional("thermInt", "tiVdwShiftCoeff", "Coeff used for generating"
01020 //     "the altered alchemical vDW interactions", &tiVdwShiftCoeff, 5.);
01021 //   opts.range("tiVdwShiftCoeff", NOT_NEGATIVE);
01022 //   
01023 //   opts.optional("thermInt", "tiElecLambdaStart", "Lambda at which to start"
01024 //      "electrostatics scaling", &tiElecLambdaStart, 0.5); 
01025 //   opts.range("tiElecLambdaStart", NOT_NEGATIVE);
01026 //   
01027 //   opts.optional("thermInt", "tiVdwLambdaEnd", "Lambda at which to end"
01028 //      "Vdw scaling", &tiVdwLambdaEnd, 0.5); 
01029 //   opts.range("tiVdwLambdaEnd", NOT_NEGATIVE);  
01030 // end TI options
01031 // REDUNDANT TI ENDS
01032 
01033    opts.optionalB("main", "alchDecouple", "Enable alchemical decoupling?",
01034      &alchDecouple, FALSE);
01035 
01036    opts.optionalB("main", "les", "Is locally enhanced sampling enabled?",
01037      &lesOn, FALSE);
01038    opts.require("les", "lesFactor", "Local enhancement factor", &lesFactor);
01039    opts.optional("les", "lesFile", "PDB file with enhancement flags "
01040      "default is the input PDB file", PARSE_STRING); 
01041    opts.optional("les", "lesCol", "Column in the lesFile with the "
01042      "enhancement flag", PARSE_STRING);
01043    opts.optionalB("les", "lesReduceTemp", "Reduce enhanced atom temperature?",
01044      &lesReduceTemp, FALSE);
01045    opts.optionalB("les", "lesReduceMass", "Reduce enhanced atom mass?",
01046      &lesReduceMass, FALSE);
01047 
01048    // Drude oscillators
01049    opts.optionalB("main", "drude", "Perform integration of Drude oscillators?",
01050        &drudeOn, FALSE);
01051    opts.require("drude", "drudeTemp", "Temperature for freezing "
01052        "Drude oscillators", &drudeTemp);
01053    opts.range("drudeTemp", NOT_NEGATIVE);
01054    opts.units("drudeTemp", N_KELVIN);
01055    opts.optional("drude", "drudeDamping", "Damping coefficient (1/ps) for "
01056        "Drude oscillators", &drudeDamping);
01057    opts.range("drudeDamping", POSITIVE);
01058    opts.optional("drude", "drudeBondLen", "Drude oscillator bond length "
01059        "beyond which to apply restraint", &drudeBondLen);
01060    opts.range("drudeBondLen", POSITIVE);
01061    opts.optional("drude", "drudeBondConst", "Drude oscillator restraining "
01062        "force constant", &drudeBondConst);
01063    opts.range("drudeBondConst", POSITIVE);
01064    opts.optional("drude", "drudeNbtholeCut", "Nonbonded Thole interactions "
01065        "interaction radius", &drudeNbtholeCut);
01066    opts.range("drudeNbtholeCut", POSITIVE);
01067    opts.optionalB("drude", "drudeHardWall", "Apply maximum Drude bond length "
01068        "restriction?", &drudeHardWallOn, FALSE);
01069 
01070    // Pair interaction calculations
01071     opts.optionalB("main", "pairInteraction", 
01072         "Are pair interactions calculated?", &pairInteractionOn, FALSE);
01073     opts.optional("pairInteraction", "pairInteractionFile", 
01074         "PDB files with interaction flags " "default is the input PDB file", 
01075         PARSE_STRING);
01076     opts.optional("pairInteraction", "pairInteractionCol", 
01077         "Column in the pairInteractionFile with the interaction flags",
01078         PARSE_STRING);
01079     opts.require("pairInteraction", "pairInteractionGroup1",
01080         "Flag for interaction group 1", &pairInteractionGroup1);
01081     opts.optional("pairInteraction", "pairInteractionGroup2",
01082         "Flag for interaction group 2", &pairInteractionGroup2, -1);
01083     opts.optionalB("pairInteraction", "pairInteractionSelf",
01084         "Compute only within-group interactions?", &pairInteractionSelf, 
01085         FALSE);
01086    // Options for CG simulations
01087    opts.optionalB("main", "cosAngles", "Are some angles cosine-based?", &cosAngles, FALSE);
01088 
01089 
01090    //  Dihedral angle dynamics
01091    opts.optionalB("main", "globalTest", "Should global integration (for development) be used?",
01092     &globalOn, FALSE);
01093    opts.optionalB("main", "dihedral", "Should dihedral angle dynamics be performed?",
01094     &dihedralOn, FALSE);
01095    COLDOn = FALSE;
01096    opts.optionalB("dihedral", "COLD", "Should overdamped Langevin dynamics be performed?",
01097     &COLDOn, FALSE);
01098    opts.require("COLD", "COLDTemp", "Temperature for heat bath in COLD",
01099     &COLDTemp);
01100    opts.range("COLDTemp", NOT_NEGATIVE);
01101    opts.units("COLDTemp", N_KELVIN);
01102    opts.require("COLD", "COLDRate", "Damping rate for COLD",
01103     &COLDRate, 3000.0);
01104    opts.range("COLDRate", NOT_NEGATIVE);
01105 
01106    //  Get the parameters for temperature coupling
01107    opts.optionalB("main", "tcouple", 
01108       "Should temperature coupling be performed?",
01109       &tCoupleOn, FALSE);
01110    opts.require("tcouple", "tCoupleTemp", 
01111     "Temperature for temperature coupling", &tCoupleTemp);
01112    opts.range("tCoupleTemp", NOT_NEGATIVE);
01113    opts.units("tCoupleTemp", N_KELVIN);
01114    opts.optional("tCouple", "tCoupleFile", "PDB file with temperature "
01115      "coupling terms (B(i)) (default is the PDB input file)",
01116      PARSE_STRING);
01117    opts.optional("tCouple", "tCoupleCol", "Column in the tCoupleFile "
01118      "containing the temperature coupling term B(i);\n"
01119      "default is 'O'", PARSE_STRING);
01120 
01121    opts.optional("main", "rescaleFreq", "Number of steps between "
01122     "velocity rescaling", &rescaleFreq);
01123    opts.range("rescaleFreq", POSITIVE);
01124    opts.optional("main", "rescaleTemp", "Target temperature for velocity rescaling",
01125     &rescaleTemp);
01126    opts.range("rescaleTemp", NOT_NEGATIVE);
01127    opts.units("rescaleTemp", N_KELVIN);
01128 
01129    opts.optional("main", "reassignFreq", "Number of steps between "
01130     "velocity reassignment", &reassignFreq);
01131    opts.range("reassignFreq", POSITIVE);
01132    opts.optional("main", "reassignTemp", "Target temperature for velocity reassignment",
01133     &reassignTemp);
01134    opts.range("reassignTemp", NOT_NEGATIVE);
01135    opts.units("reassignTemp", N_KELVIN);
01136    opts.optional("main", "reassignIncr", "Temperature increment for velocity reassignment",
01137     &reassignIncr);
01138    opts.units("reassignIncr", N_KELVIN);
01139    opts.optional("main", "reassignHold", "Final holding temperature for velocity reassignment",
01140     &reassignHold);
01141    opts.range("reassignHold", NOT_NEGATIVE);
01142    opts.units("reassignHold", N_KELVIN);
01143 
01145    opts.optionalB("main", "useGroupPressure", 
01146       "Use group rather than atomic quantities for pressure control?",
01147       &useGroupPressure, FALSE);
01148 
01150    opts.optionalB("main", "useFlexibleCell",
01151       "Use anisotropic cell fluctuation for pressure control?",
01152       &useFlexibleCell, FALSE);
01153 
01154    // Fix specific cell dimensions
01155    opts.optionalB("main", "fixCellDims",
01156       "Fix some cell dimensions?",
01157       &fixCellDims, FALSE);
01158 
01159    opts.optionalB("fixCellDims", "fixCellDimX",
01160       "Fix the X dimension?",
01161       &fixCellDimX, FALSE);
01162    opts.optionalB("fixCellDims", "fixCellDimY",
01163       "Fix the Y dimension?",
01164       &fixCellDimY, FALSE);
01165    opts.optionalB("fixCellDims", "fixCellDimZ",
01166       "Fix the Z dimension?",
01167       &fixCellDimZ, FALSE);
01168 
01170    opts.optionalB("main", "useConstantRatio",
01171       "Use constant X-Y ratio for pressure control?",
01172       &useConstantRatio, FALSE);
01173 
01175    opts.optionalB("main", "useConstantArea",
01176       "Use constant area for pressure control?",
01177       &useConstantArea, FALSE);
01178 
01180    opts.optionalB("main", "excludeFromPressure",
01181         "Should some atoms be excluded from pressure rescaling?",
01182         &excludeFromPressure, FALSE);
01183    opts.optional("excludeFromPressure", "excludeFromPressureFile",
01184         "PDB file for atoms to be excluded from pressure",
01185         PARSE_STRING);
01186    opts.optional("excludeFromPressure", "excludeFromPressureCol", 
01187         "Column in the excludeFromPressureFile"
01188         "containing the flags (nonzero means excluded);\n"
01189         "default is 'O'", PARSE_STRING);
01190 
01192    opts.optionalB("main", "BerendsenPressure", 
01193       "Should Berendsen pressure bath coupling be performed?",
01194       &berendsenPressureOn, FALSE);
01195    opts.require("BerendsenPressure", "BerendsenPressureTarget",
01196     "Target pressure for pressure coupling",
01197     &berendsenPressureTarget);
01198    // opts.units("BerendsenPressureTarget",);
01199    opts.require("BerendsenPressure", "BerendsenPressureCompressibility",
01200     "Isothermal compressibility for pressure coupling",
01201     &berendsenPressureCompressibility);
01202    // opts.units("BerendsenPressureCompressibility",);
01203    opts.require("BerendsenPressure", "BerendsenPressureRelaxationTime",
01204     "Relaxation time for pressure coupling",
01205     &berendsenPressureRelaxationTime);
01206    opts.range("BerendsenPressureRelaxationTime", POSITIVE);
01207    opts.units("BerendsenPressureRelaxationTime", N_FSEC);
01208    opts.optional("BerendsenPressure", "BerendsenPressureFreq",
01209     "Number of steps between volume rescaling",
01210     &berendsenPressureFreq, 1);
01211    opts.range("BerendsenPressureFreq", POSITIVE);
01212 
01214    opts.optionalB("main", "LangevinPiston",
01215       "Should Langevin piston pressure control be used?",
01216       &langevinPistonOn, FALSE);
01217    opts.optionalB("LangevinPiston", "LangevinPistonBarrier",
01218       "Should Langevin piston barrier be used?",
01219       &langevinPistonBarrier, TRUE);
01220    opts.require("LangevinPiston", "LangevinPistonTarget",
01221       "Target pressure for pressure control",
01222       &langevinPistonTarget);
01223    opts.require("LangevinPiston", "LangevinPistonPeriod",
01224       "Oscillation period for pressure control",
01225       &langevinPistonPeriod);
01226    opts.range("LangevinPistonPeriod", POSITIVE);
01227    opts.units("LangevinPistonPeriod", N_FSEC);
01228    opts.require("LangevinPiston", "LangevinPistonDecay",
01229       "Decay time for pressure control",
01230       &langevinPistonDecay);
01231    opts.range("LangevinPistonDecay", POSITIVE);
01232    opts.units("LangevinPistonDecay", N_FSEC);
01233    opts.require("LangevinPiston", "LangevinPistonTemp",
01234       "Temperature for pressure control piston",
01235       &langevinPistonTemp);
01236    opts.range("LangevinPistonTemp", POSITIVE);
01237    opts.units("LangevinPistonTemp", N_KELVIN);
01238    opts.optional("LangevinPiston", "StrainRate",
01239       "Initial strain rate for pressure control (x y z)",
01240       &strainRate);
01241 
01243    opts.optional("main", "SurfaceTensionTarget",
01244       "Surface tension in the x-y plane",
01245       &surfaceTensionTarget, 0);
01246 
01248    opts.optionalB("main", "pressureprofile", "Compute pressure profile?",
01249      &pressureProfileOn, FALSE);
01250    opts.require("pressureprofile", "pressureprofileslabs", 
01251      "Number of pressure profile slabs", &pressureProfileSlabs, 10);
01252    opts.optional("pressureprofile", "pressureprofilefreq",
01253      "How often to store profile data", &pressureProfileFreq, 1);
01254    opts.optional("pressureprofile", "pressureProfileAtomTypes", 
01255      "Number of pressure profile atom types", &pressureProfileAtomTypes, 1);
01256    opts.range("pressureProfileAtomTypes", POSITIVE);
01257    opts.optional("pressureProfile", "pressureProfileAtomTypesFile", 
01258         "PDB files with pressure profile atom types" "default is the input PDB file", 
01259         PARSE_STRING);
01260    opts.optional("pressureProfile", "pressureProfileAtomTypesCol", 
01261         "Column in the pressureProfileAtomTypesFile with the atom types ",
01262         PARSE_STRING);
01263    opts.optionalB("pressureProfile", "pressureProfileEwald", 
01264        "Compute Ewald contribution to pressure profile",
01265        &pressureProfileEwaldOn, FALSE);
01266    opts.optional("pressureProfile", "pressureProfileEwaldX",
01267        "Ewald grid size X", &pressureProfileEwaldX, 10);
01268    opts.range("pressureProfileEwaldX", POSITIVE);
01269    opts.optional("pressureProfile", "pressureProfileEwaldY",
01270        "Ewald grid size Y", &pressureProfileEwaldY, 10);
01271    opts.range("pressureProfileEwaldY", POSITIVE);
01272    opts.optional("pressureProfile", "pressureProfileEwaldZ",
01273        "Ewald grid size Z", &pressureProfileEwaldZ, 10);
01274    opts.range("pressureProfileEwaldZ", POSITIVE);
01275 
01277    opts.optionalB("main", "accelMD", "Perform acclerated MD?", &accelMDOn, FALSE);
01278    opts.optional("accelMD", "accelMDFirstStep", "First accelMD step", &accelMDFirstStep, 0);
01279    opts.range("accelMDFirstStep", NOT_NEGATIVE);
01280    opts.optional("accelMD", "accelMDLastStep", "Last accelMD step", &accelMDLastStep, 0);
01281    opts.range("accelMDLastStep", NOT_NEGATIVE);
01282    opts.optional("accelMD", "accelMDOutFreq", "Frequency of accelMD output", &accelMDOutFreq, 1);
01283    opts.range("accelMDOutFreq", POSITIVE);
01284    opts.optionalB("accelMD", "accelMDdihe", "Apply boost to dihedral potential", &accelMDdihe, TRUE);
01285    opts.optionalB("accelMD", "accelMDDebugOn", "Debugging accelMD", &accelMDDebugOn, FALSE);
01286    opts.require("accelMD", "accelMDE","E for AMD", &accelMDE);
01287    opts.units("accelMDE", N_KCAL);
01288    opts.require("accelMD", "accelMDalpha","alpha for AMD", &accelMDalpha);
01289    opts.units("accelMDalpha", N_KCAL);
01290    opts.range("accelMDalpha", POSITIVE);
01291    opts.optionalB("accelMD", "accelMDdual", "Apply dual boost", &accelMDdual, FALSE);
01292    opts.require("accelMDdual", "accelMDTE","E for total potential under accelMDdual mode", &accelMDTE);
01293    opts.units("accelMDTE", N_KCAL);
01294    opts.require("accelMDdual", "accelMDTalpha","alpha for total potential under accelMDdual mode", &accelMDTalpha);
01295    opts.units("accelMDTalpha", N_KCAL);
01296    opts.range("accelMDTalpha", POSITIVE);
01297 
01298    // Adaptive Temperature Sampling (adaptTemp) parameters
01299    opts.optionalB("main", "adaptTempMD", "Perform adaptive temperature sampling", &adaptTempOn, FALSE);
01300    opts.optional("adaptTempMD", "adaptTempFirstStep", "First adaptTemp step", &adaptTempFirstStep, 0);
01301    opts.range("adaptTempFirstStep", NOT_NEGATIVE);
01302    opts.optional("adaptTempMD", "adaptTempLastStep", "Last adaptTemp step", &adaptTempLastStep, 0);
01303    opts.range("adaptTempLastStep", NOT_NEGATIVE);
01304    opts.optional("adaptTempMD", "adaptTempOutFreq", "Frequency of adaptTemp output", &adaptTempOutFreq, 10);
01305    opts.range("adaptTempOutFreq", POSITIVE);
01306    opts.optional("adaptTempMD", "adaptTempFreq", "Frequency of writing average energies to adaptTempOutFile", &adaptTempFreq, 10);
01307    opts.range("adaptTempFreq", POSITIVE);
01308    opts.optionalB("adaptTempMD", "adaptTempDebug", "Print debug output for adaptTemp", &adaptTempDebug, FALSE);
01309    opts.optional("adaptTempMD", "adaptTempTmin","Minimun temperature for adaptTemp", &adaptTempTmin);
01310    opts.units("adaptTempTmin", N_KELVIN);
01311    opts.range("adaptTempTmin", POSITIVE);
01312    opts.optional("adaptTempMD", "adaptTempTmax","Maximum temperature for adaptTemp", &adaptTempTmax);
01313    opts.units("adaptTempTmax", N_KELVIN);
01314    opts.range("adaptTempTmax", POSITIVE);
01315    opts.optional("adaptTempMD", "adaptTempBins","Number of bins to store average energies", &adaptTempBins,0);
01316    opts.range("adaptTempBins", NOT_NEGATIVE);
01317    opts.optional("adaptTempMD", "adaptTempDt", "Integration timestep for Temp. updates", &adaptTempDt, 0.0001);
01318    opts.units("adaptTempDt", N_FSEC);
01319    opts.range("adaptTempDt", NOT_NEGATIVE);
01320    opts.optional("adaptTempMD", "adaptTempAutoDt", "Average temperature update in percent of temperature range", &adaptTempAutoDt, 0.0);
01321    opts.range("adaptTempAutoDt", NOT_NEGATIVE);
01322    opts.optional("adaptTempMD", "adaptTempCgamma", "Adaptive bin averaging constant", &adaptTempCgamma, 0.1);
01323    opts.range("adaptTempCgamma", NOT_NEGATIVE);
01324    opts.optionalB("adaptTempMD","adaptTempLangevin","Send adaptTemp temperature to langevin thermostat",&adaptTempLangevin,TRUE);
01325    opts.optionalB("adaptTempMD","adaptTempRescaling","Send adaptTemp temperature to velocity rescaling thermostat", &adaptTempRescale,TRUE);
01326    opts.optional("adaptTempMD", "adaptTempInFile", "File containing restart information for adaptTemp", adaptTempInFile);
01327    opts.optional("adaptTempMD", "adaptTempRestartFile", "File for writing adaptTemp restart information", adaptTempRestartFile);
01328    opts.require("adaptTempRestartFile","adaptTempRestartFreq", "Frequency of writing restart file", &adaptTempRestartFreq,0);
01329    opts.range("adaptTempRestartFreq",NOT_NEGATIVE);
01330    opts.optionalB("adaptTempMD", "adaptTempRandom", "Randomly assign a temperature if we step out of range", &adaptTempRandom, FALSE);
01331 }
01332 
01333 void SimParameters::config_parser_constraints(ParseOptions &opts) {
01334    
01336    opts.optionalB("main", "fixedatoms", "Are there fixed atoms?",
01337     &fixedAtomsOn, FALSE);
01338    opts.optionalB("fixedatoms", "fixedAtomsForces",
01339      "Calculate forces between fixed atoms?  (Required to unfix during run.)",
01340      &fixedAtomsForces, FALSE);
01341    opts.optional("fixedatoms", "fixedAtomsFile", "PDB file with flags for "
01342      "fixed atoms (default is the PDB input file)",
01343      PARSE_STRING);
01344    opts.optional("fixedatoms", "fixedAtomsCol", "Column in the fixedAtomsFile "
01345      "containing the flags (nonzero means fixed);\n"
01346      "default is 'O'", PARSE_STRING);
01347    opts.optional("fixedatoms", "fixedAtomListFile", "the text input file for fixed atoms "
01348                  "used for parallel input IO", PARSE_STRING);
01349 
01351    opts.optionalB("main", "constraints", "Are harmonic constraints active?",
01352      &constraintsOn, FALSE);
01353    opts.require("constraints", "consexp", "Exponent for harmonic potential",
01354     &constraintExp, 2);
01355    opts.range("consexp", POSITIVE);
01356 #ifndef MEM_OPT_VERSION
01357    opts.require("constraints", "consref", "PDB file containing reference "
01358     "positions",
01359     PARSE_STRING);
01360    opts.require("constraints", "conskfile", "PDB file containing force "
01361     "constaints in one of the columns", PARSE_STRING);
01362    opts.require("constraints", "conskcol", "Column of conskfile to use "
01363     "for the force constants", PARSE_STRING);
01364 #else
01365    opts.require("constraints", "consAtomListFile", "the text input file for constrained atoms "
01366                  "used for parallel input IO", PARSE_STRING);
01367 #endif
01368    opts.require("constraints", "constraintScaling", "constraint scaling factor",
01369      &constraintScaling, 1.0);
01370    opts.range("constraintScaling", NOT_NEGATIVE);
01371 
01372 
01373 
01374    //****** BEGIN selective restraints (X,Y,Z) changes
01375 
01377    opts.optionalB("constraints", "selectConstraints", 
01378    "Restrain only selected Cartesian components of the coordinates?",
01379      &selectConstraintsOn, FALSE);
01380    opts.optionalB("selectConstraints", "selectConstrX",  
01381    "Restrain X components of coordinates ", &constrXOn, FALSE);
01382    opts.optionalB("selectConstraints", "selectConstrY",  
01383    "Restrain Y components of coordinates ", &constrYOn, FALSE);
01384    opts.optionalB("selectConstraints", "selectConstrZ",  
01385    "Restrain Z components of coordinates ", &constrZOn, FALSE);
01386    //****** END selective restraints (X,Y,Z) changes
01387 
01388    // spherical constraints
01389    opts.optionalB("constraints", "sphericalConstraints", 
01390    "Restrain only radial spherical component of the coordinates?",
01391      &sphericalConstraintsOn, FALSE);
01392    opts.optional("sphericalConstraints", "sphericalConstrCenter",
01393    "Center of spherical constraints", &sphericalConstrCenter);
01394  
01395    //****** BEGIN moving constraints changes 
01396 
01398    opts.optionalB("constraints", "movingConstraints",
01399       "Are some of the constraints moving?", 
01400       &movingConstraintsOn, FALSE);
01401    opts.require("movingConstraints", "movingConsVel",
01402     "Velocity of the movement, A/timestep", &movingConsVel);
01403    //****** END moving constraints changes 
01404 
01405    // BEGIN rotating constraints changes
01406    opts.optionalB("constraints", "rotConstraints",
01407       "Are the constraints rotating?", 
01408       &rotConstraintsOn, FALSE);
01409    opts.require("rotConstraints", "rotConsAxis",
01410     "Axis of rotation", &rotConsAxis);
01411    opts.require("rotConstraints", "rotConsPivot",
01412     "Pivot point of rotation", 
01413     &rotConsPivot);
01414    opts.require("rotConstraints", "rotConsVel",
01415     "Velocity of rotation, deg/timestep", &rotConsVel);
01416 
01417    // END rotating constraints changes
01418 
01419    // external command forces
01420    opts.optionalB("main", "extForces", "External command forces?",
01421       &extForcesOn, FALSE);
01422    opts.require("extForces", "extForcesCommand",
01423       "External forces command", extForcesCommand);
01424    opts.require("extForces", "extCoordFilename",
01425       "External forces coordinate filename", extCoordFilename);
01426    opts.require("extForces", "extForceFilename",
01427       "External forces force filename", extForceFilename);
01428 
01429    //print which bad contacts are being moved downhill
01430    opts.optionalB("main", "printBadContacts", "Print atoms with huge forces?",
01431       &printBadContacts, FALSE);
01432 
01433    /* GBIS generalized born implicit solvent*/
01434 
01435    opts.optionalB("main", "GBIS", "Use GB implicit solvent?",
01436       &GBISOn, FALSE);
01437    opts.optionalB("main", "GBISSer", "Use GB implicit solvent?",
01438       &GBISserOn, FALSE);
01439 
01440    opts.optional("GBIS", "solventDielectric",
01441       "Solvent Dielectric", &solvent_dielectric, 78.5);
01442    opts.optional("GBIS", "intrinsicRadiusOffset",
01443       "Coulomb Radius Offset", &coulomb_radius_offset, 0.09);
01444    opts.optional("GBIS", "ionConcentration",
01445       "Ion Concentration", &ion_concentration, 0.2); //0.2 mol/L
01446    opts.optional("GBIS", "GBISDelta",
01447       "delta from GBOBC", &gbis_delta, 1.0); //0.8 or 1.0
01448    opts.optional("GBIS", "GBISBeta",
01449       "beta from GBOBC", &gbis_beta, 0.8);   //0.0 or 0.8
01450    opts.optional("GBIS", "GBISGamma",
01451       "gamma from GBOBC", &gbis_gamma, 4.85);//2.290912 or 4.85
01452    opts.optional("GBIS", "alphaCutoff",
01453       "cutoff for calculating effective born radius", &alpha_cutoff, 15);
01454    opts.optional("GBIS", "alphaMax",
01455       "maximum allowable born radius", &alpha_max, 30);
01456    opts.optional("GBIS", "fsMax",
01457       "maximum screened intrinsic radius", &fsMax, 1.728);
01458 
01459    opts.optionalB("main", "SASA", "Use Linear Combination of Pairwise Overlaps (LCPO) for calculating SASA",
01460       &LCPOOn, FALSE);
01461    opts.optional("SASA", "surfaceTension",
01462       "Surfce Tension for SASA (kcal/mol/Ang^2)", &surface_tension, 0.005);
01463 
01464    //****** BEGIN SMD constraints changes 
01465 
01466    // SMD constraints
01467    opts.optionalB("main", "SMD",
01468       "Do we use SMD option?", 
01469       &SMDOn, FALSE);
01470    opts.require("SMD", "SMDVel",
01471                 "Velocity of the movement, A/timestep", &SMDVel);
01472    opts.range("SMDVel", NOT_NEGATIVE);
01473    opts.require("SMD", "SMDDir",
01474                 "Direction of movement", &SMDDir);
01475    opts.require("SMD", "SMDk",
01476                 "Elastic constant for SMD", &SMDk);
01477    opts.optional("SMD", "SMDk2",
01478                 "Transverse elastic constant for SMD", &SMDk2, 0);
01479    opts.range("SMDk", NOT_NEGATIVE);
01480    opts.range("SMDk2", NOT_NEGATIVE);
01481    opts.require("SMD", "SMDFile",
01482                 "File for SMD information",
01483                  SMDFile);
01484    opts.optional("SMD", "SMDOutputFreq",
01485                  "Frequency of output",
01486                  &SMDOutputFreq, 1);
01487    opts.range("SMDOutputFreq", POSITIVE);
01488    
01489    //****** END SMD constraints changes 
01490 
01491    //****** BEGIN tabulated energies section
01492    opts.optionalB("main", "tabulatedEnergies", "Do we get energies from a table?", &tabulatedEnergies, FALSE);
01493 //   opts.require("tabulatedEnergies", "tableNumTypes","Number of types for energy tabulation", &tableNumTypes);
01494    opts.require("tabulatedEnergies", "tabulatedEnergiesFile", "File containing energy table", tabulatedEnergiesFile);
01495    opts.require("tabulatedEnergies", "tableInterpType", "Cubic or linear interpolation", tableInterpType);
01496 
01497    // TMD parameters
01498    opts.optionalB("main", "TMD", "Perform Targeted MD?", &TMDOn, FALSE);
01499    opts.optional("TMD", "TMDk", "Elastic constant for TMD", &TMDk, 0); 
01500    opts.range("TMDk", NOT_NEGATIVE);
01501    opts.require("TMD", "TMDFile", "File for TMD information", TMDFile);
01502    opts.optionalB("TMD", "TMDDiffRMSD", "Restrain Difference between the RMSD from two structures", &TMDDiffRMSD, FALSE);
01503    opts.require("TMDDiffRMSD", "TMDFile2",  "Second file for TMD information", TMDFile2); 
01504     
01505    opts.optional("TMD", "TMDOutputFreq", "Frequency of TMD output", 
01506        &TMDOutputFreq, 1);
01507    opts.range("TMDOutputFreq", POSITIVE);
01508    opts.require("TMD", "TMDLastStep", "Last TMD timestep", &TMDLastStep);
01509    opts.range("TMDLastStep", POSITIVE);
01510    opts.optional("TMD", "TMDFirstStep", "First TMD step (default 0)", &TMDFirstStep, 0);
01511    opts.optional("TMD", "TMDInitialRMSD", "Target RMSD at first TMD step (default -1 to use initial coordinates)", &TMDInitialRMSD);
01512    TMDInitialRMSD = -1;
01513    opts.optional("TMD", "TMDFinalRMSD", "Target RMSD at last TMD step (default 0 )", &TMDFinalRMSD, 0);
01514    opts.range("TMDInitialRMSD", NOT_NEGATIVE);
01515    // End of TMD parameters
01516 
01517    // Symmetry restraint parameters
01518    opts.optionalB("main", "symmetryRestraints", "Enable symmetry restraints?", &symmetryOn, FALSE); 
01519    opts.optional("symmetryRestraints", "symmetryk", "Elastic constant for symmetry restraints", &symmetryk, 0);
01520    opts.range("symmetryk", NOT_NEGATIVE);
01521    opts.optional("symmetryRestraints", "symmetrykfile", "PDB file specifying force contants on a per-atom basis", PARSE_MULTIPLES);
01522    opts.optionalB("symmetryRestraints", "symmetryScaleForces", "Scale applied forces over time?", &symmetryScaleForces, FALSE);
01523    opts.require("symmetryRestraints", "symmetryFile", "File for symmetry information", PARSE_MULTIPLES);
01524    opts.optional("symmetryRestraints", "symmetryMatrixFile", "File(s) for transfromation matrices", PARSE_MULTIPLES);
01525    opts.optional("symmetryRestraints", "symmetryLastStep", "Last symmetry timestep", &symmetryLastStep, -1);
01526    opts.optional("symmetryRestraints", "symmetryFirstStep", "First symmetry step (default 0)", &symmetryFirstStep, 0);
01527    opts.optional("symmetryRestraints", "symmetryLastFullStep", "Last full force symmetry timestep (default symmetryLastStep)", &symmetryLastFullStep, symmetryLastStep);
01528    opts.optional("symmetryRestraints", "symmetryFirstFullStep", "First full force symmetry step (default symmetryFirstStep)", &symmetryFirstFullStep, symmetryFirstStep);
01529   //End of symmetry restraint parameters.
01530 
01532    opts.optionalB("main", "tclForces", "Are Tcl global forces active?",
01533      &tclForcesOn, FALSE);
01534    opts.require("tclForces", "tclForcesScript",
01535      "Tcl script for global forces", PARSE_MULTIPLES);
01536 
01538    opts.optionalB("main", "tclBC", "Are Tcl boundary forces active?",
01539      &tclBCOn, FALSE);
01540    opts.require("tclBC", "tclBCScript",
01541      "Tcl script defining calcforces for boundary forces", PARSE_STRING);
01542    tclBCScript = 0;
01543    opts.optional("tclBC", "tclBCArgs", "Extra args for calcforces command",
01544      tclBCArgs);
01545    tclBCArgs[0] = 0;
01546 
01548    opts.optionalB("main", "miscForces", "Are misc global forces active?",
01549      &miscForcesOn, FALSE);
01550    opts.optional("miscForces", "miscForcesScript",
01551      "script for misc forces", PARSE_MULTIPLES);
01552 
01554    opts.optionalB("main", "freeEnergy", "Perform free energy perturbation?",
01555      &freeEnergyOn, FALSE);
01556    opts.require("freeEnergy", "freeEnergyConfig",
01557      "Configuration file for free energy perturbation", PARSE_MULTIPLES);
01558 
01560    opts.optionalB("main", "constantforce", "Apply constant force?",
01561      &consForceOn, FALSE);
01562    opts.optional("constantforce", "consForceFile",
01563        "Configuration file for constant forces", PARSE_STRING);
01564    opts.require("constantforce", "consForceScaling",
01565        "Scaling factor for constant forces", &consForceScaling, 1.0);
01566  
01568     opts.optionalB("main", "colvars", "Is the colvars module enabled?",
01569       &colvarsOn, FALSE);
01570     opts.require("colvars", "colvarsConfig",
01571       "configuration for the collective variables", PARSE_STRING);
01572     opts.optional("colvars", "colvarsInput",
01573       "input restart file for the collective variables", PARSE_STRING);
01574 
01575 }
01576 
01577 #ifdef OPENATOM_VERSION
01578 void SimParameters::config_parser_openatom(ParseOptions &opts) {
01579   opts.optionalB("main", "openatom", "OpenAtom active?", &openatomOn, FALSE);
01580   opts.require("openatom", "openatomDriverFile", "What config file specifies openatom input parameters", PARSE_STRING);
01581   opts.require("openatom", "openatomPhysicsFile", "What structure file specifies openatom input system", PARSE_STRING);
01582   opts.require("openatom", "openatomPdbFile", "NAMD input file defining QM and MM regions", PARSE_STRING);
01583    opts.optional("openatom", "openatomCol", "Column in the openatomPdb with the QM/MM flag", PARSE_STRING);
01584 }
01585 #endif // OPENATOM_VERSION
01586 
01587 /* BEGIN gf */
01588 void SimParameters::config_parser_mgridforce(ParseOptions &opts) {
01590     opts.optionalB("main", "mgridforce", "Is Multiple gridforce active?", 
01591                    &mgridforceOn, FALSE);
01592     opts.optional("mgridforce", "mgridforcevolts", "Is Gridforce using Volts/eV as units?",
01593                   PARSE_MULTIPLES);
01594     opts.require("mgridforce", "mgridforcescale", "Scale factor by which to multiply "
01595                  "grid forces", PARSE_MULTIPLES);
01596     opts.require("mgridforce", "mgridforcefile", "PDB file containing force "
01597                  "multipliers in one of the columns", PARSE_MULTIPLES);
01598     opts.require("mgridforce", "mgridforcecol", "Column of gridforcefile to "
01599                  "use for force multiplier", PARSE_MULTIPLES);
01600     opts.optional("mgridforce", "mgridforcechargecol", "Column of gridforcefile to "
01601                   "use for charge", PARSE_MULTIPLES);
01602     opts.require("mgridforce", "mgridforcepotfile", "Gridforce potential file",
01603                  PARSE_MULTIPLES);
01604     opts.optional("mgridforce", "mgridforcecont1", "Use continuous grid "
01605                    "in K1 direction?", PARSE_MULTIPLES);
01606     opts.optional("mgridforce", "mgridforcecont2", "Use continuous grid "
01607                    "in K2 direction?", PARSE_MULTIPLES);
01608     opts.optional("mgridforce", "mgridforcecont3", "Use continuous grid "
01609                    "in K3 direction?", PARSE_MULTIPLES);
01610     opts.optional("mgridforce", "mgridforcevoff", "Gridforce potential offsets",
01611                   PARSE_MULTIPLES);
01612     opts.optional("mgridforce", "mgridforcelite", "Use Gridforce Lite?",
01613                   PARSE_MULTIPLES);
01614 }
01615 
01616 void SimParameters::config_parser_gridforce(ParseOptions &opts) {
01618     opts.optionalB("main", "gridforce", "Is Gridforce active?", 
01619                    &gridforceOn, FALSE);
01620     opts.optionalB("gridforce", "gridforcevolts", "Is Gridforce using Volts/eV as units?",
01621                    &gridforceVolts, FALSE);
01622     opts.require("gridforce", "gridforcescale", "Scale factor by which to multiply "
01623                  "grid forces", &gridforceScale);
01624     opts.require("gridforce", "gridforcefile", "PDB file containing force "
01625                  "multipliers in one of the columns", PARSE_STRING);
01626     opts.require("gridforce", "gridforcecol", "Column of gridforcefile to "
01627                  "use for force multiplier", PARSE_STRING);
01628     opts.optional("gridforce", "gridforcechargecol", "Column of gridforcefile to "
01629                   "use for charge", PARSE_STRING);
01630     opts.require("gridforce", "gridforcepotfile", "Gridforce potential file",
01631                  PARSE_STRING);
01632     opts.optionalB("gridforce", "gridforcecont1", "Use continuous grid "
01633                    "in A1 direction?", &gridforceContA1, FALSE);
01634     opts.optionalB("gridforce", "gridforcecont2", "Use continuous grid "
01635                    "in A2 direction?", &gridforceContA2, FALSE);
01636     opts.optionalB("gridforce", "gridforcecont3", "Use continuous grid "
01637                    "in A3 direction?", &gridforceContA3, FALSE);
01638     opts.optional("gridforce", "gridforcevoff", "Gridforce potential offsets",
01639                   &gridforceVOffset);
01640     opts.optionalB("gridforce", "gridforcelite", "Use Gridforce Lite?",
01641                    &gridforceLite, FALSE);
01642 }
01643 /* END gf */
01644 
01645 void SimParameters::config_parser_movdrag(ParseOptions &opts) {
01647    opts.optionalB("main", "movDragOn", "Do we apply moving drag?",
01648       &movDragOn, FALSE);
01649    opts.require("movDragOn", "movDragFile",
01650       "Main moving drag PDB file", movDragFile);
01651    opts.require("movDragOn", "movDragCol",
01652       "Main moving drag PDB column", PARSE_STRING);
01653    opts.require("movDragOn", "movDragGlobVel",
01654       "Global moving drag velocity (A/step)", &movDragGlobVel);
01655    opts.require("movDragOn", "movDragVelFile",
01656       "Moving drag linear velocity file", movDragVelFile);
01657 }
01658 
01659 void SimParameters::config_parser_rotdrag(ParseOptions &opts) {
01661    opts.optionalB("main", "rotDragOn", "Do we apply rotating drag?",
01662       &rotDragOn, FALSE);
01663    opts.require("rotDragOn", "rotDragFile",
01664       "Main rotating drag PDB file", rotDragFile);
01665    opts.require("rotDragOn", "rotDragCol",
01666       "Main rotating drag PDB column", PARSE_STRING);
01667    opts.require("rotDragOn", "rotDragAxisFile",
01668       "Rotating drag axis file", rotDragAxisFile);
01669    opts.require("rotDragOn", "rotDragPivotFile",
01670       "Rotating drag pivot point file", rotDragPivotFile);
01671    opts.require("rotDragOn", "rotDragGlobVel",
01672       "Global rotating drag angular velocity (deg/step)", &rotDragGlobVel);
01673    opts.require("rotDragOn", "rotDragVelFile",
01674       "Rotating drag angular velocity file", rotDragVelFile);
01675    opts.require("rotDragOn", "rotDragVelCol",
01676       "Rotating drag angular velocity column", PARSE_STRING);
01677 }
01678 
01679 void SimParameters::config_parser_constorque(ParseOptions &opts) {
01681    opts.optionalB("main", "consTorqueOn", "Do we apply \"constant\" torque?",
01682       &consTorqueOn, FALSE);
01683    opts.require("consTorqueOn", "consTorqueFile",
01684       "Main \"constant\" torque PDB file", consTorqueFile);
01685    opts.require("consTorqueOn", "consTorqueCol",
01686       "Main \"constant\" torque PDB column", PARSE_STRING);
01687    opts.require("consTorqueOn", "consTorqueAxisFile",
01688       "\"Constant\" torque axis file", consTorqueAxisFile);
01689    opts.require("consTorqueOn", "consTorquePivotFile",
01690       "\"Constant\" torque pivot point file", consTorquePivotFile);
01691    opts.require("consTorqueOn", "consTorqueGlobVal",
01692       "Global \"constant\" torque value (Kcal/(mol*A^2))", &consTorqueGlobVal);
01693    opts.require("consTorqueOn", "consTorqueValFile",
01694       "\"constant\" torque factors file", consTorqueValFile);
01695    opts.require("consTorqueOn", "consTorqueValCol",
01696       "\"constant\" torque factors column", PARSE_STRING);
01697 }
01698 
01699 void SimParameters::config_parser_boundary(ParseOptions &opts) {
01700     
01702    opts.optionalB("main", "sphericalBC", "Are spherical boundary counditions "
01703       "active?", &sphericalBCOn, FALSE);
01704    opts.require("sphericalBC", "sphericalBCCenter",
01705      "Center of spherical boundaries", &sphericalCenter);
01706    opts.require("sphericalBC", "sphericalBCr1", "Radius for first sphere "
01707      "potential", &sphericalBCr1);
01708    opts.range("sphericalBCr1", POSITIVE);
01709    opts.units("sphericalBCr1", N_ANGSTROM);
01710    opts.require("sphericalBC", "sphericalBCk1", "Force constant for first "
01711     "sphere potential (+ is an inward force, - outward)",
01712     &sphericalBCk1);
01713    opts.units("sphericalBCk1", N_KCAL);
01714    opts.optional("sphericalBC", "sphericalBCexp1", "Exponent for first "
01715     "sphere potential", &sphericalBCexp1, 2);
01716    opts.range("sphericalBCexp1", POSITIVE);
01717    
01718    opts.optional("sphericalBCr1", "sphericalBCr2", "Radius for second sphere "
01719      "potential", &sphericalBCr2);
01720    opts.range("sphericalBCr2", POSITIVE);
01721    opts.units("sphericalBCr2", N_ANGSTROM);
01722    opts.require("sphericalBCr2", "sphericalBCk2", "Force constant for second "
01723     "sphere potential (+ is an inward force, - outward)",
01724     &sphericalBCk2);
01725    opts.units("sphericalBCk2", N_KCAL);
01726    opts.optional("sphericalBCr2", "sphericalBCexp2", "Exponent for second "
01727     "sphere potential", &sphericalBCexp2, 2);
01728    opts.range("sphericalBCexp2", POSITIVE);
01729 
01731    opts.optionalB("main", "cylindricalBC", "Are cylindrical boundary counditions "
01732                   "active?", &cylindricalBCOn, FALSE);
01733    opts.require("cylindricalBC", "cylindricalBCr1", "Radius for first cylinder "
01734                  "potential", &cylindricalBCr1);
01735    opts.range("cylindricalBCr1", POSITIVE);
01736    opts.units("cylindricalBCr1", N_ANGSTROM);
01737    opts.require("cylindricalBC", "cylindricalBCk1", "Force constant for first "
01738                 "cylinder potential (+ is an inward force, - outward)",
01739                 &cylindricalBCk1);
01740    opts.units("cylindricalBCk1", N_KCAL);
01741    opts.optional("cylindricalBC", "cylindricalBCexp1", "Exponent for first "
01742                 "cylinder potential", &cylindricalBCexp1, 2);
01743    opts.range("cylindricalBCexp1", POSITIVE);
01744 
01745 
01746 // additions beyond those already found in spherical parameters    JJU
01747    opts.optional("cylindricalBC", "cylindricalBCAxis", "Cylinder axis (defaults to x)",
01748     PARSE_STRING);
01749    opts.require("cylindricalBC", "cylindricalBCCenter",
01750      "Center of cylindrical boundaries", &cylindricalCenter);
01751    opts.require ("cylindricalBC", "cylindricalBCl1", "Length of first cylinder",
01752                  &cylindricalBCl1);
01753    opts.range("cylindricalBCl1", POSITIVE);
01754    opts.units("cylindricalBCl1", N_ANGSTROM);
01755    opts.optional ("cylindricalBCl1", "cylindricalBCl2", "Length of second cylinder",
01756                   &cylindricalBCl2);
01757    opts.range ("cylindricalBCl2", POSITIVE);
01758    opts.units ("cylindricalBCl2", N_ANGSTROM);
01759 // end  additions
01760 
01761    opts.optional("cylindricalBCr1", "cylindricalBCr2", "Radius for second cylinder "
01762                  "potential", &cylindricalBCr2);
01763    opts.range("cylindricalBCr2", POSITIVE);
01764    opts.units("cylindricalBCr2", N_ANGSTROM);
01765    opts.require("cylindricalBCr2", "cylindricalBCk2", "Force constant for second "
01766                 "cylinder potential (+ is an inward force, - outward)",
01767                 &cylindricalBCk2);
01768    opts.units("cylindricalBCk2", N_KCAL);
01769    opts.optional("cylindricalBCr2", "cylindricalBCexp2", "Exponent for second "
01770                 "cylinder potential", &cylindricalBCexp2, 2);
01771    opts.range("cylindricalBCexp2", POSITIVE);
01772 
01774    opts.optionalB("main", "eFieldOn", "Should an electric field be applied",
01775                  &eFieldOn, FALSE);
01776    opts.require("eFieldOn", "eField", "Electric field vector", &eField);
01777    opts.optional("eFieldOn", "eFieldFreq", "Electric field frequency", &eFieldFreq);
01778    opts.optional("eFieldOn", "eFieldPhase", "Electric field phase", &eFieldPhase);
01779 
01781    opts.optionalB("main", "stirOn", "Should stirring torque be applied",
01782                  &stirOn, FALSE);
01783    opts.optional("stirOn", "stirFilename", "PDB file with flags for "
01784      "stirred atoms (default is the PDB input file)",
01785                  PARSE_STRING);
01786    opts.optional("stirOn", "stirredAtomsCol", "Column in the stirredAtomsFile "
01787                  "containing the flags (nonzero means fixed);\n"
01788                  "default is 'O'", PARSE_STRING);
01789    opts.require("stirOn", "stirStartingTheta", "Stir starting theta offset", &stirStartingTheta);
01790    opts.require("stirOn", "stirK", "Stir force harmonic spring constant", &stirK);
01791    //should make this optional, compute from firsttimestep * stirVel
01792    opts.require("stirOn", "stirVel", "Stir angular velocity (deg/timestep)", &stirVel);
01793    opts.require("stirOn", "stirAxis", "Stir axis (direction vector)", &stirAxis);
01794    opts.require("stirOn", "stirPivot", "Stir pivot point (coordinate)", &stirPivot);
01795 
01797    opts.optionalB("main", "extraBonds",
01798                 "Should extra bonded forces be applied",
01799                  &extraBondsOn, FALSE);
01800    opts.optional("extraBonds", "extraBondsFile",
01801                 "file with list of extra bonds",
01802                  PARSE_MULTIPLES);
01803 
01804 }
01805 
01806 void SimParameters::config_parser_misc(ParseOptions &opts) {
01807    
01809    opts.optional("main", "ldBalancer", "Load balancer",
01810      loadBalancer);
01811    opts.optional("main", "ldbStrategy", "Load balancing strategy",
01812      loadStrategy);
01813    opts.optional("main", "ldbPeriod", "steps between load balancing", 
01814      &ldbPeriod);
01815    opts.range("ldbPeriod", POSITIVE);
01816    opts.optional("main", "firstLdbStep", "when to start load balancing",
01817      &firstLdbStep);
01818    opts.range("firstLdbStep", POSITIVE);
01819    opts.optional("main", "lastLdbStep", "when to stop load balancing",
01820      &lastLdbStep);
01821    opts.range("lastLdbStep", POSITIVE);
01822    opts.optional("main", "hybridGroupSize", "Hybrid load balancing group size",
01823      &hybridGroupSize);
01824    opts.optional("main", "ldbBackgroundScaling",
01825      "background load scaling", &ldbBackgroundScaling);
01826    opts.range("ldbBackgroundScaling", NOT_NEGATIVE);
01827    opts.optional("main", "ldbPMEBackgroundScaling",
01828      "PME node background load scaling", &ldbPMEBackgroundScaling);
01829    opts.range("ldbPMEBackgroundScaling", NOT_NEGATIVE);
01830    opts.optional("main", "ldbHomeBackgroundScaling",
01831      "home node background load scaling", &ldbHomeBackgroundScaling);
01832    opts.range("ldbHomeBackgroundScaling", NOT_NEGATIVE);
01833    opts.optional("main", "ldbRelativeGrainsize",
01834      "fraction of average load per compute", &ldbRelativeGrainsize, 0.);
01835    opts.range("ldbRelativeGrainsize", NOT_NEGATIVE);
01836    
01837    opts.optional("main", "traceStartStep", "when to start tracing", &traceStartStep);
01838    opts.range("traceStartStep", POSITIVE);
01839    opts.optional("main", "numTraceSteps", "the number of timesteps to be traced", &numTraceSteps);
01840    opts.range("numTraceSteps", POSITIVE);
01841  
01842 #ifdef MEASURE_NAMD_WITH_PAPI
01843    opts.optionalB("main", "papiMeasure", "whether use PAPI to measure performacne", &papiMeasure, FALSE);
01844    opts.optional("main", "papiMeasureStartStep", "when to measure performacne using PAPI", &papiMeasureStartStep);
01845    opts.range("papiMeasureStartStep", POSITIVE);
01846    opts.optional("main", "numPapiMeasureSteps", "the number of timesteps to be measured using PAPI", &numPapiMeasureSteps);
01847    opts.range("numPapiMeasureSteps", POSITIVE);
01848 #endif
01849 
01850    opts.optionalB("main", "outputMaps", "whether to dump compute map and patch map for analysis just before load balancing", &outputMaps, FALSE);
01851    opts.optionalB("main", "benchTimestep", "whether to do benchmarking timestep in which case final file output is disabled", &benchTimestep, FALSE);
01852    opts.optional("main", "useCkLoop", "whether to use CkLoop library to parallelize a loop in a function like OpenMP", &useCkLoop, 0);
01853    opts.range("useCkLoop", NOT_NEGATIVE);
01854 
01855    opts.optionalB("main", "simulateInitialMapping", "whether to study the initial mapping scheme", &simulateInitialMapping, FALSE);
01856    opts.optional("main", "simulatedPEs", "the number of PEs to be used for studying initial mapping", &simulatedPEs);
01857    opts.range("simulatedPEs", POSITIVE);
01858    opts.optional("main", "simulatedNodeSize", "the node size to be used for studying initial mapping", &simulatedNodeSize);
01859    opts.range("simulatedNodeSize", POSITIVE);
01860 
01861    opts.optionalB("main", "ldbUnloadPME", "no load on PME nodes",
01862      &ldbUnloadPME, FALSE);
01863    opts.optionalB("main", "ldbUnloadZero", "no load on pe zero",
01864      &ldbUnloadZero, FALSE);
01865    opts.optionalB("main", "ldbUnloadOne", "no load on pe one",
01866      &ldbUnloadOne, FALSE);
01867    opts.optionalB("main", "ldbUnloadOutputPEs", "no load on output PEs",
01868      &ldbUnloadOutputPEs, FALSE);
01869    opts.optionalB("main", "noPatchesOnZero", "no patches on pe zero",
01870      &noPatchesOnZero, FALSE);
01871    opts.optionalB("main", "noPatchesOnOutputPEs", "no patches on Output PEs",
01872      &noPatchesOnOutputPEs, FALSE);
01873    opts.optionalB("main", "noPatchesOnOne", "no patches on pe one",
01874      &noPatchesOnOne, FALSE);
01875    opts.optionalB("main", "useCompressedPsf", "The structure file psf is in the compressed format",
01876                   &useCompressedPsf, FALSE);
01877    opts.optionalB("main", "genCompressedPsf", "Generate the compressed version of the psf file",
01878                   &genCompressedPsf, FALSE);
01879    opts.optionalB("main", "usePluginIO", "Use the plugin I/O to load the molecule system", 
01880                   &usePluginIO, FALSE);   
01881    opts.optionalB("main", "mallocTest", "test how much memory all PEs can allocate", 
01882                   &mallocTest, FALSE);   
01883    opts.optionalB("main", "printExclusions", "print exclusion lists to stdout", 
01884                   &printExclusions, FALSE);   
01885    opts.optional("main", "proxySendSpanningTree", "using spanning tree to send proxies",
01886                   &proxySendSpanningTree, -1);
01887    opts.optional("main", "proxyRecvSpanningTree", "using spanning tree to receive proxies",
01888                   &proxyRecvSpanningTree, 0);  // default off due to memory leak -1);
01889    opts.optional("main", "proxyTreeBranchFactor", "the branch factor when building a spanning tree",
01890                   &proxyTreeBranchFactor, 4);
01891    opts.optionalB("main", "twoAwayX", "half-size patches in 1st dimension",
01892      &twoAwayX, -1);
01893    opts.optionalB("main", "twoAwayY", "half-size patches in 2nd dimension",
01894      &twoAwayY, -1);
01895    opts.optionalB("main", "twoAwayZ", "half-size patches in 3rd dimension",
01896      &twoAwayZ, -1);
01897    opts.optional("main", "maxPatches", "maximum patch count", &maxPatches, -1);
01898 
01900    opts.optional("main", "firsttimestep", "Timestep to start simulation at",
01901      &firstTimestep, 0);
01902    opts.range("firsttimestep", NOT_NEGATIVE);
01903  
01905    opts.optionalB("main", "test", "Perform self-tests rather than simulation",
01906                 &testOn, FALSE);
01907    opts.optionalB("main", "commOnly", "Do not evaluate forces or integrate",
01908                 &commOnly, FALSE);
01909 
01910    opts.optionalB("main", "statsOn", "counters in machine layer",
01911                 &statsOn, FALSE);
01913    opts.optionalB("main", "hbonds", "Use explicit hydrogen bond term",
01914                  &HydrogenBonds, FALSE);
01915    opts.optionalB("hbonds","hbAntecedents","Include Antecedent in hbond term",
01916                  &useAntecedent, TRUE);
01917    opts.optional("hbonds","hbAAexp","Hbond AA-A-H angle cos exponential",
01918                  &aaAngleExp, 2);
01919    opts.optional("hbonds","hbHAexp","Hbond D-H-A angle cos exponential",
01920                  &haAngleExp, 4);
01921    opts.optional("hbonds","hbDistAexp","Hbond A-D dist attractive exponential",
01922                  &distAttExp, 4);
01923    opts.optional("hbonds","hbDistRexp","Hbond A-D dist repulstive exponential",
01924                  &distRepExp, 6);
01925    opts.optional("hbonds","hbCutoffAngle","Hbond D-H-A cutoff angle",
01926                  &dhaCutoffAngle, 100.0);
01927    opts.range("hbCutoffAngle", NOT_NEGATIVE);
01928    opts.optional("hbonds","hbOnAngle","Hbond D-H-A switch function on angle",
01929                  &dhaOnAngle, 60.0);
01930    opts.range("hbOnAngle", NOT_NEGATIVE);
01931    opts.optional("hbonds","hbOffAngle","Hbond D-H-A switch function off angle",
01932                  &dhaOffAngle, 80.0);
01933    opts.range("hbOffAngle", NOT_NEGATIVE);
01934    opts.optional("hbonds","hbCutoffDist","Hbond A-D cutoff distance",
01935                  &daCutoffDist, 7.5);
01936    opts.range("hbCutoffDist", POSITIVE);
01937    opts.units("hbCutoffDist", N_ANGSTROM);
01938    opts.optional("hbonds","hbOnDist","Hbond A-D switch function on distance",
01939                  &daOnDist, 5.5);
01940    opts.range("hbOnDist", POSITIVE);
01941    opts.units("hbOnDist", N_ANGSTROM);
01942    opts.optional("hbonds","hbOffDist","Hbond A-D switch function off distance",
01943                  &daOffDist, 6.5);
01944    opts.range("hbOffDist", POSITIVE);
01945    opts.units("hbOffDist", N_ANGSTROM);
01946 
01947    // IMD options
01948    opts.optionalB("main","IMDon","Connect using IMD?",&IMDon, FALSE);
01949    opts.require("IMDon","IMDport", "Port to which to bind", &IMDport);
01950    opts.range("IMDport",POSITIVE);
01951    opts.require("IMDon","IMDfreq", "Frequency at which to report", &IMDfreq);
01952    opts.range("IMDfreq",POSITIVE);
01953    opts.optionalB("IMDon","IMDwait","Pause until IMD connection?",&IMDwait,
01954      FALSE);
01955    opts.optionalB("IMDon","IMDignore","Ignore forces, etc.?",&IMDignore,
01956      FALSE);
01957 
01958    // Maximum Partition options
01959    opts.optional("ldBalancer", "maxSelfPart", 
01960      "maximum number of self partitions in one patch", &maxSelfPart, 20);
01961    opts.range("maxSelfPart",POSITIVE);
01962    opts.optional("ldBalancer", "maxPairPart", 
01963      "maximum number of pair partitions in one patch", &maxPairPart, 8);
01964    opts.range("maxPairPart",POSITIVE);
01965    opts.optional("ldBalancer", "numAtomsSelf", 
01966                  "maximum number of atoms in one self compute distribution", 
01967                  &numAtomsSelf, 154);
01968    opts.range("numAtomsSelf",NOT_NEGATIVE);
01969 
01970    opts.optional("ldBalancer", "numAtomsSelf2", 
01971                  "maximum number of atoms in one self compute distribution", 
01972                  &numAtomsSelf2, 154);
01973    opts.range("numAtomsSelf2",NOT_NEGATIVE);
01974 
01975    opts.optional("ldBalancer", "numAtomsPair", 
01976                  "maximum number of atoms in one pair compute distribution", 
01977                  &numAtomsPair, 318);
01978    opts.range("numAtomsPair",NOT_NEGATIVE);
01979    opts.optional("ldBalancer", "numAtomsPair2", 
01980                "maximum number of atoms in one pair compute distribution", 
01981                &numAtomsPair2, 637);
01982    opts.range("numAtomsPair2",NOT_NEGATIVE);
01983    opts.optional("main", "minAtomsPerPatch", 
01984                "minimum average atoms per patch", 
01985                &minAtomsPerPatch, 40);
01986    opts.range("minAtomsPerPatch",NOT_NEGATIVE);
01987 
01988    // Maximum exclusion flags per atom
01989    opts.optional("main", "maxExclusionFlags", 
01990      "maximum number of exclusion flags per atom", &maxExclusionFlags, 256);
01991    opts.range("maxExclusionFlags",POSITIVE);
01992 
01993 }
01994 
01995 #ifdef MEM_OPT_VERSION
01996 //This global var is defined in mainfunc.C
01997 extern char *gWorkDir;
01998 #endif
01999 
02000 void SimParameters::check_config(ParseOptions &opts, ConfigList *config, char *&cwd) {
02001    
02002    int len;    //  String length
02003    StringList *current; //  Pointer to config option list
02004 
02005 #ifdef MEM_OPT_VERSION
02006    char *namdWorkDir = NULL;
02007 #endif
02008 
02009   if ( opts.defined("obsolete") ) {
02010     iout << iWARN <<
02011       "\"obsolete\" defined, silently ignoring obsolete options\n" << endi;
02012   }
02013 
02014    //  Take care of cwd processing
02015    if (opts.defined("cwd"))
02016    {
02017     //  First allocate and get the cwd value
02018     current = config->find("cwd");
02019 
02020     len = strlen(current->data);
02021 
02022     if ( CHDIR(current->data) )
02023     {
02024       NAMD_die("chdir() to given cwd failed!");
02025     } else {
02026       iout << iINFO << "Changed directory to " << current->data << "\n" << endi;
02027     }
02028 
02029     if (current->data[len-1] != PATHSEP)
02030       len++;
02031 
02032     cwd = new char[len+1];
02033 
02034     strcpy(cwd, current->data);
02035 
02036     if (current->data[strlen(current->data)-1] != PATHSEP)
02037       strcat(cwd, PATHSEPSTR);
02038    }
02039 
02040 #ifdef MEM_OPT_VERSION
02041    if(cwd!=NULL)namdWorkDir = cwd;     
02042    else namdWorkDir = gWorkDir;
02043    int dirlen = strlen(namdWorkDir);
02044    //only support the path representation on UNIX-like platforms
02045    char *tmpDir;
02046    if(namdWorkDir[dirlen-1]=='/'){
02047      tmpDir = new char[dirlen+1];
02048      tmpDir[dirlen] = 0;
02049    }else{
02050      tmpDir = new char[dirlen+2];
02051      tmpDir[dirlen]='/';
02052      tmpDir[dirlen+1]=0;
02053    } 
02054    memcpy(tmpDir, namdWorkDir, dirlen);
02055    namdWorkDir = tmpDir;
02056  //finished recording the per atom files, free the space for gWorkDir
02057    delete [] gWorkDir;
02058 #endif
02059 
02060 
02061    // Don't try to specify coordinates with pluginIO
02062    if ( usePluginIO && opts.defined("coordinates") ) {
02063      NAMD_die("Separate coordinates file not allowed with plugin IO, coordinates will be taken from structure file.");
02064    }
02065 
02066    // If it's not AMBER||GROMACS, then "coordinates", "structure"
02067    // and "parameters" must be specified.
02068    if (!amberOn && !gromacsOn) {
02069 #ifndef MEM_OPT_VERSION
02070      if (useCompressedPsf)
02071        NAMD_die("useCompressedPsf requires memory-optimized build!");
02072      if (!usePluginIO && !genCompressedPsf && !opts.defined("coordinates"))
02073        NAMD_die("coordinates not found in the configuration file!");
02074 #else
02075      if(!usePluginIO && !opts.defined("bincoordinates")) {
02076        NAMD_die("bincoordinates not found in the configuration file for the memory optimized version!");
02077      }
02078      if(!usePluginIO && opts.defined("coordinates")) {
02079        NAMD_die("coordinates not allowed in the configuration file for the memory optimized version!");
02080      }
02081 #endif
02082      if (!opts.defined("structure"))
02083        NAMD_die("structure not found in the configuration file!");
02084      if (!opts.defined("parameters"))
02085        NAMD_die("parameters not found in the configuration file!");
02086    }
02087    
02088    // In any case, there should be either "coordinates" or
02089    // "ambercoor", but not both
02090    if (opts.defined("coordinates") && opts.defined("ambercoor"))
02091      NAMD_die("Cannot specify both coordinates and ambercoor!");
02092 #ifndef MEM_OPT_VERSION
02093    if (!genCompressedPsf && !opts.defined("coordinates") && !opts.defined("ambercoor")
02094        && !opts.defined("grocoorfile") && !usePluginIO)
02095      NAMD_die("Coordinate file not found!");
02096 #endif
02097 
02098    //  Make sure that both a temperature and a velocity PDB were
02099    //  specified
02100    if (opts.defined("temperature") &&
02101        (opts.defined("velocities") || opts.defined("binvelocities")) ) 
02102    {
02103       NAMD_die("Cannot specify both an initial temperature and a velocity file");
02104    }
02105 
02106 #ifdef MEM_OPT_VERSION
02107 //record the absolute file name for binAtomFile, binCoorFile and binVelFile etc.
02108    binAtomFile = NULL;
02109    binCoorFile = NULL;
02110    binVelFile = NULL;   
02111    binRefFile = NULL;
02112 
02113    char *curfile = NULL;
02114    dirlen = strlen(namdWorkDir);
02115    current = config->find("structure");;
02116    curfile = current->data;
02117    int filelen = strlen(curfile);
02118    if(*curfile == '/' || *curfile=='~') {
02119      //check whether it is an absolute path
02120      //WARNING: Only works on Unix-like platforms!
02121      //Needs to fix on Windows platform.
02122      //-Chao Mei     
02123      //adding 5 because of ".bin"+"\0"
02124      binAtomFile = new char[filelen+5];
02125      memcpy(binAtomFile, curfile, filelen);
02126      memcpy(binAtomFile+filelen, ".bin", 4);
02127      binAtomFile[filelen+4] = 0;
02128    }else{
02129      binAtomFile = new char[dirlen+filelen+5];
02130      memcpy(binAtomFile, namdWorkDir, dirlen);
02131      memcpy(binAtomFile+dirlen, curfile, filelen);
02132      memcpy(binAtomFile+dirlen+filelen, ".bin", 4);
02133      binAtomFile[dirlen+filelen+4] = 0;
02134    }
02135 
02136    current = config->find("bincoordinates");
02137    curfile = current->data;
02138    filelen = strlen(curfile);
02139    if(*curfile == '/' || *curfile=='~') {
02140      binCoorFile = new char[filelen+1];
02141      memcpy(binCoorFile, curfile, filelen);
02142      binCoorFile[filelen] = 0;
02143    }else{
02144      binCoorFile = new char[dirlen+filelen+1];
02145      memcpy(binCoorFile, namdWorkDir, dirlen);
02146      memcpy(binCoorFile+dirlen, curfile, filelen);
02147      binCoorFile[dirlen+filelen] = 0;
02148    }
02149 
02150    if(opts.defined("binvelocities")){
02151      current = config->find("binvelocities");
02152      curfile = current->data;
02153      filelen = strlen(curfile);
02154      if(*curfile == '/' || *curfile=='~') {
02155        binVelFile = new char[filelen+1];
02156        memcpy(binVelFile, curfile, filelen);
02157        binVelFile[filelen] = 0;
02158      }else{
02159        binVelFile = new char[dirlen+filelen+1];
02160        memcpy(binVelFile, namdWorkDir, dirlen);
02161        memcpy(binVelFile+dirlen, curfile, filelen);
02162        binVelFile[dirlen+filelen] = 0;
02163      }
02164    }
02165 
02166    if(opts.defined("binrefcoords")){
02167      current = config->find("binrefcoords");
02168      curfile = current->data;
02169      filelen = strlen(curfile);
02170      if(*curfile == '/' || *curfile=='~') {
02171        binRefFile = new char[filelen+1];
02172        memcpy(binRefFile, curfile, filelen);
02173        binRefFile[filelen] = 0;
02174      }else{
02175        binRefFile = new char[dirlen+filelen+1];
02176        memcpy(binRefFile, namdWorkDir, dirlen);
02177        memcpy(binRefFile+dirlen, curfile, filelen);
02178        binRefFile[dirlen+filelen] = 0;
02179      }
02180    }
02181 
02182    //deal with output file name to make it absolute path for parallel output
02183    if(outputFilename[0] != '/' && outputFilename[0]!='~') {
02184      filelen = strlen(outputFilename);
02185      char *tmpout = new char[filelen];
02186      memcpy(tmpout, outputFilename, filelen);
02187      CmiAssert(filelen+dirlen <= 120); //leave 8 chars for file suffix
02188      memcpy(outputFilename, namdWorkDir, dirlen);
02189      memcpy(outputFilename+dirlen, tmpout, filelen);
02190      outputFilename[filelen+dirlen] = 0;     
02191      delete [] tmpout;
02192    }
02193 
02194    if ( dcdFrequency && opts.defined("dcdfile") &&
02195         dcdFilename[0] != '/' && dcdFilename[0]!='~' ) {
02196      filelen = strlen(dcdFilename);
02197      char *tmpout = new char[filelen];
02198      memcpy(tmpout, dcdFilename, filelen);
02199      CmiAssert(filelen+dirlen <= 120); //leave 8 chars for file suffix
02200      memcpy(dcdFilename, namdWorkDir, dirlen);
02201      memcpy(dcdFilename+dirlen, tmpout, filelen);
02202      dcdFilename[filelen+dirlen] = 0;     
02203      delete [] tmpout;
02204    }
02205 
02206    if ( velDcdFrequency && opts.defined("veldcdfile") &&
02207         velDcdFilename[0] != '/' && velDcdFilename[0]!='~' ) {
02208      filelen = strlen(velDcdFilename);
02209      char *tmpout = new char[filelen];
02210      memcpy(tmpout, velDcdFilename, filelen);
02211      CmiAssert(filelen+dirlen <= 120); //leave 8 chars for file suffix
02212      memcpy(velDcdFilename, namdWorkDir, dirlen);
02213      memcpy(velDcdFilename+dirlen, tmpout, filelen);
02214      velDcdFilename[filelen+dirlen] = 0;     
02215      delete [] tmpout;
02216    }
02217 
02218    if ( forceDcdFrequency && opts.defined("forcedcdfile") &&
02219         forceDcdFilename[0] != '/' && forceDcdFilename[0]!='~' ) {
02220      filelen = strlen(forceDcdFilename);
02221      char *tmpout = new char[filelen];
02222      memcpy(tmpout, forceDcdFilename, filelen);
02223      CmiAssert(filelen+dirlen <= 120); //leave 8 chars for file suffix
02224      memcpy(forceDcdFilename, namdWorkDir, dirlen);
02225      memcpy(forceDcdFilename+dirlen, tmpout, filelen);
02226      forceDcdFilename[filelen+dirlen] = 0;     
02227      delete [] tmpout;
02228    }
02229 
02230    if ( restartFrequency && opts.defined("restartname") &&
02231         restartFilename[0] != '/' && restartFilename[0]!='~' ) {
02232      filelen = strlen(restartFilename);
02233      char *tmpout = new char[filelen];
02234      memcpy(tmpout, restartFilename, filelen);
02235      CmiAssert(filelen+dirlen <= 120); //leave 8 chars for file suffix
02236      memcpy(restartFilename, namdWorkDir, dirlen);
02237      memcpy(restartFilename+dirlen, tmpout, filelen);
02238      restartFilename[filelen+dirlen] = 0;     
02239      delete [] tmpout;
02240    }
02241 
02242    delete [] namdWorkDir;
02243 
02244    if (opts.defined("numinputprocs")) { 
02245      if(numinputprocs > CkNumPes()) {
02246        iout << iWARN << "The number of input processors exceeds the total number of processors. Resetting to half of the number of total processors.\n" << endi;
02247        numinputprocs = (CkNumPes()>>1)+(CkNumPes()&1);
02248      }
02249    }
02250 
02251    if (opts.defined("numoutputprocs")) {        
02252      if(numoutputprocs > CkNumPes()) {
02253        iout << iWARN << "The number of output processors exceeds the total number of processors. Resetting to half of the number of total processors.\n" << endi;
02254        numoutputprocs = (CkNumPes()>>1)+(CkNumPes()&1);
02255      }
02256    }
02257 
02258    #if !OUTPUT_SINGLE_FILE
02259    //create directories for multi-file output scheme   
02260    create_output_directories("coor");
02261    create_output_directories("vel");
02262    if(dcdFrequency) {
02263            create_output_directories("dcd");
02264            if(opts.defined("dcdfile")){
02265                    iout << iWARN << "The dcd file output has been changed to directory: " << outputFilename << ".\n" << endi; 
02266            }
02267    }
02268    if (velDcdFrequency) {
02269            create_output_directories("veldcd");
02270            if(opts.defined("veldcdfile")){       
02271                    iout << iWARN << "The veldcd file output has been changed to directory: " << outputFilename << ".\n" << endi;
02272            }
02273    }
02274    if (forceDcdFrequency) {
02275            create_output_directories("forcedcd");
02276            if(opts.defined("forcedcdfile")){       
02277                    iout << iWARN << "The forcedcd file output has been changed to directory: " << outputFilename << ".\n" << endi;
02278            }
02279    }
02280    #endif
02281 #endif
02282 
02283    if (! opts.defined("auxFile")) {
02284      strcpy(auxFilename,outputFilename);
02285      strcat(auxFilename,".aux");
02286    }
02287 
02288    //  Check for frequencies
02289    if (dcdFrequency) {
02290      if (! opts.defined("dcdfile")) {
02291        strcpy(dcdFilename,outputFilename);
02292        strcat(dcdFilename,".dcd");
02293      }
02294    } else {
02295      dcdFilename[0] = STRINGNULL;
02296    }
02297 
02298    if (velDcdFrequency) {
02299      if (! opts.defined("veldcdfile")) {
02300        strcpy(velDcdFilename,outputFilename);
02301        strcat(velDcdFilename,".veldcd");
02302      }
02303    } else {
02304      velDcdFilename[0] = STRINGNULL;
02305    }
02306    
02307    if (forceDcdFrequency) {
02308      if (! opts.defined("forcedcdfile")) {
02309        strcpy(forceDcdFilename,outputFilename);
02310        strcat(forceDcdFilename,".forcedcd");
02311      }
02312    } else {
02313      forceDcdFilename[0] = STRINGNULL;
02314    }
02315    
02316    if (xstFrequency) {
02317      if (! opts.defined("xstfile")) {
02318        strcpy(xstFilename,outputFilename);
02319        strcat(xstFilename,".xst");
02320      }
02321    } else {
02322      xstFilename[0] = STRINGNULL;
02323    }
02324 
02325    if (restartFrequency) {
02326      if (! opts.defined("restartname")) {
02327        strcpy(restartFilename,outputFilename);
02328        strcat(restartFilename,".restart");
02329      }
02330    } else {
02331      restartFilename[0] = STRINGNULL;
02332      restartSave = FALSE;
02333      binaryRestart = FALSE;
02334    }
02335 
02336    if (storeComputeMap || loadComputeMap) {
02337      if (! opts.defined("computeMapFile")) {
02338        strcpy(computeMapFilename,"computeMapFile");
02339        strcat(computeMapFilename,".txt");
02340      }
02341    }
02342 
02343 
02344    if (!amberOn)
02345    { //****** BEGIN CHARMM/XPLOR type changes
02347      if (!paraTypeXplorOn && !paraTypeCharmmOn) 
02348      {
02349        paraTypeXplorOn = TRUE;
02350      }
02352      if (paraTypeXplorOn && paraTypeCharmmOn) 
02353      {
02354        NAMD_die("Please specify either XPLOR or CHARMM format for parameters!");
02355      }
02356      //****** END CHARMM/XPLOR type changes
02357    }
02358 
02359    
02360    //  If minimization isn't on, must have a temp or velocity
02361    if (!(minimizeOn||minimizeCGOn) && !opts.defined("temperature") && 
02362        !opts.defined("velocities") && !opts.defined("binvelocities") ) 
02363    {
02364       NAMD_die("Must have either an initial temperature or a velocity file");
02365    }
02366 
02367    if (minimizeOn||minimizeCGOn) { initialTemp = 0.0; }
02368    if (opts.defined("velocities") || opts.defined("binvelocities") )
02369    {
02370   initialTemp = -1.0;
02371    }
02372 
02374 
02375    if ( opts.defined("extendedSystem") )
02376    {
02377      current = config->find("extendedSystem");
02378 
02379      iout << iINFO << "EXTENDED SYSTEM FILE   "
02380         << current->data << "\n" << endi;
02381 
02382      ifstream xscFile(current->data);
02383      if ( ! xscFile ) NAMD_die("Unable to open extended system file.\n");
02384 
02385      char labels[1024];
02386      do {
02387        if ( ! xscFile ) NAMD_die("Error reading extended system file.\n");
02388        xscFile.getline(labels,1023);
02389      } while ( strncmp(labels,"#$LABELS ",9) );
02390 
02391      int a_x, a_y, a_z, b_x, b_y, b_z, c_x, c_y, c_z;
02392      a_x = a_y = a_z = b_x = b_y = b_z = c_x = c_y = c_z = -1;
02393      int o_x, o_y, o_z, s_u, s_v, s_w, s_x, s_y, s_z;
02394      o_x = o_y = o_z = s_u = s_v = s_w = s_x = s_y = s_z = -1;
02395 
02396      int pos = 0;
02397      char *l_i = labels + 8;
02398      while ( *l_i ) {
02399        if ( *l_i == ' ' ) { ++l_i; continue; }
02400        char *l_i2;
02401        for ( l_i2 = l_i; *l_i2 && *l_i2 != ' '; ++l_i2 );
02402        if ( (l_i2 - l_i) == 3 && (l_i[1] == '_') ) {
02403          if (l_i[0] == 'a' && l_i[2] == 'x') a_x = pos;
02404          if (l_i[0] == 'a' && l_i[2] == 'y') a_y = pos;
02405          if (l_i[0] == 'a' && l_i[2] == 'z') a_z = pos;
02406          if (l_i[0] == 'b' && l_i[2] == 'x') b_x = pos;
02407          if (l_i[0] == 'b' && l_i[2] == 'y') b_y = pos;
02408          if (l_i[0] == 'b' && l_i[2] == 'z') b_z = pos;
02409          if (l_i[0] == 'c' && l_i[2] == 'x') c_x = pos;
02410          if (l_i[0] == 'c' && l_i[2] == 'y') c_y = pos;
02411          if (l_i[0] == 'c' && l_i[2] == 'z') c_z = pos;
02412          if (l_i[0] == 'o' && l_i[2] == 'x') o_x = pos;
02413          if (l_i[0] == 'o' && l_i[2] == 'y') o_y = pos;
02414          if (l_i[0] == 'o' && l_i[2] == 'z') o_z = pos;
02415          if (l_i[0] == 's' && l_i[2] == 'u') s_u = pos;
02416          if (l_i[0] == 's' && l_i[2] == 'v') s_v = pos;
02417          if (l_i[0] == 's' && l_i[2] == 'w') s_w = pos;
02418          if (l_i[0] == 's' && l_i[2] == 'x') s_x = pos;
02419          if (l_i[0] == 's' && l_i[2] == 'y') s_y = pos;
02420          if (l_i[0] == 's' && l_i[2] == 'z') s_z = pos;
02421        }
02422        ++pos;
02423        l_i = l_i2;
02424      }
02425      int numpos = pos;
02426 
02427      for ( pos = 0; pos < numpos; ++pos ) {
02428        double tmp;
02429        xscFile >> tmp;
02430        if ( ! xscFile ) NAMD_die("Error reading extended system file.\n");
02431        if ( pos == a_x ) cellBasisVector1.x = tmp;
02432        if ( pos == a_y ) cellBasisVector1.y = tmp;
02433        if ( pos == a_z ) cellBasisVector1.z = tmp;
02434        if ( pos == b_x ) cellBasisVector2.x = tmp;
02435        if ( pos == b_y ) cellBasisVector2.y = tmp;
02436        if ( pos == b_z ) cellBasisVector2.z = tmp;
02437        if ( pos == c_x ) cellBasisVector3.x = tmp;
02438        if ( pos == c_y ) cellBasisVector3.y = tmp;
02439        if ( pos == c_z ) cellBasisVector3.z = tmp;
02440        if ( pos == o_x ) cellOrigin.x = tmp;
02441        if ( pos == o_y ) cellOrigin.y = tmp;
02442        if ( pos == o_z ) cellOrigin.z = tmp;
02443        if ( pos == s_u ) strainRate2.x = tmp;
02444        if ( pos == s_v ) strainRate2.y = tmp;
02445        if ( pos == s_w ) strainRate2.z = tmp;
02446        if ( pos == s_x ) strainRate.x = tmp;
02447        if ( pos == s_y ) strainRate.y = tmp;
02448        if ( pos == s_z ) strainRate.z = tmp;
02449      }
02450 
02451    }
02452 
02453    if ( LJcorrection && ! cellBasisVector3.length2() ) {
02454      NAMD_die("Can't use LJ tail corrections without periodic boundary conditions!");
02455    }
02456 
02457    if ( cellBasisVector3.length2() && ! cellBasisVector2.length2() ) {
02458      NAMD_die("Used cellBasisVector3 without cellBasisVector2!");
02459    }
02460 
02461    if ( cellBasisVector2.length2() && ! cellBasisVector1.length2() ) {
02462      NAMD_die("Used cellBasisVector2 without cellBasisVector1!");
02463    }
02464 
02465    if ( cellOrigin.length2() && ! cellBasisVector1.length2() ) {
02466      NAMD_die("Used cellOrigin without cellBasisVector1!");
02467    }
02468 
02469    lattice.set(cellBasisVector1,cellBasisVector2,cellBasisVector3,cellOrigin);
02470 
02471    if (! opts.defined("DCDunitcell")) {
02472       dcdUnitCell = lattice.a_p() && lattice.b_p() && lattice.c_p();
02473    }
02474 
02475    char s[129];
02476 
02478    if ( ! opts.defined("cylindricalBCAxis") )
02479    {
02480       cylindricalBCAxis = 'x';
02481    }
02482    else
02483    {
02484      opts.get("cylindricalBCAxis", s);
02485 
02486      if (!strcasecmp(s, "x"))
02487      {
02488       cylindricalBCAxis = 'x';
02489      }
02490      else if (!strcasecmp(s, "y"))
02491      {
02492       cylindricalBCAxis = 'y';
02493      }
02494      else if (!strcasecmp(s, "z"))
02495      {
02496       cylindricalBCAxis = 'z';
02497      }
02498    else
02499      {
02500       char err_msg[128];
02501 
02502       sprintf(err_msg, "Illegal value '%s' for 'cylindricalBCAxis' in configuration file", s);
02503       NAMD_die(err_msg);
02504      }
02505    }
02506 
02507    if (!opts.defined("splitPatch"))
02508    {
02509      splitPatch = SPLIT_PATCH_HYDROGEN;
02510    }
02511    else
02512    {
02513      opts.get("splitPatch", s);
02514      if (!strcasecmp(s, "position"))
02515        splitPatch = SPLIT_PATCH_POSITION;
02516      else if (!strcasecmp(s,"hydrogen"))
02517        splitPatch = SPLIT_PATCH_HYDROGEN;
02518      else
02519      {
02520        char err_msg[129];
02521        sprintf(err_msg, 
02522           "Illegal value '%s' for 'splitPatch' in configuration file", 
02523        s);
02524        NAMD_die(err_msg);
02525      }
02526    }
02527 
02529    opts.get("exclude", s);
02530 
02531    if (!strcasecmp(s, "none"))
02532    {
02533       exclude = NONE;
02534       splitPatch = SPLIT_PATCH_POSITION;
02535    }
02536    else if (!strcasecmp(s, "1-2"))
02537    {
02538       exclude = ONETWO;
02539       splitPatch = SPLIT_PATCH_POSITION;
02540    }
02541    else if (!strcasecmp(s, "1-3"))
02542    {
02543       exclude = ONETHREE;
02544    }
02545    else if (!strcasecmp(s, "1-4"))
02546    {
02547       exclude = ONEFOUR;
02548    }
02549    else if (!strcasecmp(s, "scaled1-4"))
02550    {
02551       exclude = SCALED14;
02552    }
02553    else
02554    {
02555       char err_msg[128];
02556 
02557       sprintf(err_msg, "Illegal value '%s' for 'exclude' in configuration file",
02558    s);
02559       NAMD_die(err_msg);
02560    }
02561 
02562    if (scale14 != 1.0 && exclude != SCALED14)
02563    {
02564       iout << iWARN << "Exclude is not scaled1-4; 1-4scaling ignored.\n" << endi;
02565    }
02566 
02567    // water model stuff
02568    if (!opts.defined("waterModel")) {
02569      watmodel = WAT_TIP3;
02570    } else {
02571      opts.get("waterModel", s);
02572      if (!strncasecmp(s, "tip4", 4)) {
02573        iout << iINFO << "Using TIP4P water model.\n" << endi;
02574        watmodel = WAT_TIP4;
02575      } else if (!strncasecmp(s, "tip3", 4)) {
02576        iout << iINFO << "Using TIP3P water model.\n" << endi;
02577      } else if (!strncasecmp(s, "swm4", 4)) {
02578        iout << iINFO << "Using SWM4-DP water model.\n" << endi;
02579      } else {
02580        char err_msg[128];
02581        sprintf(err_msg,
02582            "Illegal value %s for 'waterModel' in configuration file", s);
02583        NAMD_die(err_msg);
02584      }
02585    }
02586    if (watmodel == WAT_SWM4 && !drudeOn) {
02587      NAMD_die("Must have 'drudeOn' enabled to use SWM4-DP water model.");
02588    }
02589    if (drudeOn && watmodel != WAT_SWM4) {
02590      watmodel = WAT_SWM4;
02591      iout << iWARN
02592        << "Setting water model to 'swm4' (SWM4-DP) for Drude polarization.\n"
02593        << endi;
02594    }
02595 
02596    // Drude water model uses "lonepairs"
02597    if (watmodel == WAT_SWM4) {
02598      lonepairs = TRUE;
02599    }
02600 
02601    // Added by JLai -- 8.2.11 -- Checks if Go method is defined
02602    if (goForcesOn) {
02603      iout << iINFO << "Go forces are on\n" << endi;
02604    // Added by JLai -- 6.3.11 -- Checks if Go method is defined
02605      int * gomethod = &goMethod;
02606      if (!opts.defined("GoMethod")) {
02607        *gomethod = 0;
02608        //     printf("GO METHOD IS NOT DEFINED SO WE'LL SET IT TO SOME WEIRD VALUE\n");
02609      } else {
02610        opts.get("GoMethod",s);
02611        // printf("GO METHOD IS DEFINED SO WE'LL PRINT IT OUT: %s\n",s);
02612        *gomethod = atoi(s);
02613      }
02614      if (!strcasecmp(s, "matrix")) {
02615        goMethod = 1;
02616        //GoMethod = GO_MATRIX;
02617      } else if (!strcasecmp(s, "faster")) {
02618        goMethod = 2;
02619        //GoMethod = GO_FASTER;
02620      } else if (!strcasecmp(s, "lowmem")) {
02621        goMethod = 3;
02622        //GoMethod = GO_LOWMEM;
02623      }
02624      else {
02625        char err_msg[129];     
02626        sprintf(err_msg,
02627                "Illegal value '%s' for 'GoMethod' in configuration file",
02628                s);
02629        NAMD_die(err_msg);
02630      }
02631    }  // End of NAMD code to check goMethod
02632    // End of Port -- JL
02633 
02634    //  Get multiple timestep integration scheme
02635    if (!opts.defined("MTSAlgorithm"))
02636    {
02637   MTSAlgorithm = VERLETI;
02638    }
02639    else
02640    {
02641   opts.get("MTSAlgorithm", s);
02642 
02643   if (!strcasecmp(s, "naive"))
02644   {
02645     MTSAlgorithm = NAIVE;
02646   }
02647   else if (!strcasecmp(s, "constant"))
02648   {
02649     MTSAlgorithm = NAIVE;
02650   }
02651   else if (!strcasecmp(s, "impulse"))
02652   {
02653     MTSAlgorithm = VERLETI;
02654   }
02655   else if (!strcasecmp(s, "verleti"))
02656   {
02657     MTSAlgorithm = VERLETI;
02658   }
02659   else
02660   {
02661     char err_msg[129];
02662 
02663     sprintf(err_msg, 
02664        "Illegal value '%s' for 'MTSAlgorithm' in configuration file", 
02665        s);
02666     NAMD_die(err_msg);
02667   }
02668    }
02669 
02670    //  Get the long range force splitting specification
02671    if (!opts.defined("longSplitting"))
02672    {
02673   longSplitting = C1;
02674    }
02675    else
02676    {
02677   opts.get("longSplitting", s);
02678   if (!strcasecmp(s, "sharp"))
02679     longSplitting = SHARP;
02680   else if (!strcasecmp(s, "xplor"))
02681     longSplitting = XPLOR;
02682   else if (!strcasecmp(s, "c1"))
02683     longSplitting = C1;
02684   else if (!strcasecmp(s, "c2"))
02685     longSplitting = C2;
02686   else
02687   {
02688     char err_msg[129];
02689 
02690     sprintf(err_msg, 
02691        "Illegal value '%s' for 'longSplitting' in configuration file", 
02692        s);
02693     NAMD_die(err_msg);
02694   }
02695    }
02696 
02697    // take care of rigid bond options
02698    if (!opts.defined("rigidBonds"))
02699    {
02700       rigidBonds = RIGID_NONE;
02701    }
02702    else
02703    {
02704       opts.get("rigidBonds", s); 
02705       if (!strcasecmp(s, "all"))
02706       {
02707           rigidBonds = RIGID_ALL;
02708       }
02709       else if (!strcasecmp(s, "water"))
02710       {
02711            rigidBonds = RIGID_WATER;
02712       } 
02713       else if (!strcasecmp(s, "none"))
02714       {
02715            rigidBonds = RIGID_NONE;
02716       } 
02717       else
02718       {
02719         char err_msg[256];
02720         sprintf(err_msg, 
02721           "Illegal value '%s' for 'rigidBonds' in configuration file", s);
02722         NAMD_die(err_msg);
02723       }
02724    }
02725    
02726    //  Take care of switching stuff
02727    if (switchingActive)
02728    {
02729 
02730      if (!opts.defined("switchDist")) {
02731        NAMD_die("switchDist must be defined when switching is enabled");
02732      }
02733 
02734      if ( (switchingDist>cutoff) || (switchingDist<0) )
02735      {
02736        char err_msg[129];
02737 
02738        sprintf(err_msg, 
02739          "switchDist muct be between 0 and cutoff, which is %f", cutoff);
02740        NAMD_die(err_msg);
02741      }
02742 
02743    }
02744 
02745    if ( martiniSwitching )
02746    {
02747      if ( ! switchingActive ) 
02748      { 
02749        NAMD_die("martiniSwitching requires switching");
02750      }
02751      if ( vdwForceSwitching ) 
02752      { 
02753        NAMD_die("martiniSwitching and vdwForceSwitching are exclusive to one another. Select only one."); 
02754      }
02755      if ( dielectric != 15.0 && ! martiniDielAllow ) 
02756      {
02757        iout << iWARN << "USE DIELECTRIC OF 15.0 WITH MARTINI.\n";
02758        iout << iWARN << "SETTING dielectric 15.0\n";
02759        iout << iWARN << "FOR NON-STANDARD DIELECTRIC WITH MARTINI, SET: martiniDielAllow on\n";
02760        dielectric = 15.0;
02761      }
02762      if ( ! cosAngles )
02763      {
02764        iout << iWARN << "USE COSINE BASED ANGLES WITH MARTINI.\n";
02765        iout << iWARN << "SETTING cosAngles on\n";
02766        cosAngles = TRUE;
02767      }
02768      if ( PMEOn )
02769      {
02770        NAMD_die("Do not use Particle Mesh Ewald with Martini.  Set: PME off");
02771      }
02772      if ( MSMOn )
02773      {
02774        NAMD_die("Do not use Multilevel Summation Method with Martini.  Set: MSM off");
02775      }
02776      if ( FMMOn )
02777      {
02778        NAMD_die("Do not use Fast Multipole Method with Martini.  Set: FMM off");
02779      }
02780 
02781    }
02782 
02783 
02784    if (!opts.defined("pairlistDist"))
02785    {
02786   pairlistDist = cutoff;
02787    }
02788    else if (pairlistDist < cutoff)
02789    {
02790   NAMD_die("pairlistDist must be >= cutoff distance");
02791    }
02792 
02793    patchDimension = pairlistDist;
02794 
02795    if ( splitPatch == SPLIT_PATCH_HYDROGEN ) {
02796      patchDimension += hgroupCutoff;
02797    }
02798 
02799    BigReal defaultMargin = 0.0;
02800    if (berendsenPressureOn || langevinPistonOn) {
02801       defaultMargin = ( useFlexibleCell ? 0.06 : 0.03 ) * patchDimension;
02802    }
02803    if ( margin == XXXBIGREAL ) {
02804      margin = defaultMargin;
02805    }
02806    if ( defaultMargin != 0.0 && margin == 0.0 ) {
02807      margin = defaultMargin;
02808      iout << iWARN << "ALWAYS USE NON-ZERO MARGIN WITH CONSTANT PRESSURE!\n";
02809      iout << iWARN << "CHANGING MARGIN FROM 0 to " << margin << "\n" << endi;
02810    }
02811 
02812    patchDimension += margin;
02813 
02814     //ensure patch can handle alpha_cutoff for gbis
02815     if (GBISOn) {
02816       //Check compatibility
02817       if (fullDirectOn) {
02818         NAMD_die("GBIS not compatible with FullDirect");
02819       }
02820       if (PMEOn) {
02821         NAMD_die("GBIS not compatible with PME");
02822       }
02823       if (MSMOn) {
02824         NAMD_die("GBIS not compatible with MSM");
02825       }
02826       if (FMMOn) {
02827         NAMD_die("GBIS not compatible with FMM");
02828       }
02829       if (alchOn) {
02830         NAMD_die("GBIS not compatible with Alchemical Transformations");
02831       }
02832       if (lesOn) {
02833         NAMD_die("GBIS not compatible with Locally Enhanced Sampling");
02834       }
02835       if (FMAOn) {
02836         NAMD_die("GBIS not compatible with FMA");
02837       }
02838       if (drudeOn) {
02839         NAMD_die("GBIS not compatible with Drude Polarization");
02840       }
02841       if (fixedAtomsOn) {
02842 #ifdef NAMD_CUDA
02843         NAMD_die("GBIS CUDA not yet compatible with fixed atoms");
02844 #endif
02845       }
02846 
02847       if (alpha_cutoff > patchDimension) {
02848         patchDimension = alpha_cutoff; 
02849       }
02850       //calculate kappa
02851       BigReal tmp = (initialTemp > 0) ? initialTemp : 300;
02852       kappa = 50.29216*sqrt(ion_concentration/solvent_dielectric/tmp);
02853       /*magic number = 1/sqrt(eps0*kB/(2*nA*e^2*1000))*/
02854     } // GBISOn
02855 
02856     if (LCPOOn) {
02857 #ifdef MEM_OPT_VERSION
02858       NAMD_die("SASA not yet available for memory optimized builds");
02859 #endif
02860       if ( lattice.volume() > 0 ) {
02861         NAMD_die("SASA does not yet support periodic boundary conditions.");
02862       }
02863       //LCPO requires patches to be at least 16.2Ang in each dimension
02864       // twoAway[XYZ} is ignored for now
02865     }
02866 
02867    //  Turn on global integration if not explicitly specified
02868 
02869    if ( dihedralOn ) globalOn = TRUE;
02870 
02871 #ifdef NAMD_CUDA
02872    if (loweAndersenOn) {
02873        NAMD_die("Lowe-Andersen dynamics not compatible with CUDA at this time");
02874    }
02875 #endif
02876    // END LA
02877 
02878    // BEGIN LA
02879    if (loweAndersenOn && (langevinOn || tCoupleOn))
02880    {
02881       NAMD_die("Lowe-Andersen dynamics, Langevin dynamics and temperature coupling are mutually exclusive dynamics modes");
02882    }
02883    // END LA
02884 
02885    if (tCoupleOn && opts.defined("rescaleFreq") )
02886    {
02887       NAMD_die("Temperature coupling and temperature rescaling are mutually exclusive");
02888    }
02889 
02890    if (globalOn && CkNumPes() > 1)
02891    {
02892       NAMD_die("Global integration does not run in parallel (yet).");
02893    }
02894 
02895    if (COLDOn && langevinOn)
02896    {
02897       NAMD_die("COLD and Langevin dynamics are mutually exclusive dynamics modes");
02898    }
02899    if (COLDOn && minimizeOn)
02900    {
02901       NAMD_die("COLD and minimization are mutually exclusive dynamics modes");
02902    }
02903    if (COLDOn && tCoupleOn)
02904    {
02905       NAMD_die("COLD and temperature coupling are mutually exclusive dynamics modes");
02906    }
02907    if (COLDOn && opts.defined("rescaleFreq"))
02908    {
02909       NAMD_die("COLD and velocity rescaling are mutually exclusive dynamics modes");
02910    }
02911 
02912    if (splitPatch == SPLIT_PATCH_POSITION && mollyOn )
02913    {
02914       NAMD_die("splitPatch hydrogen is required for MOLLY");
02915    }
02916 
02917    if (splitPatch == SPLIT_PATCH_POSITION && rigidBonds != RIGID_NONE)
02918    {
02919       NAMD_die("splitPatch hydrogen is required for rigidBonds");
02920    }
02921 
02922    //  Set the default value for the maximum movement parameter
02923    //  for minimization
02924    if (minimizeOn && (maximumMove == 0.0)) 
02925    {
02926       maximumMove = 0.75 * pairlistDist/stepsPerCycle;
02927    }
02928    if (adaptTempOn) {
02929      if (!adaptTempRescale && !adaptTempLangevin) 
02930         NAMD_die("Adaptive tempering needs to be coupled to either the Langevin thermostat or velocity rescaling.");
02931      if (opts.defined("adaptTempInFile") && (opts.defined("adaptTempTmin") ||
02932                                              opts.defined("adaptTempTmax") ||
02933                                              adaptTempBins != 0)) 
02934         NAMD_die("cannot simultaneously specify adaptTempInFile and any of {adaptTempTmin, adaptTempTmax,adaptTempBins} as these are read from the input file");
02935      if (!opts.defined("adaptTempInFile") && !(opts.defined("adaptTempTmin") &&
02936                                              opts.defined("adaptTempTmax") &&
02937                                              adaptTempBins != 0 ))  
02938         NAMD_die("Need to specify either adaptTempInFile or all of {adaptTempTmin, adaptTempTmax,adaptTempBins} if adaptTempMD is on.");
02939    }
02940    if (langevinOn) {
02941      if ( ! opts.defined("langevinDamping") ) langevinDamping = 0.0;
02942      if ( ! opts.defined("langevinHydrogen") ) langevinHydrogen = TRUE;
02943      if ( (opts.defined("langevinDamping") || opts.defined("langevinHydrogen"))
02944        && (opts.defined("langevinFile") || opts.defined("langevinCol")) )
02945        NAMD_die("To specify Langevin dynamics parameters, use either langevinDamping and langevinHydrogen or langevinFile and langevinCol.  Do not combine them.");
02946      if ( opts.defined("langevinHydrogen") && langevinDamping == 0.0 )
02947        NAMD_die("langevinHydrogen requires langevinDamping to be set.");
02948    }
02949    
02950    // BEGIN LA
02951    if (loweAndersenOn) {
02952        if (!opts.defined("loweAndersenRate")) loweAndersenRate = 100;
02953        if (!opts.defined("loweAndersenCutoff")) loweAndersenCutoff = 2.7;
02954    }
02955    // END LA
02956 
02957    if (opts.defined("rescaleFreq"))
02958    {
02959   if (!opts.defined("rescaleTemp"))
02960   {
02961     if (opts.defined("temperature"))
02962     {
02963       rescaleTemp = initialTemp;
02964     }
02965     else
02966     {
02967       NAMD_die("Must give a rescale temperature if rescaleFreq is defined");
02968     }
02969   }
02970    }
02971    else
02972    {
02973   rescaleFreq = -1;
02974   rescaleTemp = 0.0;
02975    }
02976 
02977    if (opts.defined("rescaleTemp"))
02978    {
02979   if (!opts.defined("rescaleFreq"))
02980   {
02981     NAMD_die("Must give a rescale freqency if rescaleTemp is given");
02982   }
02983    }
02984 
02985    if (opts.defined("reassignFreq"))
02986    {
02987   if (!opts.defined("reassignTemp"))
02988   {
02989     if (opts.defined("temperature"))
02990     {
02991       reassignTemp = initialTemp;
02992     }
02993     else
02994     {
02995       NAMD_die("Must give a reassign temperature if reassignFreq is defined");
02996     }
02997   }
02998    }
02999    else
03000    {
03001   reassignFreq = -1;
03002   reassignTemp = 0.0;
03003    }
03004 
03005    if (opts.defined("reassignTemp"))
03006    {
03007   if (!opts.defined("reassignFreq"))
03008   {
03009     NAMD_die("Must give a reassignment freqency if reassignTemp is given");
03010   }
03011    }
03012 
03013    if (opts.defined("reassignIncr"))
03014    {
03015   if (!opts.defined("reassignFreq"))
03016   {
03017     NAMD_die("Must give a reassignment freqency if reassignIncr is given");
03018   }
03019    }
03020    else
03021    {
03022   reassignIncr = 0.0;
03023    }
03024 
03025    if (opts.defined("reassignHold"))
03026    {
03027   if (!opts.defined("reassignIncr"))
03028   {
03029     NAMD_die("Must give a reassignment increment if reassignHold is given");
03030   }
03031    }
03032    else
03033    {
03034   reassignHold = 0.0;
03035    }
03036 
03037    if (!opts.defined("seed")) 
03038    {
03039       randomSeed = (unsigned int) time(NULL);
03040    }
03041 
03042 //Modifications for alchemical fep
03043 
03044    alchFepOn = FALSE;
03045    alchThermIntOn = FALSE;
03046 
03047    if (alchOn) {
03048 
03049      if (!opts.defined("alchType")) 
03050      {
03051        NAMD_die("Must define type of alchemical simulation: fep or ti\n");
03052      }
03053      else
03054      {
03055        opts.get("alchType",s);
03056        if (!strcasecmp(s, "fep"))
03057        {
03058          alchFepOn = TRUE;
03059        }
03060        else if (!strcasecmp(s, "ti"))
03061        {
03062          alchThermIntOn = TRUE;
03063        }
03064        else
03065        {
03066          NAMD_die("Unknown type of alchemical simulation; choices are fep or ti\n");
03067        }
03068      }
03069 
03070      if      (rescaleFreq > 0)  alchTemp = rescaleTemp;
03071      else if (reassignFreq > 0) alchTemp = reassignTemp;
03072      else if (langevinOn)       alchTemp = langevinTemp;
03073      else if (tCoupleOn)        alchTemp = tCoupleTemp;
03074      else NAMD_die("Alchemical FEP can be performed only in constant temperature simulations\n");
03075 
03076      if (reassignFreq > 0 && reassignIncr != 0)
03077         NAMD_die("reassignIncr cannot be used in alchemical simulations\n");
03078 
03079      if (alchLambda < 0.0 || alchLambda > 1.0 || alchLambda2 < 0.0 || alchLambda2 > 1.0)
03080         NAMD_die("Alchemical lambda values should be in the range [0.0, 1.0]\n");
03081 
03082      if ( alchOn && alchVdwLambdaEnd > 1.0) 
03083         NAMD_die("Gosh tiny Elvis, you kicked soft-core in the van der Waals! alchVdwLambdaEnd should be in the range [0.0, 1.0]\n");
03084 
03085      if ( alchOn && alchElecLambdaStart > 1.0) 
03086         NAMD_die("alchElecLambdaStart should be in the range [0.0, 1.0]\n");
03087      
03088      if (alchFepOn)
03089      {
03090        if (!opts.defined("alchoutfile")) {
03091        strcpy(alchOutFile, outputFilename);
03092        strcat(alchOutFile, ".fep");
03093        }
03094 
03095            if( (!alchFepWhamOn) && ( (!opts.defined("alchLambda")) || (!opts.defined("alchLambda2"))) ) {
03096                    NAMD_die("alchFepOn is on, but alchLambda or alchLambda2 is not set.");
03097              }
03098        
03099        if(alchRepLambda > 1.0)  NAMD_die("alchRepLambda should be in the range [0.0, 1.0].");
03100        else if(alchRepLambda >= 0.0)    alchFepWCARepuOn = true;
03101        else     alchFepWCARepuOn = false;
03102  
03103        if(alchDispLambda > 1.0) NAMD_die("alchDispLambda should be in the range [0.0, 1.0].");
03104        else if(alchDispLambda >= 0.0)   alchFepWCADispOn = true;
03105        else     alchFepWCADispOn = false;
03106  
03107        if(alchElecLambda > 1.0) NAMD_die("alchElecLambda should be in the range [0.0, 1.0].");
03108        else if(alchElecLambda >= 0.0)   alchFepElecOn = true;
03109        else     alchFepElecOn = false;
03110        
03111        if( (alchFepWCARepuOn || alchFepWCADispOn || alchFepElecOn) && (!alchFepWhamOn))
03112           NAMD_die("alchFepWhamOn has to be on if one of alchFepWCARepuOn/alchFepWCADispOn/alchFepElecOn is set.");
03113        if (alchFepWCARepuOn && alchFepWCADispOn)
03114           NAMD_die("With WCA decomposition, repulsion and dispersion can NOT be in the same FEP stage");
03115        if (alchFepWCARepuOn && alchFepElecOn)
03116           NAMD_die("With WCA decomposition, repulsion and electrostatic perturbation can NOT be in the same FEP stage");
03117        if (alchFepWCADispOn && alchFepElecOn)
03118           NAMD_die("With WCA decomposition, dispersion and electrostatic perturbation can NOT be in the same FEP stage");
03119        if (alchFepWCARepuOn && (!opts.defined("alchFepWCArcut1")||!opts.defined("alchFepWCArcut2")||!opts.defined("alchFepWCArcut3")))
03120           NAMD_die("When using WCA repulsion,  alchFepWCArcut1, alchFepWCArcut2, and alchFepWCArcut3 must be defined!");
03121        if (alchFepWCARepuOn && ((alchFepWCArcut1 > alchFepWCArcut2) || (alchFepWCArcut2 > alchFepWCArcut3) ))
03122            NAMD_die("When using WCA repulsion,  alchFepWCArcut2 must be larger than alchFEPWCArcut1, alchFepWCArcut3 must be larger than alchFEPWCArcut2!");
03123 //       if ((alchFepWCARepuOn || alchFepWCADispOn) && (alchElecLambdaStart < 1.0) )
03124 //           NAMD_die("When using WCA decomposition,  repulsion, dispersion and electrostatic must be in 3 different stages!");
03125        if(alchFepWhamOn && (alchRepLambda < 0.0) && (alchDispLambda < 0.0) && (alchElecLambda < 0.0) )  
03126            NAMD_die("One of alchRepLambda, alchDispLambda and alchElecLambda should be set up when alchFepWhamOn is true!");
03127        if(alchFepWhamOn && (!alchFepElecOn) )   {
03128          alchElecLambda = 0.0;
03129          ComputeNonbondedUtil::alchElecLambda = alchElecLambda;
03130 //               ComputeNonbondedUtil::select();
03131        }
03132      }
03133      else if (alchThermIntOn)
03134      {
03135        if (!opts.defined("alchoutfile")) { 
03136          strcpy(alchOutFile, outputFilename); 
03137          strcat(alchOutFile, ".ti"); 
03138        } 
03139      }
03140 
03141    } else {
03142      alchLambda = alchLambda2 = 0;
03143      alchElecLambdaStart = 0;
03144      alchOutFile[0] = STRINGNULL;
03145    }
03146 
03147 
03148 //fepe
03149 
03150    if ( alchOn && alchFepOn && alchThermIntOn )
03151      NAMD_die("Sorry, combined TI and FEP is not implemented.\n");
03152    if ( alchOn && lesOn )
03153      NAMD_die("Sorry, combined LES with FEP or TI is not implemented.\n");
03154    if ( alchOn && alchThermIntOn && lesOn )
03155      NAMD_die("Sorry, combined LES and TI is not implemented.\n");
03156    if ( alchDecouple && (! (alchFepOn || alchThermIntOn) ) ) {
03157          iout << iWARN << "Alchemical decoupling was requested but \
03158            alchemical free energy calculation is not active. Setting \
03159            alchDecouple to off.\n" << endi;
03160          alchDecouple = FALSE;
03161      //NAMD_die("Alchemcial decoupling was requested but alchemical free energy calculation is not active.\n");
03162    }
03163 
03164    if ( lesOn && ( lesFactor < 1 || lesFactor > 255 ) ) {
03165      NAMD_die("lesFactor must be positive and less than 256");
03166    }
03167    if ((pairInteractionOn && alchFepOn) || (pairInteractionOn && lesOn) || (pairInteractionOn && alchThermIntOn) ) 
03168      NAMD_die("Sorry, pair interactions may not be calculated when LES, FEP or TI is enabled.");
03169 
03170    // Drude model
03171    if (drudeOn) {
03172      if ( ! langevinOn ) {
03173        NAMD_die("Drude model requires use of Langevin thermostat.");
03174      }
03175      if ( ! opts.defined("drudeDamping")) {
03176        drudeDamping = langevinDamping;
03177        iout << iWARN << "Undefined 'drudeDamping' will be set to "
03178          "value of 'langevinDamping'\n" << endi;
03179      }
03180      if ( ! opts.defined("drudeBondConst")) {
03181        drudeBondConst = 0;
03182        if (opts.defined("drudeBondLen")) {
03183          iout << iWARN << "Resetting 'drudeBondLen' to 0 "
03184            "since 'drudeBondConst' is unset\n" << endi;
03185        }
03186        drudeBondLen = 0;
03187      }
03188      if ( ! opts.defined("drudeNbtholeCut")) {
03189        drudeNbtholeCut = 0;
03190          iout << iWARN << "Resetting 'drudeNbtholeCut' to 0 "
03191            "since 'drudeNbtholeCut' is unset\n" << endi;
03192      }
03193    }
03194 
03195    //  Set up load balancing variables
03196    if (opts.defined("ldBalancer")) {
03197      if (strcasecmp(loadBalancer, "none") == 0)
03198        ldBalancer = LDBAL_NONE;
03199      else if (strcasecmp(loadBalancer, "hybrid") == 0)
03200        ldBalancer = LDBAL_HYBRID;
03201      else
03202        NAMD_die("Unknown ldBalancer selected");
03203    } else {
03204      ldBalancer = LDBAL_CENTRALIZED;
03205    }
03206 
03207    if (opts.defined("ldbStrategy")) {
03208      //  Assign the load balancing strategy
03209      if (strcasecmp(loadStrategy, "comprehensive") == 0)
03210        ldbStrategy = LDBSTRAT_COMPREHENSIVE;
03211      else if (strcasecmp(loadStrategy, "refineonly") == 0)
03212        ldbStrategy = LDBSTRAT_REFINEONLY;
03213      else if (strcasecmp(loadStrategy, "old") == 0)
03214        ldbStrategy = LDBSTRAT_OLD;
03215      else
03216        NAMD_die("Unknown ldbStrategy selected");
03217    } else {
03218      ldbStrategy = LDBSTRAT_DEFAULT;
03219    }
03220 
03221   if (!opts.defined("ldbPeriod")) {
03222     ldbPeriod=200*stepsPerCycle;
03223   }
03224 
03225   //  Set default values
03226   if (!opts.defined("firstLdbStep")) {
03227     firstLdbStep=5*stepsPerCycle;
03228   }
03229 
03230   if (ldbPeriod <= firstLdbStep) {
03231     NAMD_die("ldbPeriod must greater than firstLdbStep.");
03232   }
03233 
03234   if (!opts.defined("lastLdbStep")) {
03235     lastLdbStep = -1;
03236   }
03237 
03238   if (!opts.defined("hybridGroupSize")) {
03239     hybridGroupSize = 512;
03240   }
03241 
03242   // tracing will be done if trace is available and user says +traceOff
03243   // in that case we set nice values for some functions
03244   bool specialTracing = traceAvailable() && (traceIsOn() == 0);
03245 
03246   if(!opts.defined("traceStartStep")) {
03247     traceStartStep = 4 * firstLdbStep + 2 * ldbPeriod;
03248   }
03249   if(!opts.defined("numTraceSteps")) {
03250     numTraceSteps = 100;
03251   }
03252 
03253   if(specialTracing) {
03254     if (!opts.defined("firstLdbStep")) firstLdbStep = 20;
03255     if (!opts.defined("ldbPeriod")) ldbPeriod = 100;
03256 
03257     if(!opts.defined("traceStartStep")) {
03258       traceStartStep = 4 * firstLdbStep + 2 * ldbPeriod; // 380
03259     }
03260 
03261     if(!opts.defined("numTraceSteps")) {
03262       numTraceSteps = 80;
03263     }
03264   }
03265 
03266 #ifdef MEASURE_NAMD_WITH_PAPI
03267   if(papiMeasure){
03268           if(!opts.defined("papiMeasureStartStep")) {
03269                   papiMeasureStartStep = 3 * firstLdbStep;
03270           }
03271           if(!opts.defined("numPapiMeasureSteps")) {
03272                   numPapiMeasureSteps = 8; //including two pme steps
03273           }
03274   }
03275 #endif
03276 
03277   if(simulateInitialMapping) {
03278           if(!opts.defined("simulatedPEs")){
03279                   simulatedPEs = CkNumPes();
03280           }
03281           if(!opts.defined("simulatedNodeSize")){
03282                   simulatedNodeSize = CkMyNodeSize();
03283           }
03284   }
03285 
03286 #ifdef MEM_OPT_VERSION
03287   //Some constraints on the values of load balancing parameters.
03288   //The reason is related to communication schemes used in sending proxy
03289   //data. If the step immediately after the load balancing is not a step
03290   //for atom migration, then it's possible there are some necessary information 
03291   // missing inside the ProxyPatch which will crash the program. Therefore,
03292   // It's better that the step immediately after the load balancing be a step
03293   // for atom migration so that the some overhead in Proxy msgs are removed. 
03294   // --Chao Mei
03295   if(ldbPeriod%stepsPerCycle!=0 || firstLdbStep%stepsPerCycle!=0) {
03296       iout << iWARN << "In memory optimized version, the ldbPeriod parameter or firstLdbStep parameter is better set to be a multiple of stepsPerCycle parameter!\n";
03297   }
03298 #endif
03299 
03300    if (N < firstTimestep) { N = firstTimestep; }
03301 
03302    if ( (firstTimestep%stepsPerCycle) != 0)
03303    {
03304   NAMD_die("First timestep must be a multiple of stepsPerCycle!!");
03305    }
03306 
03307    //  Make sure only one full electrostatics algorithm is selected
03308    {
03309      int i = 0;
03310      if ( FMAOn ) ++i;
03311      if ( PMEOn ) ++i;
03312      if ( MSMOn ) ++i;
03313      if ( FMMOn ) ++i;
03314      if ( fullDirectOn ) ++i;
03315      if ( i > 1 )
03316         NAMD_die("More than one full electrostatics algorithm selected!!!");
03317    }
03318 
03319    if (!opts.defined("ldbBackgroundScaling")) {
03320      ldbBackgroundScaling = 1.0;
03321    }
03322    if (!opts.defined("ldbPMEBackgroundScaling")) {
03323      ldbPMEBackgroundScaling = ldbBackgroundScaling;
03324    }
03325    if (!opts.defined("ldbHomeBackgroundScaling")) {
03326      ldbHomeBackgroundScaling = ldbBackgroundScaling;
03327    }
03328 
03329    //  Check on PME parameters
03330    if (PMEOn) {  // idiot checking
03331      if ( lattice.volume() == 0. ) {
03332         NAMD_die("PME requires periodic boundary conditions.");
03333      }
03334      if ( PMEGridSpacing == 0. ) {
03335        if ( PMEGridSizeX * PMEGridSizeY * PMEGridSizeZ == 0 )
03336          NAMD_die("Either PMEGridSpacing or PMEGridSizeX, PMEGridSizeY, and PMEGridSizeZ must be specified.");
03337        else PMEGridSpacing = 1.5;  // only exit in very bad cases
03338      }
03339 #ifndef TEST_PME_GRID
03340      for ( int idim = 0; idim < 3; ++idim ) {
03341         int *gridSize;
03342         BigReal cellLength;
03343         const char *direction;
03344         switch ( idim ) {
03345         case 0:  direction = "X";
03346            gridSize = &PMEGridSizeX;  cellLength = lattice.a().length();
03347            break;
03348         case 1:  direction = "Y";
03349            gridSize = &PMEGridSizeY;  cellLength = lattice.b().length();
03350            break;
03351         case 2:  direction = "Z";
03352            gridSize = &PMEGridSizeZ;  cellLength = lattice.c().length();
03353            break;
03354         }
03355         int minSize = (int) ceil(cellLength/PMEGridSpacing);
03356 #else
03357   for ( int minSize = 1; minSize < 300; ++minSize ) {
03358 #endif
03359         int bestSize = 10 * (minSize + 10);  // make sure it's big
03360         int max2, max3, ts;
03361         for ( max2=2, ts=1; ts < minSize; ++max2 ) ts *= 2;
03362         for ( max3=2, ts=1; ts < minSize; ++max3 ) ts *= 3;
03363         int max5 = 2;
03364         int max7 = 1;
03365         int max11 = 1;
03366         for ( int i2 = 0; i2 <= max2; ++i2 ) {
03367         for ( int i3 = 0; i3 <= max3; ++i3 ) {
03368         for ( int i5 = 0; i5 <= max5; ++i5 ) {
03369         for ( int i7 = 0; i7 <= max7; ++i7 ) {
03370         for ( int i11 = 0; i11 <= max11; ++i11 ) {
03371            if ( i5 + i7 + i11 > i2 ) continue;
03372            int testSize = 2;  // must be even
03373            for ( int j2 = 0; j2 < i2; ++j2 ) testSize *= 2;
03374            if ( testSize > bestSize ) continue;
03375            for ( int j3 = 0; j3 < i3; ++j3 ) testSize *= 3;
03376            if ( testSize > bestSize ) continue;
03377            for ( int j5 = 0; j5 < i5; ++j5 ) testSize *= 5;
03378            if ( testSize > bestSize ) continue;
03379            for ( int j7 = 0; j7 < i7; ++j7 ) testSize *= 7;
03380            if ( testSize > bestSize ) continue;
03381            for ( int j11 = 0; j11 < i11; ++j11 ) testSize *= 11;
03382            if ( testSize > bestSize ) continue;
03383            if ( testSize >= minSize ) bestSize = testSize;
03384         } } } } }
03385 #ifdef TEST_PME_GRID
03386   iout << minSize << " " << bestSize << "\n" << endi;
03387 #else
03388         if ( ! *gridSize ) {   // set it
03389            *gridSize = bestSize;
03390         }
03391         if ( *gridSize * PMEGridSpacing < cellLength ) {
03392            char errmsg[512];
03393            sprintf(errmsg, "PMEGridSize%s %d is too small for cell length %f and PMEGridSpacing %f\n",
03394                 direction, *gridSize, cellLength, PMEGridSpacing);
03395            NAMD_die(errmsg);
03396         }
03397 #endif
03398      }
03399      if ( PMEGridSizeX < 5 ) {
03400         NAMD_die("PMEGridSizeX (number of grid points) is very small.");
03401      }
03402      if ( PMEGridSizeY < 5 ) {
03403         NAMD_die("PMEGridSizeY (number of grid points) is very small.");
03404      }
03405      if ( PMEGridSizeZ < 5 ) {
03406         NAMD_die("PMEGridSizeZ (number of grid points) is very small.");
03407      }
03408      BigReal tolerance = PMETolerance;
03409      BigReal ewaldcof = 1.0;
03410      while ( erfc(ewaldcof*cutoff)/cutoff >= tolerance ) ewaldcof *= 2.0;
03411      BigReal ewaldcof_lo = 0.;
03412      BigReal ewaldcof_hi = ewaldcof;
03413      for ( int i = 0; i < 100; ++i ) {
03414        ewaldcof = 0.5 * ( ewaldcof_lo + ewaldcof_hi );
03415        if ( erfc(ewaldcof*cutoff)/cutoff >= tolerance ) {
03416          ewaldcof_lo = ewaldcof;
03417        } else {
03418          ewaldcof_hi = ewaldcof;
03419        }
03420      }
03421      PMEEwaldCoefficient = ewaldcof;
03422    } else {  // initialize anyway
03423      useDPME = 0;
03424      PMEGridSizeX = 0;
03425      PMEGridSizeY = 0;
03426      PMEGridSizeZ = 0;
03427      PMEGridSpacing = 1000.;
03428      PMEEwaldCoefficient = 0;
03429    }
03430 
03431 
03432    //  Take care of initializing FMA values to something if FMA is not
03433    //  active
03434    if (!FMAOn)
03435    {
03436      FMALevels = 0;
03437      FMAMp = 0;
03438      FMAFFTOn = FALSE;
03439      FMAFFTBlock = 0;
03440    }
03441    else
03442    {
03443   // idiot checking: frm bug reported by Tom Bishop.
03444   // DPMTA requires: (#terms in fma)/(fft blocking factor) = integer.
03445   if (FMAFFTBlock != 4)
03446     NAMD_die("FMAFFTBlock: Block length must be 4 for short FFT's");
03447   if (FMAMp % FMAFFTBlock != 0)
03448     NAMD_die("FMAMp: multipole term must be multiple of block length (FMAFFTBlock)");
03449     }
03450 
03451    if ( (nonbondedFrequency > stepsPerCycle) || ( (stepsPerCycle % nonbondedFrequency) != 0) )
03452    {
03453      NAMD_die("stepsPerCycle must be a multiple of nonbondedFreq");
03454    }
03455 
03456    if (!LCPOOn && !GBISOn && !GBISserOn && !FMAOn && !PMEOn && !MSMOn && !fullDirectOn && !FMMOn)
03457    {
03458      fullElectFrequency = 0;
03459    }
03460    else
03461    {
03462      if (!opts.defined("fullElectFrequency"))
03463      {
03464        if (opts.defined("fmaFrequency")) {
03465          iout << iWARN << "The parameter fmaFrequency has been renamed fullElectFrequency.\n" << endi;
03466          fullElectFrequency = fmaFrequency;
03467        } else {
03468          iout << iWARN << "The parameter fullElectFrequency now defaults to nonbondedFreq (" << nonbondedFrequency << ") rather than stepsPerCycle.\n" << endi;
03469          fullElectFrequency = nonbondedFrequency;
03470        }
03471      }
03472      else
03473      {
03474        if (opts.defined("fmaFrequency")) {
03475          iout << iWARN << "Ignoring redundant parameter fmaFrequency in favor of fullElectFrequency.\n" << endi;
03476        }
03477        if ( (fullElectFrequency > stepsPerCycle) || ( (stepsPerCycle % fullElectFrequency) != 0) )
03478        {
03479          NAMD_die("stepsPerCycle must be a multiple of fullElectFrequency");
03480        }
03481      }
03482 
03483      if ( (nonbondedFrequency > fullElectFrequency) || ( (fullElectFrequency % nonbondedFrequency) != 0) )
03484      {
03485        NAMD_die("fullElectFrequency must be a multiple of nonbondedFreq");
03486      }
03487 
03488 
03489      if (!opts.defined("fmaTheta"))
03490      fmaTheta=0.715;  /* Suggested by Duke developers */
03491    }
03492 
03493    if ( lesOn && ( FMAOn || useDPME || fullDirectOn ) ) {
03494      NAMD_die("Sorry, LES is only implemented for PME full electrostatics.");
03495    }
03496    if ( alchFepOn && ( FMAOn || useDPME || fullDirectOn ) ) {
03497      NAMD_die("Sorry, FEP is only implemented for PME full electrostatics.");
03498    }
03499    if ( alchThermIntOn && ( FMAOn || useDPME || fullDirectOn ) ) {
03500      NAMD_die("Sorry, TI is only implemented for PME full electrostatics.");
03501    }
03502    if ( pairInteractionOn && FMAOn ) {
03503      NAMD_die("Sorry, pairInteraction not implemented for FMA.");
03504    }
03505    if ( pairInteractionOn && useDPME ) {
03506      NAMD_die("Sorry, pairInteraction not implemented for DPME.");
03507    }
03508    if ( pairInteractionOn && fullDirectOn ) {
03509      NAMD_die("Sorry, pairInteraction not implemented for full direct electrostatics.");
03510    }
03511    if ( ! pairInteractionOn ) {
03512      pairInteractionSelf = 0;
03513    }
03514    if ( pairInteractionOn && !pairInteractionSelf && !config->find("pairInteractionGroup2")) 
03515      NAMD_die("pairInteractionGroup2 must be specified");
03516 
03517    if ( ! fixedAtomsOn ) {
03518      fixedAtomsForces = 0;
03519    }
03520 
03521    if ( gridforceOn || mgridforceOn ) {
03522      parse_mgrid_params(config);
03523    }
03524       
03525    if (!opts.defined("constraints"))
03526    {
03527      constraintExp = 0;     
03528      constraintScaling = 1.0;     
03529 
03530      //****** BEGIN selective restraints (X,Y,Z) changes
03531      selectConstraintsOn = FALSE;
03532      //****** END selective restraints (X,Y,Z) changes
03533  
03534      //****** BEGIN moving constraints changes 
03535      movingConstraintsOn = FALSE;
03536      //****** END moving constraints changes 
03537      //****** BEGIN rotating constraints changes 
03538      rotConstraintsOn = FALSE;
03539     //****** END rotating constraints changes 
03540    } 
03541    //****** BEGIN rotating constraints changes 
03542    else {
03543      if (rotConstraintsOn) {
03544        rotConsAxis = rotConsAxis.unit();
03545      }
03546    }
03547    if(opts.defined("rotConstraints") 
03548       && opts.defined("movingConstraints")) {
03549      NAMD_die("Rotating and moving constraints are mutually exclusive!");
03550    }
03551    //****** END rotating constraints changes 
03552 
03553    //****** BEGIN selective restraints (X,Y,Z) changes
03554    if(opts.defined("selectConstraints") && !opts.defined("selectConstrX")
03555       && !opts.defined("selectConstrY") && !opts.defined("selectConstrZ")) {
03556      NAMD_die("selectConstraints was specified, but no Cartesian components were defined!");
03557    }
03558    if (!opts.defined("selectConstraints")) {
03559        constrXOn = FALSE;
03560        constrYOn = FALSE;
03561        constrZOn = FALSE;
03562    }
03563    //****** END selective restraints (X,Y,Z) changes
03564 
03565 
03566    //****** BEGIN SMD constraints changes 
03567    
03568    if (!opts.defined("SMD")) {     
03569      SMDOn = FALSE;
03570    }
03571 
03572    if (SMDOn) {
03573      // normalize direction
03574      if (SMDDir.length2() == 0) {
03575        NAMD_die("SMD direction vector must be non-zero");
03576      }
03577      else {
03578        SMDDir = SMDDir.unit();
03579      }
03580 
03581      if (SMDOutputFreq > 0 && SMDOutputFreq < stepsPerCycle
03582          || SMDOutputFreq % stepsPerCycle != 0) {
03583        NAMD_die("SMDOutputFreq must be a multiple of stepsPerCycle");
03584      }
03585    }
03586      
03587    //****** END SMD constraints changes 
03588    
03589    if (!sphericalBCOn)
03590      {
03591   sphericalBCr1 = 0.0;
03592   sphericalBCk1 = 0.0;
03593   sphericalBCexp1 = 0;
03594   sphericalBCr2 = 0.0;
03595   sphericalBCk2 = 0.0;
03596   sphericalBCexp2 = 0;
03597    }
03598    else if (!opts.defined("sphericalBCr2"))
03599    {
03600       sphericalBCr2 = -1.0;
03601       sphericalBCk2 = 0.0;
03602       sphericalBCexp2 = 0;
03603    }
03604 
03605    if (!cylindricalBCOn)
03606    {
03607     cylindricalBCr1 = 0.0;
03608     cylindricalBCk1 = 0.0;
03609     cylindricalBCexp1 = 0;
03610     cylindricalBCr2 = 0.0;
03611     cylindricalBCk2 = 0.0;
03612     cylindricalBCexp2 = 0;
03613     cylindricalBCl1 = 0.0;
03614     cylindricalBCl2 = 0.0;
03615    }
03616    else if (!opts.defined("cylindricalBCr2"))
03617    {
03618     cylindricalBCr2 = -1.0;
03619     cylindricalBCk2 = 0.0;
03620     cylindricalBCexp2 = 0;
03621     cylindricalBCl2 = 0.0;
03622    }
03623 
03624    if (!eFieldOn)
03625    {
03626         eField.x = 0.0;
03627         eField.y = 0.0;
03628         eField.z = 0.0;
03629         eFieldFreq = 0.0;
03630         eFieldPhase = 0.0;
03631    }
03632    else
03633    {
03634         if (!opts.defined("eFieldFreq")) eFieldFreq = 0.0;
03635         if (!opts.defined("eFieldPhase")) eFieldPhase = 0.0;
03636    }
03637 
03638    if (!stirOn)
03639    { 
03640      stirFilename[0] = STRINGNULL;
03641      stirStartingTheta = 0.0;
03642      stirVel = 0.0;
03643      stirK = 0.0;
03644      stirAxis.x = 0.0;
03645      stirAxis.y = 0.0;
03646      stirAxis.z = 0.0;
03647      stirPivot.x = 0.0;
03648      stirPivot.y = 0.0;
03649      stirPivot.z = 0.0;
03650    }
03651   
03652    if (!opts.defined("langevin"))
03653    {
03654   langevinTemp = 0.0;
03655    }
03656    
03657    // BEGIN LA
03658    if (!opts.defined("loweAndersen"))
03659    {
03660        loweAndersenTemp = 0.0;
03661    }
03662    // END LA
03663 
03664    if (!opts.defined("tcouple"))
03665    {
03666   tCoupleTemp = 0.0;
03667    }
03668 
03669    if (HydrogenBonds)
03670    {
03671      if (daCutoffDist > pairlistDist)
03672        NAMD_die("Hydrogen bond cutoff distance must be <= pairlist distance");
03673    }
03674 
03675    // If we're doing pair interaction, set 
03676    // outputEnergies to 1 to make NAMD not die (the other nonbonded code paths 
03677    // aren't defined when these options are enabled), and set nonbondedFreq to 
03678    // 1 to avoid getting erroneous output.  Warn the user of what we're doing.
03679    if (pairInteractionOn) {
03680            if (outputEnergies != 1) {
03681                    iout << iWARN << "Setting outputEnergies to 1 due to\n";
03682                    iout << iWARN << "pairInteraction calculations\n" << endi;
03683                    outputEnergies  = 1;
03684            }
03685    }
03686    if (pairInteractionOn || pressureProfileOn) {
03687            if (nonbondedFrequency != 1) {
03688                    iout << iWARN << "Setting nonbondedFreq to 1 due to\n";
03689                    iout << iWARN << "pairInteraction or pressure profile calculations\n" << endi;
03690            }
03691    }
03692 
03693    // print timing at a reasonable interval by default
03694    if (!opts.defined("outputTiming"))
03695    {
03696       outputTiming = firstLdbStep;
03697       int ot2 = 10 * outputEnergies;
03698       if ( outputTiming < ot2 ) outputTiming = ot2;
03699    }
03700 
03701 }
03702 
03703 void SimParameters::print_config(ParseOptions &opts, ConfigList *config, char *&cwd) {
03704 
03705    StringList *current; //  Pointer to config option list
03706 
03707    //  Now that we have read everything, print it out so that
03708    //  the user knows what is going on
03709    iout << iINFO << "SIMULATION PARAMETERS:\n";
03710    iout << iINFO << "TIMESTEP               " << dt << "\n" << endi;
03711    iout << iINFO << "NUMBER OF STEPS        " << N << "\n";
03712    iout << iINFO << "STEPS PER CYCLE        " << stepsPerCycle << "\n";
03713    iout << endi;
03714 
03715    if ( lattice.a_p() || lattice.b_p() || lattice.c_p() ) {
03716      if ( lattice.a_p() )
03717        iout << iINFO << "PERIODIC CELL BASIS 1  " << lattice.a() << "\n";
03718      if ( lattice.b_p() )
03719        iout << iINFO << "PERIODIC CELL BASIS 2  " << lattice.b() << "\n";
03720      if ( lattice.c_p() )
03721        iout << iINFO << "PERIODIC CELL BASIS 3  " << lattice.c() << "\n";
03722      iout << iINFO << "PERIODIC CELL CENTER   " << lattice.origin() << "\n";
03723      if (wrapWater) {
03724        iout << iINFO << "WRAPPING WATERS AROUND PERIODIC BOUNDARIES ON OUTPUT.\n";
03725      }
03726      if (wrapAll) {
03727        iout << iINFO << "WRAPPING ALL CLUSTERS AROUND PERIODIC BOUNDARIES ON OUTPUT.\n";
03728      }
03729      if (wrapNearest) {
03730        iout << iINFO << "WRAPPING TO IMAGE NEAREST TO PERIODIC CELL CENTER.\n";
03731      }
03732      iout << endi;
03733    }
03734 
03735    if ( CkNumPes() > 512 ) ldbUnloadOne = TRUE;
03736    if ( ldbUnloadOne || CkNumPes() > 128 ) ldbUnloadZero = TRUE;
03737 
03738    if (ldBalancer == LDBAL_NONE) {
03739      iout << iINFO << "LOAD BALANCER  None\n" << endi;
03740    } else {
03741      if (ldBalancer == LDBAL_CENTRALIZED) {
03742        iout << iINFO << "LOAD BALANCER  Centralized\n" << endi;
03743      } else if (ldBalancer == LDBAL_HYBRID) {
03744        iout << iINFO << "LOAD BALANCER  Hybrid\n" << endi;
03745      }
03746 
03747      if (ldbStrategy == LDBSTRAT_DEFAULT) {
03748        iout << iINFO << "LOAD BALANCING STRATEGY  New Load Balancers -- DEFAULT\n";
03749      } else if (ldbStrategy == LDBSTRAT_REFINEONLY) {
03750        iout << iINFO << "LOAD BALANCING STRATEGY  Refinement Only\n";
03751      } else if (ldbStrategy == LDBSTRAT_COMPREHENSIVE) {
03752        iout << iINFO << "LOAD BALANCING STRATEGY  Comprehensive\n";
03753      } else if (ldbStrategy == LDBSTRAT_OLD) {
03754        iout << iINFO << "LOAD BALANCING STRATEGY  Old Load Balancers\n";
03755      }
03756 
03757      iout << iINFO << "LDB PERIOD             " << ldbPeriod << " steps\n";
03758      iout << iINFO << "FIRST LDB TIMESTEP     " << firstLdbStep << "\n";
03759      if (ldBalancer == LDBAL_HYBRID)
03760        iout << iINFO << "HYBRIDLB GROUP SIZE     " << hybridGroupSize << "\n";
03761      iout << iINFO << "LAST LDB TIMESTEP     " << lastLdbStep << "\n";
03762      if ( ldbRelativeGrainsize > 0. )
03763        iout << iINFO << "LDB RELATIVE GRAINSIZE " << ldbRelativeGrainsize << "\n";
03764      iout << iINFO << "LDB BACKGROUND SCALING " << ldbBackgroundScaling << "\n";
03765      iout << iINFO << "HOM BACKGROUND SCALING " << ldbHomeBackgroundScaling << "\n";
03766      if ( PMEOn ) {
03767        iout << iINFO << "PME BACKGROUND SCALING "
03768                                 << ldbPMEBackgroundScaling << "\n";
03769      if ( ldbUnloadPME )
03770      iout << iINFO << "REMOVING LOAD FROM PME NODES" << "\n";
03771      }
03772      if ( ldbUnloadZero ) iout << iINFO << "REMOVING LOAD FROM NODE 0\n";
03773      if ( ldbUnloadOne ) iout << iINFO << "REMOVING LOAD FROM NODE 1\n";
03774      if ( ldbUnloadOutputPEs ) iout << iINFO << "REMOVING LOAD FROM OUTPUT PES\n";
03775      iout << endi;
03776    }
03777 
03778    if ( ldbUnloadOne || CkNumPes() > 256 ) noPatchesOnOne = TRUE;
03779    if ( ldbUnloadZero || noPatchesOnOne ||
03780           CkNumPes() > 64 || ( IMDon && CkNumPes() > 8 ) ) {
03781      noPatchesOnZero = TRUE;
03782    }
03783    if ( noPatchesOnZero ) iout << iINFO << "REMOVING PATCHES FROM PROCESSOR 0\n";
03784    if ( noPatchesOnOne ) iout << iINFO << "REMOVING PATCHES FROM PROCESSOR 1\n";     
03785    iout << endi;
03786 
03787 #ifdef NAMD_CUDA
03788     maxSelfPart = maxPairPart = 1;
03789 #endif
03790 
03791    if (ldBalancer == LDBAL_HYBRID) {
03792      iout << iINFO << "MAX SELF PARTITIONS    " << maxSelfPart << "\n"
03793           << iINFO << "MAX PAIR PARTITIONS    " << maxPairPart << "\n"
03794           << iINFO << "SELF PARTITION ATOMS   " << numAtomsSelf << "\n"
03795           << iINFO << "SELF2 PARTITION ATOMS   " << numAtomsSelf2 << "\n"
03796           << iINFO << "PAIR PARTITION ATOMS   " << numAtomsPair << "\n"
03797           << iINFO << "PAIR2 PARTITION ATOMS  " << numAtomsPair2 << "\n";
03798    }
03799    iout << iINFO << "MIN ATOMS PER PATCH    " << minAtomsPerPatch << "\n"
03800         << endi;
03801    
03802    if (initialTemp < 0)
03803    {
03804   current = config->find("velocities");
03805 
03806   if (current == NULL)
03807   {
03808     current = config->find("binvelocities");
03809   }
03810 
03811   iout << iINFO << "VELOCITY FILE          " << current->data << "\n";
03812    }
03813    else
03814    {
03815   iout << iINFO << "INITIAL TEMPERATURE    " 
03816      << initialTemp << "\n";
03817    }
03818    iout << endi;
03819 
03820    iout << iINFO << "CENTER OF MASS MOVING INITIALLY? ";
03821 
03822    if (comMove)
03823    {
03824      iout << "YES\n";
03825    }
03826    else
03827    {
03828      iout << "NO\n";
03829    }
03830    iout << endi;
03831 
03832    if ( zeroMomentum ) {
03833      iout << iINFO << "REMOVING CENTER OF MASS DRIFT DURING SIMULATION";
03834      if ( zeroMomentumAlt ) iout << " (ALT METHOD)";
03835      iout << "\n" << endi;
03836    }
03837 
03838    iout << iINFO << "DIELECTRIC             " 
03839       << dielectric << "\n";
03840 
03841    if ( nonbondedScaling != 1.0 )
03842    {
03843      iout << iINFO << "NONBONDED SCALING    " << nonbondedScaling << "\n" << endi;
03844    }
03845    iout << iINFO << "EXCLUDE                ";
03846 
03847    switch (exclude)
03848    {
03849      case NONE:
03850        iout << "NONE\n";
03851        break;
03852      case ONETWO:
03853        iout << "ONETWO\n";
03854        break;
03855      case ONETHREE:
03856        iout << "ONETHREE\n";
03857        break;
03858      case ONEFOUR:
03859        iout << "ONE-FOUR\n";
03860        break;
03861      default:
03862        iout << "SCALED ONE-FOUR\n";
03863        break;
03864    }
03865    iout << endi;
03866 
03867    if (exclude == SCALED14)
03868    {
03869      iout << iINFO << "1-4 ELECTROSTATICS SCALED BY " << scale14 << "\n";
03870      iout << iINFO << "MODIFIED 1-4 VDW PARAMETERS WILL BE USED\n" << endi;
03871    } else {
03872      iout << iWARN << "MODIFIED 1-4 VDW PARAMETERS WILL BE IGNORED\n" << endi;
03873    }
03874 
03875 #ifdef SPEC_DISABLED_VERSION
03876    if (dcdFrequency > 0) {
03877      dcdFrequency = 0;
03878      iout << iWARN << "DCD TRAJECTORY OUTPUT IS DISABLED IN SPEC RELEASE\n";
03879    }
03880 #endif
03881 
03882    if (dcdFrequency > 0)
03883    {
03884      iout << iINFO << "DCD FILENAME           " 
03885         << dcdFilename << "\n";
03886      iout << iINFO << "DCD FREQUENCY          " 
03887         << dcdFrequency << "\n";
03888      iout << iINFO << "DCD FIRST STEP         " 
03889         << ( ((firstTimestep + dcdFrequency)/dcdFrequency)*dcdFrequency ) << "\n";
03890      if ( dcdUnitCell ) {
03891        iout << iINFO << "DCD FILE WILL CONTAIN UNIT CELL DATA\n";
03892      }
03893    }
03894    else
03895    {
03896      iout << iINFO << "NO DCD TRAJECTORY OUTPUT\n";
03897    }
03898    iout << endi;
03899    
03900    if (xstFrequency > 0)
03901    {
03902      iout << iINFO << "XST FILENAME           " 
03903         << xstFilename << "\n";
03904      iout << iINFO << "XST FREQUENCY          " 
03905         << xstFrequency << "\n";
03906    }
03907    else
03908    {
03909      iout << iINFO << "NO EXTENDED SYSTEM TRAJECTORY OUTPUT\n";
03910    }
03911    iout << endi;
03912    
03913    if (velDcdFrequency > 0)
03914    {
03915      iout << iINFO << "VELOCITY DCD FILENAME    " 
03916         << velDcdFilename << "\n";
03917      iout << iINFO << "VELOCITY DCD FREQUENCY   " 
03918         << velDcdFrequency << "\n";
03919      iout << iINFO << "VELOCITY DCD FIRST STEP  " 
03920         << ( ((firstTimestep + velDcdFrequency)/velDcdFrequency)*velDcdFrequency ) << "\n";
03921    }
03922    else
03923    {
03924      iout << iINFO << "NO VELOCITY DCD OUTPUT\n";
03925    }
03926    iout << endi;
03927    
03928    if (forceDcdFrequency > 0)
03929    {
03930      iout << iINFO << "FORCE DCD FILENAME     " 
03931         << forceDcdFilename << "\n";
03932      iout << iINFO << "FORCE DCD FREQUENCY    " 
03933         << forceDcdFrequency << "\n";
03934      iout << iINFO << "FORCE DCD FIRST STEP   " 
03935         << ( ((firstTimestep + forceDcdFrequency)/forceDcdFrequency)*forceDcdFrequency ) << "\n";
03936    }
03937    else
03938    {
03939      iout << iINFO << "NO FORCE DCD OUTPUT\n";
03940    }
03941    iout << endi;
03942    
03943    iout << iINFO << "OUTPUT FILENAME        " 
03944       << outputFilename << "\n" << endi;
03945    if (binaryOutput)
03946    {
03947      iout << iINFO << "BINARY OUTPUT FILES WILL BE USED\n" << endi;
03948    }
03949 #ifdef MEM_OPT_VERSION
03950     if(!binaryOutput){
03951         iout << iWARN <<"SINCE MEMORY OPTIMIZED VERSION IS USED, OUTPUT IN TEXT FORMAT IS DISABLED!\n" << endi;
03952         binaryOutput = TRUE;
03953     }
03954 #endif
03955 
03956    if (! restartFrequency)
03957    {
03958      iout << iINFO << "NO RESTART FILE\n";
03959    }
03960    else
03961    {
03962      iout << iINFO << "RESTART FILENAME       "
03963         << restartFilename << "\n";
03964      iout << iINFO << "RESTART FREQUENCY      " 
03965         << restartFrequency << "\n";
03966   if (restartSave) {
03967     iout << iINFO << "RESTART FILES WILL NOT BE OVERWRITTEN\n";
03968   }
03969 
03970   if (binaryRestart)
03971   {
03972     iout << iINFO << "BINARY RESTART FILES WILL BE USED\n";
03973   }
03974    }
03975    iout << endi;
03976    
03977    if (switchingActive)
03978    {
03979       iout << iINFO << "SWITCHING ACTIVE\n";
03980       if ( vdwForceSwitching ) {
03981         iout << iINFO << "VDW FORCE SWITCHING ACTIVE\n";
03982       }
03983       if ( martiniSwitching ) { 
03984         iout << iINFO << "MARTINI RESIDUE-BASED COARSE-GRAIN SWITCHING ACTIVE\n";
03985       }
03986       iout << iINFO << "SWITCHING ON           "
03987                << switchingDist << "\n";
03988       iout << iINFO << "SWITCHING OFF          "
03989                << cutoff << "\n";
03990    }
03991    else
03992    {
03993       iout << iINFO << "CUTOFF                 " 
03994          << cutoff << "\n";
03995    }
03996 
03997    iout << iINFO << "PAIRLIST DISTANCE      " << pairlistDist << "\n";
03998    iout << iINFO << "PAIRLIST SHRINK RATE   " << pairlistShrink << "\n";
03999    iout << iINFO << "PAIRLIST GROW RATE     " << pairlistGrow << "\n";
04000    iout << iINFO << "PAIRLIST TRIGGER       " << pairlistTrigger << "\n";
04001    iout << iINFO << "PAIRLISTS PER CYCLE    " << pairlistsPerCycle << "\n";
04002    if ( outputPairlists )
04003      iout << iINFO << "PAIRLIST OUTPUT STEPS  " << outputPairlists << "\n";
04004    iout << endi;
04005 
04006    if ( pairlistMinProcs > 1 )
04007      iout << iINFO << "REQUIRING " << pairlistMinProcs << " PROCESSORS FOR PAIRLISTS\n";
04008    usePairlists = ( CkNumPes() >= pairlistMinProcs );
04009 
04010 #ifdef OPENATOM_VERSION
04011 if ( openatomOn ) 
04012 {
04013   iout << iINFO << "OPENATOM QM/MM CAR-PARINELLO ACTIVE\n";
04014   iout << iINFO << "OPENATOM CONFIG FILE:  " << openatomConfig << "\n";
04015   iout << iINFO << "OPENATOM STRUCT FILE:  " << openatomStruct << "\n";
04016   iout << iINFO << "OPENATOM PDB FILE:     " << openatomPDB << "\n";
04017 }
04018 #endif // OPENATOM_VERSION
04019 
04020    // FB - FEP and TI are now dependent on pairlists - disallow usePairlists=0
04021    if ( (alchOn) && (!usePairlists)) {
04022      NAMD_die("Sorry, Alchemical simulations require pairlists to be enabled\n");
04023    }
04024 #ifdef NAMD_CUDA
04025    if ( ! usePairlists ) {
04026      usePairlists = 1;
04027      iout << iINFO << "CUDA ACCELERATION REQUIRES PAIRLISTS\n";
04028    }
04029 #endif
04030      
04031    iout << iINFO << "PAIRLISTS " << ( usePairlists ? "ENABLED" : "DISABLED" )
04032                                                         << "\n" << endi;
04033 
04034    iout << iINFO << "MARGIN                 " << margin << "\n";
04035 
04036    if ( splitPatch == SPLIT_PATCH_HYDROGEN ) {
04037       iout << iINFO << "HYDROGEN GROUP CUTOFF  " << hgroupCutoff << "\n";
04038    }
04039    
04040    iout << iINFO << "PATCH DIMENSION        "
04041             << patchDimension << "\n";
04042 
04043    iout << endi;
04044 
04045    if (outputEnergies != 1)
04046    {
04047       iout << iINFO << "ENERGY OUTPUT STEPS    "
04048          << outputEnergies << "\n";
04049       iout << endi;
04050    }
04051 
04052    if (mergeCrossterms) {
04053       iout << iINFO << "CROSSTERM ENERGY INCLUDED IN DIHEDRAL\n" << endi;
04054    }
04055    
04056    if (outputMomenta != 0)
04057    {
04058       iout << iINFO << "MOMENTUM OUTPUT STEPS  "
04059          << outputMomenta << "\n";
04060       iout << endi;
04061    }
04062    
04063    if (outputTiming != 0)
04064    {
04065       iout << iINFO << "TIMING OUTPUT STEPS    "
04066          << outputTiming << "\n";
04067       iout << endi;
04068    }
04069    
04070    if (outputCudaTiming != 0)
04071    {
04072       iout << iINFO << "CUDA TIMING OUTPUT STEPS    "
04073          << outputCudaTiming << "\n";
04074       iout << endi;
04075    }
04076    
04077    if (outputPressure != 0)
04078    {
04079       iout << iINFO << "PRESSURE OUTPUT STEPS  "
04080          << outputPressure << "\n";
04081       iout << endi;
04082    }
04083    
04084    if (fixedAtomsOn)
04085    {
04086       iout << iINFO << "FIXED ATOMS ACTIVE\n";
04087 #ifdef NAMD_CUDA
04088       iout << iWARN << "FOR CUDA VERSION PRESSURE WILL BE INCORRECT WITHOUT FIXEDATOMSFORCES\n";
04089       if (berendsenPressureOn || langevinPistonOn) {
04090         fixedAtomsForces = TRUE;
04091         iout << iWARN << "ENABLING FIXEDATOMSFORCES DUE TO PRESSURE CONTROL\n";
04092       }
04093 #endif
04094       if ( fixedAtomsForces )
04095         iout << iINFO << "FORCES BETWEEN FIXED ATOMS ARE CALCULATED\n";
04096       iout << endi;
04097    }
04098 
04099    if (constraintsOn)
04100    {
04101       iout << iINFO << "HARMONIC CONSTRAINTS ACTIVE\n";
04102 
04103       iout << iINFO << "HARMONIC CONS EXP      "
04104          << constraintExp << "\n";
04105 
04106       if (constraintScaling != 1.0) {
04107         iout << iINFO << "HARMONIC CONS SCALING  "
04108          << constraintScaling << "\n";
04109       }
04110 
04111       //****** BEGIN selective restraints (X,Y,Z) changes 
04112 
04113       if (selectConstraintsOn) {
04114         iout << iINFO << "SELECTED CARTESIAN COMPONENTS OF HARMONIC RESTRAINTS ACTIVE\n";
04115 
04116         if (constrXOn)
04117         iout << iINFO << "RESTRAINING X-COMPONENTS OF CARTESIAN COORDINATES!\n";
04118 
04119         if (constrYOn)
04120         iout << iINFO << "RESTRAINING Y-COMPONENTS OF CARTESIAN COORDINATES!\n";
04121 
04122         if (constrZOn)
04123         iout << iINFO << "RESTRAINING Z-COMPONENTS OF CARTESIAN COORDINATES!\n";
04124       }
04125       //****** END selective restraints (X,Y,Z) changes 
04126 
04127       if (sphericalConstraintsOn) {
04128         iout << iINFO << "SPHERICAL HARMONIC CONSTRAINTS ACTIVE\n";
04129         iout << iINFO << "RESTRAINING DISTANCE TO " << sphericalConstrCenter <<"\n";
04130       }
04131       iout << endi;
04132 
04133       //****** BEGIN moving constraints changes 
04134 
04135       if (movingConstraintsOn) {
04136         iout << iINFO << "MOVING HARMONIC CONSTRAINTS ACTIVE\n";
04137 
04138         iout << iINFO << "MOVING CONSTRAINT VELOCITY    "
04139              << movingConsVel << " ANGSTROM/TIMESTEP\n";
04140         
04141         iout << iINFO << "ALL CONSTRAINED ATOMS WILL MOVE\n";
04142       }
04143       //****** END moving constraints changes 
04144       iout << endi;
04145 
04146       //****** BEGIN rotating constraints changes 
04147 
04148       if (rotConstraintsOn) {
04149         iout << iINFO << "ROTATING HARMONIC CONSTRAINTS ACTIVE\n";
04150 
04151         iout << iINFO << "AXIS OF ROTATION    "
04152              << rotConsAxis << "\n";
04153         
04154         iout << iINFO << "PIVOT OF ROTATION   "
04155              << rotConsPivot << "\n";
04156 
04157         iout << iINFO << "ROTATING CONSTRAINT VELOCITY    "
04158              << rotConsVel << " DEGREES/TIMESTEP\n";
04159       }
04160       iout << endi;
04161       //****** END rotating constraints changes 
04162    }
04163 
04164    // moving drag
04165    if (movDragOn) {
04166      iout << iINFO << "MOVING DRAG ACTIVE.\n";
04167      
04168      iout << iINFO << "MOVING DRAG MAIN PDB FILE "
04169           << movDragFile << "\n";
04170      
04171      iout << iINFO << "MOVING DRAG GLOBAL VELOCITY (A/step) "
04172           << movDragGlobVel << "\n";
04173      
04174      iout << iINFO << "MOVING DRAG LINEAR VELOCITY FILE " 
04175           << movDragVelFile << "\n";
04176      
04177      iout << endi;
04178    }
04179    
04180    // rotating drag
04181    if (rotDragOn) {
04182      iout << iINFO << "ROTATING DRAG ACTIVE.\n";
04183      
04184      iout << iINFO << "ROTATING DRAG MAIN PDB FILE "
04185           << rotDragFile << "\n";
04186      
04187      iout << iINFO << "ROTATING DRAG AXIS FILE " 
04188           << rotDragAxisFile << "\n";
04189      
04190      iout << iINFO << "ROTATING DRAG PIVOT POINT FILE " 
04191           << rotDragPivotFile << "\n";
04192      
04193      iout << iINFO << "ROTATING DRAG GLOBAL ANGULAR VELOCITY (deg/step) "
04194           << rotDragGlobVel << "\n";
04195      
04196      iout << iINFO << "ROTATING DRAG ANGULAR VELOCITY FILE " 
04197           << rotDragVelFile << "\n";
04198      
04199      iout << endi;
04200    }
04201    
04202 
04203    // "constant" torque
04204    if (consTorqueOn) {
04205      iout << iINFO << "\"CONSTANT\" TORQUE ACTIVE.\n";
04206      
04207      iout << iINFO << "\"CONSTANT\" TORQUE MAIN PDB FILE "
04208           << consTorqueFile << "\n";
04209      
04210      iout << iINFO << "\"CONSTANT\" TORQUE AXIS FILE " 
04211           << consTorqueAxisFile << "\n";
04212      
04213      iout << iINFO << "\"CONSTANT\" TORQUE PIVOT POINT FILE " 
04214           << consTorquePivotFile << "\n";
04215      
04216      iout << iINFO << "\"CONSTANT\" TORQUE GLOBAL VALUE (Kcal/(mol*A^2)) "
04217           << consTorqueGlobVal << "\n";
04218      
04219      iout << iINFO << "\"CONSTANT\" TORQUE DACTORS FILE " 
04220           << consTorqueValFile << "\n";
04221      
04222      iout << endi;
04223    }
04224 
04225    if (mgridforceOn) {
04226      iout << iINFO << "GRID FORCE ACTIVE\n";
04227      iout << iINFO << " Please include this reference in published work using\n";
04228      iout << iINFO << " the Gridforce module of NAMD: David Wells, Volha Abramkina,\n";
04229      iout << iINFO << " and Aleksei Aksimentiev, J. Chem. Phys. 127:125101-10 (2007).\n";
04230      print_mgrid_params();
04231    }
04232    
04233    //****** BEGIN SMD constraints changes 
04234    
04235    if (SMDOn) {
04236      iout << iINFO << "SMD ACTIVE\n";
04237      
04238      iout << iINFO << "SMD VELOCITY    "
04239           << SMDVel << " ANGSTROM/TIMESTEP\n";
04240         
04241      iout << iINFO << "SMD DIRECTION   "
04242           << SMDDir << "\n";
04243  
04244      iout << iINFO << "SMD K   " 
04245           << SMDk << "\n";
04246 
04247      iout << iINFO << "SMD K2  " 
04248           << SMDk2 << "\n";
04249 
04250      iout << iINFO << "SMD OUTPUT FREQUENCY   "
04251           << SMDOutputFreq << " TIMESTEPS\n";
04252     
04253      iout << iINFO << "SMD FILE " << SMDFile << "\n"; 
04254 
04255      iout << endi;
04256    }
04257    
04258    //****** END SMD constraints changes 
04259 
04260    if (TMDOn) {
04261      iout << iINFO << "TMD ACTIVE BETWEEN STEPS " << TMDFirstStep 
04262           << " and " << TMDLastStep << "\n";
04263      iout << iINFO << "TMD K  " << TMDk << "\n";
04264      iout << iINFO << "TMD FILE  " << TMDFile << "\n";
04265      iout << iINFO << "TMD OUTPUT FREQUENCY  " << TMDOutputFreq << "\n";
04266      if (TMDInitialRMSD) {
04267        iout << iINFO << "TMD TARGET RMSD AT FIRST STEP  " << TMDInitialRMSD << "\n";
04268      } else {
04269        iout << iINFO << "TMD TARGET RMSD AT FIRST STEP COMPUTED FROM INITIAL COORDINATES\n";
04270      }
04271      iout << iINFO << "TMD TARGET RMSD AT FINAL STEP  " << TMDFinalRMSD << "\n";
04272      iout << endi;
04273    }
04274 
04275    if (symmetryOn) {
04276      if (symmetryLastStep == -1){
04277        iout << iINFO << "SYMMETRY RESTRAINTS ACTIVE BETWEEN STEPS " << symmetryFirstStep << " and " << "INFINITY" << "\n";
04278      }
04279      else{
04280        iout << iINFO << "SYMMETRY RESTRAINTS ACTIVE BETWEEN STEPS " << symmetryFirstStep << " and " << symmetryLastStep << "\n";
04281      }
04282     // iout << iINFO << "SYMMETRY FILE " << symmetryFile << "\n";
04283 
04284      current = config->find("symmetryFile");
04285      for ( ; current; current = current->next ) {
04286        iout << iINFO << "SYMMETRY FILE  " << current->data << "\n";
04287      }
04288 
04289      current = config->find("symmetryMatrixFile");
04290      for ( ; current; current = current->next ) {
04291       iout << iINFO << "SYMMETRY MATRIX FILE " << current->data << "\n";
04292      }
04293      iout << iINFO << "SYMMETRY FORCE CONSTANT " << symmetryk << "\n";
04294      if (symmetryScaleForces){
04295       iout << iINFO << "SYMMETRY SCALE FORCES ON\n";
04296      }
04297      iout << iINFO << "SYMMETRY FIRST FULL STEP " << symmetryFirstFullStep << "\n";
04298      if (symmetryLastFullStep == -1){
04299       iout << iINFO << "SYMMETRY LAST FULL STEP " << "INFINITY" << "\n"; 
04300       //iout << iINFO << "FULL SYMMETRY FORCE BETWEEN STEPS " << symmetryFirstFullStep << " and " << "INFINITY" << "\n";
04301      }
04302      else {
04303       iout << iINFO << "SYMMETRY LAST FULL STEP " << symmetryLastFullStep << "\n"; 
04304      // iout << iINFO << "FULL SYMMETRY FORCE BETWEEN STEPS " << symmetryFirstFullStep << " and " << symmetryLastFullStep << "\n";
04305      }
04306 
04307      iout << endi;
04308    }
04309 //Modifications for alchemical fep
04310 //  Alchemical FEP status
04311 
04312 //   current = config->find("alchOutFile");
04313    if (alchFepOn)
04314    {
04315      iout << iINFO << "ALCHEMICAL FEP ON\n";
04316      iout << iINFO << "FEP CURRENT LAMBDA VALUE     "
04317           << alchLambda << "\n";
04318      iout << iINFO << "FEP COMPARISON LAMBDA VALUE  "
04319           << alchLambda2 << "\n";
04320      iout << iINFO << "FEP VDW SHIFTING COEFFICIENT "
04321           << alchVdwShiftCoeff << "\n";
04322      iout << iINFO << "FEP ELEC. ACTIVE FOR ANNIHILATED "
04323           << "PARTICLES BETWEEN LAMBDA = 0 AND LAMBDA = "
04324           << (1 - alchElecLambdaStart) << "\n";
04325      iout << iINFO << "FEP ELEC. ACTIVE FOR EXNIHILATED "
04326           << "PARTICLES BETWEEN LAMBDA = "
04327           << alchElecLambdaStart << " AND LAMBDA = 1\n";
04328      iout << iINFO << "FEP VDW ACTIVE FOR ANNIHILATED "
04329           << "PARTICLES BETWEEN LAMBDA = "
04330           << (1 - alchVdwLambdaEnd) << " AND LAMBDA = 1\n";
04331      iout << iINFO << "FEP VDW ACTIVE FOR EXNIHILATED "
04332           << "PARTICLES BETWEEN LAMBDA = 0 AND LAMBDA = "
04333           << alchVdwLambdaEnd << "\n";
04334      if (alchFepWCADispOn)
04335      {
04336        iout << iINFO << "FEP WEEKS-CHANDLER-ANDERSEN DECOMPOSITION (DISPERSION) ON\n";
04337      }
04338      if (alchFepWCARepuOn)
04339      {
04340        iout << iINFO << "FEP WEEKS-CHANDLER-ANDERSEN DECOMPOSITION (REPULSION) ON\n";
04341        iout << iINFO << "FEP WEEKS-CHANDLER-ANDERSEN RCUT1 = " 
04342             << alchFepWCArcut1 << " , RCUT2 = " 
04343             << alchFepWCArcut2  << " AND RCUT3 = " << alchFepWCArcut3 << "\n";
04344      }
04345    }
04346 //fepe
04347 
04348    if (alchThermIntOn)
04349    {
04350      iout << iINFO << "THERMODYNAMIC INTEGRATION (TI) ON\n";
04351      iout << iINFO << "TI LAMBDA VALUE     "
04352           << alchLambda << "\n";
04353      iout << iINFO << "TI VDW SHIFTING COEFFICIENT "
04354           << alchVdwShiftCoeff << "\n";
04355      iout << iINFO << "TI ELEC. ACTIVE FOR ANNIHILATED "
04356           << "PARTICLES BETWEEN LAMBDA = 0 AND LAMBDA = "
04357           << (1 - alchElecLambdaStart) << "\n";
04358      iout << iINFO << "TI ELEC. ACTIVE FOR EXNIHILATED "
04359           << "PARTICLES BETWEEN LAMBDA = "
04360           << alchElecLambdaStart << " AND LAMBDA = 1\n";
04361      iout << iINFO << "TI VDW ACTIVE FOR ANNIHILATED "
04362           << "PARTICLES BETWEEN LAMBDA = "
04363           << (1 - alchVdwLambdaEnd) << " AND LAMBDA = 1\n";
04364      iout << iINFO << "TI VDW ACTIVE FOR EXNIHILATED "
04365           << "PARTICLES BETWEEN LAMBDA = 0 AND LAMBDA = "
04366           << alchVdwLambdaEnd << "\n";
04367    }
04368 
04369 
04370    if ( lesOn ) {
04371      iout << iINFO << "LOCALLY ENHANCED SAMPLING ACTIVE\n";
04372      iout << iINFO << "LOCAL ENHANCEMENT FACTOR IS "
04373           << lesFactor << "\n";
04374      if ( lesReduceTemp ) iout << iINFO
04375        << "SCALING ENHANCED ATOM TEMPERATURE BY 1/" << lesFactor << "\n";
04376      if ( lesReduceMass ) iout << iINFO
04377        << "SCALING ENHANCED ATOM MASS BY 1/" << lesFactor << "\n";
04378    }
04379    
04380    if ( pairInteractionOn ) {
04381      iout << iINFO << "PAIR INTERACTION CALCULATIONS ACTIVE\n";
04382      iout << iINFO << "USING FLAG " << pairInteractionGroup1 
04383           << " FOR GROUP 1\n";
04384      if (pairInteractionSelf) {
04385        iout << iINFO << "COMPUTING ONLY SELF INTERACTIONS FOR GROUP 1 ATOMS\n";
04386      } else {
04387        iout << iINFO << "USING FLAG " << pairInteractionGroup2 
04388             << " FOR GROUP 2\n";
04389      }
04390    }
04391 
04392    if (consForceOn) {
04393      iout << iINFO << "CONSTANT FORCE ACTIVE\n";
04394      if ( consForceScaling != 1.0 ) {
04395        iout << iINFO << "CONSTANT FORCE SCALING   "
04396                                 << consForceScaling << "\n" << endi;
04397      }
04398    }
04399 
04400 
04401    // external command forces
04402 
04403    if (extForcesOn) {
04404      iout << iINFO << "EXTERNAL COMMAND FORCES ACTIVE\n";
04405      iout << iINFO << "EXT FORCES COMMAND: " << extForcesCommand << "\n";
04406      iout << iINFO << "EXT COORD FILENAME: " << extCoordFilename << "\n";
04407      iout << iINFO << "EXT FORCE FILENAME: " << extForceFilename << "\n";
04408      iout << endi;
04409    }
04410 
04411    // gbis gbobc implicit solvent parameters
04412 
04413   if (GBISserOn) {
04414       GBISOn = 0;//turning gbis-ser on turns gbis-parallel off
04415      iout << iINFO<< "GBIS GENERALIZED BORN IMPLICIT SOLVENT ACTIVE (SERIAL)\n";
04416   }
04417   if (GBISOn) {
04418      iout << iINFO << "GBIS GENERALIZED BORN IMPLICIT SOLVENT ACTIVE\n";
04419   }
04420   if (GBISOn || GBISserOn) {
04421      iout << iINFO << "GBIS SOLVENT DIELECTRIC: " << solvent_dielectric<< "\n";
04422      iout << iINFO << "GBIS PROTEIN DIELECTRIC: " << dielectric<< "\n";
04423      iout <<iINFO<<"GBIS COULOMB RADIUS OFFSET: "<< coulomb_radius_offset<<" Ang\n";
04424      iout << iINFO << "GBIS ION CONCENTRATION: " << ion_concentration << " M\n";
04425      iout << iINFO << "GBIS DEBYE SCREENING LENGTH: " << 1.0/kappa << " Ang\n";
04426      iout << iINFO << "GBIS DELTA: " << gbis_delta << "\n";
04427      iout << iINFO << "GBIS BETA: " << gbis_beta << "\n";
04428      iout << iINFO << "GBIS GAMMA: " << gbis_gamma << "\n";
04429      iout << iINFO << "GBIS BORN RADIUS CUTOFF: " << alpha_cutoff << " Ang\n";
04430      iout << iINFO << "GBIS MAX BORN RADIUS: " << alpha_max << " Ang\n";
04431      iout << endi;
04432   }
04433 
04434   if (LCPOOn) {
04435     iout << iINFO << "SASA SURFACE TENSION: " << surface_tension<< " kcal/mol/Ang^2\n";
04436   }
04437 
04438    tclBCScript = 0;
04439    if (tclBCOn) {
04440      iout << iINFO << "TCL BOUNDARY FORCES ACTIVE\n";
04441      current = config->find("tclBCScript");
04442      if ( current ) {
04443        tclBCScript = current->data;
04444        iout << iINFO << "TCL BOUNDARY FORCES SCRIPT   " << current->data << "\n";
04445      }
04446        iout << iINFO << "TCL BOUNDARY FORCES ARGS     " << tclBCArgs << "\n";
04447      iout << endi;
04448    }
04449    
04450    // Global forces configuration
04451 
04452    globalForcesOn = ( tclForcesOn || freeEnergyOn || miscForcesOn ||
04453                       (IMDon && ! IMDignore) || SMDOn || TMDOn || colvarsOn || symmetryOn );
04454 
04455 
04456    if (tclForcesOn)
04457    {
04458      iout << iINFO << "TCL GLOBAL FORCES ACTIVE\n";
04459 
04460      current = config->find("tclForcesScript");
04461 
04462      for ( ; current; current = current->next ) {
04463 
04464      if ( strstr(current->data,"\n") ) {
04465        iout << iINFO << "TCL GLOBAL FORCES SCRIPT INLINED IN CONFIG FILE\n";
04466        continue;
04467      }
04468 
04469      iout << iINFO << "TCL GLOBAL FORCES SCRIPT   " << current->data << "\n";
04470 
04471      }
04472      iout << endi;
04473    }
04474 
04475    if (miscForcesOn)
04476    {
04477      iout << iINFO << "MISC FORCES ACTIVE\n";
04478 
04479      current = config->find("miscForcesScript");
04480 
04481      for ( ; current; current = current->next ) {
04482 
04483      if ( strstr(current->data,"\n") ) {
04484        iout << iINFO << "MISC FORCES SCRIPT INLINED IN CONFIG FILE\n";
04485        continue;
04486      }
04487 
04488      iout << iINFO << "MISC FORCES SCRIPT   " << current->data << "\n";
04489 
04490      }
04491      iout << endi;
04492    }
04493 
04494    if (freeEnergyOn)
04495    {
04496      iout << iINFO << "FREE ENERGY PERTURBATION ACTIVE\n";
04497 
04498      current = config->find("freeEnergyConfig");
04499 
04500      for ( ; current; current = current->next ) {
04501 
04502      if ( strstr(current->data,"\n") ) {
04503        iout << iINFO << "FREE ENERGY PERTURBATION SCRIPT INLINED IN CONFIG FILE\n";
04504        continue;
04505      }
04506 
04507      iout << iINFO << "FREE ENERGY PERTURBATION SCRIPT   " << current->data << "\n";
04508 
04509      }
04510      iout << endi;
04511    }
04512 
04513    if (colvarsOn)
04514    {
04515      iout << iINFO << "COLLECTIVE VARIABLES CALCULATION REQUESTED\n";
04516 
04517      current = config->find ("colvarsConfig");
04518      for ( ; current; current = current->next ) {
04519        if ( strstr(current->data,"\n") ) {
04520          iout << iINFO << "COLLECTIVE VARIABLES CONFIGURATION INLINED IN CONFIG FILE\n";
04521          continue;
04522        }
04523        iout << iINFO << "COLLECTIVE VARIABLES CONFIGURATION   " << current->data << "\n";
04524      }
04525 
04526      current = config->find ("colvarsInput");
04527      for ( ; current; current = current->next ) {
04528        if ( strstr(current->data,"\n") ) {
04529          iout << iINFO << "COLLECTIVE VARIABLES RESTART INFORMATION INLINED IN CONFIG FILE\n";
04530          continue;
04531        }
04532        iout << iINFO << "COLLECTIVE VARIABLES RESTART INFORMATION   " << current->data << "\n";
04533      }
04534 
04535      iout << endi;
04536    }
04537 
04538    if (IMDon)
04539    {
04540      iout << iINFO << "INTERACTIVE MD ACTIVE\n";
04541      iout << iINFO << "INTERACTIVE MD PORT    " << IMDport << "\n";
04542      iout << iINFO << "INTERACTIVE MD FREQ    " << IMDfreq << "\n";
04543      if (IMDignore) {
04544         iout << iINFO << "INTERACTIVE MD WILL NOT INFLUENCE SIMULATION\n";
04545      } else {
04546        if (IMDwait) iout << iINFO << "WILL AWAIT INTERACTIVE MD CONNECTION\n";
04547      }
04548      iout << endi;
04549    }
04550 
04551    if (globalOn && !dihedralOn)
04552    {
04553       iout << iINFO << "GLOBAL INTEGRATION TEST MODE ACTIVE\n";
04554    }
04555 
04556 
04557    if (dihedralOn)
04558    {
04559       iout << iINFO << "DIHEDRAL ANGLE DYNAMICS ACTIVE\n";
04560       if (!COLDOn)
04561       {
04562          iout << iINFO << "*** DIHEDRAL ANGLE DYNAMICS IS HIGHLY EXPERIMENTAL ***\n";
04563          iout << iINFO << "PLEASE CONSIDER USING THE COLD OPTION AS WELL\n";
04564       }
04565    }
04566 
04567    // This function is so long that it exceeds the emacs brace
04568    // matching default max length of 25600 (a bit before this comment). We
04569    // should take this is a not too subtle hint about scale of this
04570    // violation of good coding practices.  Especially considering the
04571    // fact that this function still has about a thousand lines to go
04572    // before its done, and is doomed to grow with new features. 
04573 
04574    if (COLDOn)
04575    {
04576       iout << iINFO << "COLD (CONSTRAINED OVERDAMPED LANGEVIN DYNAMICS) ACTIVE\n";
04577 
04578       iout << iINFO << "COLD TARGET TEMP       "
04579          << COLDTemp << "\n";
04580 
04581       iout << iINFO << "COLD COLLISION RATE    "
04582          << COLDRate << "\n";
04583    }
04584 
04585    if (cylindricalBCOn)
04586    {
04587     iout << iINFO << "CYLINDRICAL BOUNDARY CONDITIONS ACTIVE\n";
04588     iout << iINFO << "AXIS                     " << cylindricalBCAxis << "\n";
04589     iout << iINFO << "RADIUS #1                " << cylindricalBCr1 << "\n";
04590     iout << iINFO << "FORCE CONSTANT #1        " << cylindricalBCk1 << "\n";
04591     iout << iINFO << "EXPONENT #1              " << cylindricalBCexp1 << "\n";
04592     iout << iINFO << "LENGTH #1                " << cylindricalBCl1 << "\n";
04593     if (cylindricalBCr2 > 0.0)
04594     {
04595      iout << iINFO << "RADIUS #2               " << cylindricalBCr2 << "\n";
04596      iout << iINFO << "FORCE CONSTANT #2       " << cylindricalBCk2 << "\n";
04597      iout << iINFO << "EXPONENT #2             " << cylindricalBCexp2 << "\n";
04598      iout << iINFO << "LENGTH #2               " << cylindricalBCl2 << "\n";
04599     }
04600     iout << iINFO << "CYLINDER BOUNDARY CENTER(" << cylindricalCenter.x << ", "
04601              << cylindricalCenter.y << ", " << cylindricalCenter.z << ")\n";
04602     iout << endi;
04603   }
04604 
04605    if (sphericalBCOn)
04606    {
04607       iout << iINFO << "SPHERICAL BOUNDARY CONDITIONS ACTIVE\n";
04608 
04609       iout << iINFO << "RADIUS #1              "
04610          << sphericalBCr1 << "\n";
04611       iout << iINFO << "FORCE CONSTANT #1      "
04612          << sphericalBCk1 << "\n";
04613       iout << iINFO << "EXPONENT #1            "
04614          << sphericalBCexp1 << "\n";
04615 
04616       if (sphericalBCr2 > 0)
04617       {
04618         iout << iINFO << "RADIUS #2              "
04619               << sphericalBCr2 << "\n";
04620         iout << iINFO << "FORCE CONSTANT #2      "
04621             << sphericalBCk2 << "\n";
04622         iout << iINFO << "EXPONENT #2            "
04623         << sphericalBCexp2 << "\n";
04624       }
04625 
04626       iout << iINFO << "SPHERE BOUNDARY CENTER(" << sphericalCenter.x << ", "
04627                << sphericalCenter.y << ", " << sphericalCenter.z << ")\n";
04628       iout << endi;
04629    }
04630 
04631    if (eFieldOn)
04632    {
04633       iout << iINFO << "ELECTRIC FIELD ACTIVE\n";
04634       
04635       iout << iINFO << "E-FIELD VECTOR         ("
04636          << eField.x << ", " << eField.y
04637          << ", " << eField.z << ")\n";
04638       iout << iINFO << "E-FIELD FREQUENCY IS (1/ps) " << eFieldFreq << "\n";
04639       iout << iINFO << "E-FIELD PHASE IS     (deg)  " << eFieldPhase << "\n";
04640 
04641       iout << endi;
04642    }
04643 
04644       if (stirOn)
04645    {
04646       iout << iINFO << "STIRRING TORQUES ACTIVE\n";
04647       
04648       iout << iINFO << "STIR STARTING THETA   (deg)  "<< stirStartingTheta << "\n";
04649       iout << iINFO << "STIR ANGULAR VELOCITY (deg/ts)   " << stirVel <<"\n";
04650       iout << iINFO << "STIR FORCE HARMONIC SPRING CONSTANT "<< stirK << "\n";
04651       iout << iINFO << "STIR AXIS OF ROTATION (DIRECTION)      ("
04652          << stirAxis.x << ", " << stirAxis.y
04653          << ", " << stirAxis.z << ")\n";
04654                   iout << iINFO << "STIR PIVOT POINT (COORDINATE)           ("
04655          << stirPivot.x << ", " << stirPivot.y
04656          << ", " << stirPivot.z << ")\n";
04657       current = config->find("stirFilename");
04658                   
04659       iout << iINFO << "STIR ATOMS AND ORIGINAL POSITIONS FROM FILE    " <<current ->data << '\n';
04660       current = config->find("stirredAtomsCol");
04661       iout << iINFO <<"STIR FILE COLUMN " << current ->data << '\n';
04662       iout << endi;
04663    }
04664 
04665    if (drudeOn)
04666    {
04667       iout << iINFO << "DRUDE MODEL DUAL THERMOSTAT IS ACTIVE\n";
04668       iout << iINFO << "DRUDE BOND TEMPERATURE " << drudeTemp << "\n";
04669       if (drudeDamping > 0.0) {
04670         iout << iINFO << "DRUDE DAMPING COEFFICIENT IS "
04671              << drudeDamping << " INVERSE PS\n";
04672       }
04673       if (drudeBondConst > 0.0) {
04674         iout << iINFO << "DRUDE QUARTIC RESTRAINT IS ACTIVE FOR DRUDE BONDS\n";
04675         iout << iINFO << "DRUDE MAXIMUM BOND LENGTH BEFORE RESTRAINT IS   "
04676              << drudeBondLen << "\n";
04677         iout << iINFO << "DRUDE BOND RESTRAINT CONSTANT IS                "
04678              << drudeBondConst << "\n";
04679       }
04680       if (drudeNbtholeCut > 0.0) {
04681         iout << iINFO << "DRUDE NBTHOLE IS ACTIVE\n";
04682         iout << iINFO << "DRUDE NBTHOLE RADIUS IS   "
04683              << drudeNbtholeCut << "\n";
04684       }
04685    }
04686 
04687    if (langevinOn)
04688    {
04689       iout << iINFO << "LANGEVIN DYNAMICS ACTIVE\n";
04690       iout << iINFO << "LANGEVIN TEMPERATURE   "
04691          << langevinTemp << "\n";
04692       if (! langevin_useBAOAB) iout << iINFO << "LANGEVIN USING BBK INTEGRATOR\n";
04693       else  iout << iINFO << "LANGEVIN USING BAOAB INTEGRATOR\n"; // [!!] Info file
04694       if (langevinDamping > 0.0) {
04695         iout << iINFO << "LANGEVIN DAMPING COEFFICIENT IS "
04696                 << langevinDamping << " INVERSE PS\n";
04697         if (langevinHydrogen)
04698                 iout << iINFO << "LANGEVIN DYNAMICS APPLIED TO HYDROGENS\n";
04699         else
04700                 iout << iINFO << "LANGEVIN DYNAMICS NOT APPLIED TO HYDROGENS\n";
04701       } else {
04702         iout << iINFO << "LANGEVIN DAMPING COEFFICIENTS DETERMINED FROM FILES\n";
04703         current = config->find("langevinFile");
04704         if ( current ) iout << iINFO << "LANGEVIN DAMPING FILE:  " <<
04705           current->data << "\n";
04706         else iout << iINFO << "LANGEVIN DAMPING FILE IS COORDINATE PDB\n";
04707         current = config->find("langevinCol");
04708         if ( current ) iout << iINFO << "LANGEVIN DAMPING COLUMN:  " <<
04709           current->data << "\n";
04710         else iout << iINFO << "LANGEVIN DAMPING COLUMN:  DEFAULT (4TH, O)\n";
04711       }
04712       iout << endi;
04713    }
04714    
04715    // BEGIN LA
04716    if (loweAndersenOn)
04717    {
04718       iout << iINFO << "LOWE-ANDERSEN DYNAMICS ACTIVE\n";
04719       iout << iINFO << "LOWE-ANDERSEN TEMPERATURE     "
04720          << loweAndersenTemp << " K\n";
04721       iout << iINFO << "LOWE-ANDERSEN RATE            "
04722          << loweAndersenRate << " INVERSE PS\n";
04723       iout << iINFO << "LOWE-ANDERSEN CUTOFF          "
04724          << loweAndersenCutoff << " ANGSTROMS\n";
04725       iout << endi;
04726    }
04727    // END LA
04728 
04729    if (tCoupleOn)
04730    {
04731       iout << iINFO << "TEMPERATURE COUPLING ACTIVE\n";
04732       iout << iINFO << "COUPLING TEMPERATURE   "
04733          << tCoupleTemp << "\n";
04734       iout << endi;
04735    }
04736 
04737    if (minimizeOn)
04738    {
04739       iout << iINFO << "OLD STYLE MINIMIZATION ACTIVE\n";
04740       iout << endi;
04741    }
04742 
04743    if (minimizeCGOn)
04744    {
04745       iout << iINFO << "CONJUGATE GRADIENT MINIMIZATION ACTIVE\n";
04746       iout << iINFO << "LINE MINIMIZATION GOAL = " << minLineGoal << "\n";
04747       iout << iINFO << "BABY STEP SIZE = " << minBabyStep << "\n";
04748       iout << iINFO << "TINY STEP SIZE = " << minTinyStep << "\n";
04749       iout << endi;
04750    }
04751 
04752    if (maximumMove)
04753    {
04754       iout << iINFO << "MAXIMUM MOVEMENT       "
04755          << maximumMove << "\n";
04756       iout << endi;
04757    }
04758 
04759    if (rescaleFreq > 0)
04760    {
04761      iout << iINFO << "VELOCITY RESCALE FREQ  "
04762         << rescaleFreq << "\n";
04763      iout << iINFO << "VELOCITY RESCALE TEMP  "
04764         << rescaleTemp << "\n";
04765      iout << endi;
04766    }
04767 
04768    if (reassignFreq > 0)
04769    {
04770      iout << iINFO << "VELOCITY REASSIGNMENT FREQ  "
04771         << reassignFreq << "\n";
04772      iout << iINFO << "VELOCITY REASSIGNMENT TEMP  "
04773         << reassignTemp << "\n";
04774      if ( reassignIncr != 0. )
04775        iout << iINFO << "VELOCITY REASSIGNMENT INCR  "
04776         << reassignIncr << "\n";
04777      if ( reassignHold != 0. )
04778        iout << iINFO << "VELOCITY REASSIGNMENT HOLD  "
04779         << reassignHold << "\n";
04780      iout << endi;
04781    }
04782 
04783    if (berendsenPressureOn && langevinPistonOn)
04784    {
04785       NAMD_die("Multiple pressure control algorithms selected!\n");
04786    }
04787 
04788    if (excludeFromPressure) {
04789      iout << iINFO << "EXCLUDE FROM PRESSURE ACTIVE\n";
04790    }
04791    if (useConstantArea && useConstantRatio) {
04792      NAMD_die("useConstantArea and useConstantRatio are mutually exclusive.\n");
04793    }
04794    if (useConstantRatio && !useFlexibleCell) {
04795      NAMD_die("useConstantRatio requires useFlexibleCell.\n");
04796    }
04797    if (useConstantArea && surfaceTensionTarget) {
04798      NAMD_die("surfaceTensionTarget and useConstantArea are mutually exclusive.\n");
04799    }
04800    if (useConstantArea && !useFlexibleCell) {
04801      NAMD_die("useConstantArea requires useFlexibleCell.\n");
04802    }
04803 
04804    if (berendsenPressureOn || langevinPistonOn) {
04805      if (rigidBonds != RIGID_NONE && useGroupPressure == FALSE) {
04806        useGroupPressure = TRUE;
04807        iout << iWARN << "Option useGroupPressure is being enabled "
04808             << "due to pressure control with rigidBonds.\n" << endi;
04809      }
04810    }
04811 
04812    if (berendsenPressureOn)
04813    {
04814      if ( ! opts.defined("BerendsenPressureFreq") ) {
04815         berendsenPressureFreq = nonbondedFrequency;
04816         if ( fullElectFrequency )
04817                 berendsenPressureFreq = fullElectFrequency;
04818      }
04819      if ( (berendsenPressureFreq % nonbondedFrequency) || ( fullElectFrequency
04820                 && (berendsenPressureFreq % fullElectFrequency) ) )
04821         NAMD_die("berendsenPressureFreq must be a multiple of both fullElectFrequency and nonbondedFrequency\n");
04822      iout << iINFO << "BERENDSEN PRESSURE COUPLING ACTIVE\n";
04823      iout << iINFO << "    TARGET PRESSURE IS "
04824         << berendsenPressureTarget << " BAR\n";
04825      iout << iINFO << "    COMPRESSIBILITY ESTIMATE IS "
04826         << berendsenPressureCompressibility << " BAR^(-1)\n";
04827      iout << iINFO << "    RELAXATION TIME IS "
04828         << berendsenPressureRelaxationTime << " FS\n";
04829      iout << iINFO << "    APPLIED EVERY "
04830         << berendsenPressureFreq << " STEPS\n";
04831      iout << iINFO << "    PRESSURE CONTROL IS "
04832         << (useGroupPressure?"GROUP":"ATOM") << "-BASED\n";
04833      iout << endi;
04834      berendsenPressureTarget /= PRESSUREFACTOR;
04835      berendsenPressureCompressibility *= PRESSUREFACTOR;
04836    }
04837 
04838    if (langevinPistonOn)
04839    {
04840      iout << iINFO << "LANGEVIN PISTON PRESSURE CONTROL ACTIVE\n";
04841      iout << iINFO << "       TARGET PRESSURE IS "
04842         << langevinPistonTarget << " BAR\n";
04843      iout << iINFO << "    OSCILLATION PERIOD IS "
04844         << langevinPistonPeriod << " FS\n";
04845      iout << iINFO << "            DECAY TIME IS "
04846         << langevinPistonDecay << " FS\n";
04847      iout << iINFO << "    PISTON TEMPERATURE IS "
04848         << langevinPistonTemp << " K\n";
04849      iout << iINFO << "      PRESSURE CONTROL IS "
04850         << (useGroupPressure?"GROUP":"ATOM") << "-BASED\n";
04851      iout << iINFO << "   INITIAL STRAIN RATE IS "
04852         << strainRate << "\n";
04853      iout << endi;
04854      langevinPistonTarget /= PRESSUREFACTOR;
04855    }
04856 
04857    if (berendsenPressureOn || langevinPistonOn) {
04858      iout << iINFO << "      CELL FLUCTUATION IS "
04859             << (useFlexibleCell?"AN":"") << "ISOTROPIC\n";
04860      if (useConstantRatio) 
04861        iout << iINFO << "    SHAPE OF CELL IS CONSTRAINED IN X-Y PLANE\n";
04862      if (useConstantArea) 
04863        iout << iINFO << "    CONSTANT AREA PRESSURE CONTROL ACTIVE\n";
04864    }
04865 
04866    if (surfaceTensionTarget != 0)
04867    {
04868      iout << iINFO << "SURFACE TENSION CONTROL ACTIVE\n";
04869      iout << iINFO << "      TARGET SURFACE TENSION IS "
04870           << surfaceTensionTarget << " DYN/CM\n";
04871      iout << endi;
04872      // multiply by 100 to convert from dyn/cm to bar-Angstroms, then divide
04873      // by PRESSURE factor to convert bar to NAMD internal pressure units. 
04874      surfaceTensionTarget *= 100.0 / PRESSUREFACTOR;
04875    }
04876 
04877    if (pressureProfileOn) {
04878      if ((berendsenPressureOn || langevinPistonOn) && !dcdUnitCell) {
04879 #if 1
04880        iout << iWARN << "Turning on dcdUnitCell so that trajectory files contain unit cell data.\n" << endi;
04881        dcdUnitCell = 1;
04882 #else
04883        NAMD_die("Sorry, pressure profile not implemented for constant pressure.");
04884 #endif
04885      }
04886      // if Ewald is on, only calculate Ewald
04887      if (pressureProfileEwaldOn)
04888        pressureProfileOn = 0;
04889 
04890      if (pressureProfileSlabs < 1) 
04891        NAMD_die("pressureProfileSlabs must be positive.");
04892      iout << iINFO << "PRESSURE PROFILE CALCULATIONS ACTIVE\n";
04893      iout << iINFO << "      NUMBER OF SLABS: " << pressureProfileSlabs << "\n";
04894      iout << iINFO << "      SLAB THICKNESS: " << cellBasisVector3.z / pressureProfileSlabs
04895                    << "\n";
04896      iout << iINFO << "      TIMESTEPS BETWEEN DATA OUTPUT: " 
04897                    << pressureProfileFreq << "\n";
04898      iout << iINFO << "      NUMBER OF ATOM TYPES: " << pressureProfileAtomTypes << "\n";
04899      iout << endi;
04900    } else {
04901      pressureProfileEwaldOn = 0;
04902      pressureProfileAtomTypes = 1;
04903    }
04904 
04905    if (accelMDOn) {
04906      iout << iINFO << "ACCELERATED MD ACTIVE\n";
04907 
04908      if ( accelMDdual) {
04909         accelMDdihe = FALSE;
04910         iout << iINFO << "APPLYING DUAL BOOST\n";
04911      }
04912      else if ( accelMDdihe ) {
04913         iout << iINFO << "BOOSTING DIHEDRAL POTENTIAL\n";
04914      } else {
04915         iout << iINFO << "BOOSTING TOTAL POTENTIAL\n";
04916      }
04917 
04918      iout << iINFO << "accelMDE: " << accelMDE << " KCAL/MOL, accelMDalpha: " << accelMDalpha << " KCAL/MOL\n";
04919      if (accelMDdual) {
04920         iout << iINFO << "accelMDTE: " << accelMDTE << " KCAL/MOL, accelMDTalpha: " << accelMDTalpha << " KCAL/MOL\n";
04921      }
04922      if ( accelMDLastStep > 0) {
04923         iout << iINFO << "accelMD WILL BE DONE FROM STEP " << accelMDFirstStep << " TO STEP " << accelMDLastStep << "\n";
04924      } else {
04925         iout << iINFO << "accelMD WILL BE DONE FROM STEP " << accelMDFirstStep << " TO THE END OF THE SIMULATION \n";
04926      }        
04927      iout << iINFO << "accelMD OUTPUT FREQUENCY " << accelMDOutFreq << "\n";
04928      iout << endi;
04929    }
04930 
04931    if (adaptTempOn) {
04932      iout << iINFO << "ADAPTIVE TEMPERING ACTIVE:\n";
04933      iout << iINFO << "      OUTPUT FREQUENCY: " << adaptTempOutFreq << "\n";
04934      iout << iINFO << "      TEMPERATURE UPDATE FREQUENCY: " << adaptTempFreq << "\n";
04935      if ( adaptTempLastStep > 0 )
04936         iout << iINFO << "      ADAPTIVE TEMPERING WILL BE DONE FROM STEP " << adaptTempFirstStep  << " TO " << adaptTempLastStep << "\n";
04937      else
04938         iout << iINFO << "      ADAPTIVE TEMPERING WILL BE DONE FROM STEP " << adaptTempFirstStep << "\n";
04939      if ( adaptTempLangevin )
04940         iout << iINFO << "      ADAPTIVE TEMPERING COUPLED TO LANGEVIN THERMOSTAT\n";
04941      if ( adaptTempRescale )
04942         iout << iINFO << "      ADAPTIVE TEMPERING COUPLED TO VELOCITY RESCALING\n";
04943      if (adaptTempRestartFreq > 0) {
04944         iout << iINFO << "      WRITING RESTART INFORMATION TO " << adaptTempRestartFile << " EVERY " << adaptTempRestartFreq << " STEPS\n";
04945      }
04946         
04947    }
04948 
04949    if (FMAOn)
04950    {
04951      iout << iINFO << "FMA ACTIVE\n";
04952      iout << iINFO << "FMA THETA              "
04953         << fmaTheta << "\n";
04954      iout << endi;
04955    }
04956 
04957    FFTWWisdomString = 0;
04958    if (PMEOn)
04959    {
04960      iout << iINFO << "PARTICLE MESH EWALD (PME) ACTIVE\n";
04961      iout << iINFO << "PME TOLERANCE               "
04962         << PMETolerance << "\n";
04963      iout << iINFO << "PME EWALD COEFFICIENT       "
04964         << PMEEwaldCoefficient << "\n";
04965      iout << iINFO << "PME INTERPOLATION ORDER     "
04966         << PMEInterpOrder << "\n";
04967      iout << iINFO << "PME GRID DIMENSIONS         "
04968         << PMEGridSizeX << " "
04969         << PMEGridSizeY << " "
04970         << PMEGridSizeZ << "\n";
04971      iout << iINFO << "PME MAXIMUM GRID SPACING    "
04972         << PMEGridSpacing << "\n";
04973      if ( PMEBarrier ) {
04974        iout << iINFO << "PME BARRIER ENABLED\n";
04975      }
04976      iout << endi;
04977      if ( useDPME ) iout << iINFO << "USING OLD DPME CODE\n";
04978 #ifdef NAMD_FFTW
04979      else if ( FFTWUseWisdom ) {  // handle FFTW wisdom
04980 #ifdef NAMD_FFTW_3
04981        iout << iINFO << "Attempting to read FFTW data from system" <<"\n" <<endi;
04982        fftwf_import_system_wisdom();
04983 #endif
04984        if (! opts.defined("FFTWWisdomFile")) {
04985          strcpy(FFTWWisdomFile,"FFTW_NAMD_");
04986          strcat(FFTWWisdomFile,NAMD_VERSION);
04987          strcat(FFTWWisdomFile,"_");
04988          strcat(FFTWWisdomFile,NAMD_PLATFORM);
04989 #ifdef NAMD_FFTW_3
04990          strcat(FFTWWisdomFile,"_FFTW3");
04991 #endif
04992          strcat(FFTWWisdomFile,".txt");
04993        }
04994 
04995        iout << iINFO << "Attempting to read FFTW data from "
04996                 << FFTWWisdomFile << "\n" << endi;
04997        FILE *wisdom_file = fopen(FFTWWisdomFile,"r");
04998        if ( wisdom_file ) {
04999 #ifdef NAMD_FFTW_3
05000          fftwf_import_wisdom_from_file(wisdom_file);
05001 #else
05002          fftw_import_wisdom_from_file(wisdom_file);
05003 #endif
05004          fclose(wisdom_file);
05005        }
05006        int nrp = 1;
05007 
05008        // rules based on work available
05009        int minslices = PMEMinSlices;
05010        int dimx = PMEGridSizeX;
05011        int nrpx = ( dimx + minslices - 1 ) / minslices;
05012        if ( nrpx > nrp ) nrp = nrpx;
05013        int dimy = PMEGridSizeY;
05014        int nrpy = ( dimy + minslices - 1 ) / minslices;
05015        if ( nrpy > nrp ) nrp = nrpy;
05016 
05017        // rules based on processors available
05018        int nrpp = CkNumPes();
05019        // if ( nrpp > 32 ) nrpp = 32;  // cap to limit messages
05020        if ( nrpp < nrp ) nrp = nrpp;
05021 
05022        // user override
05023        int nrps = PMEProcessors;
05024        if ( nrps > CkNumPes() ) nrps = CkNumPes();
05025        if ( nrps > 0 ) nrp = nrps;
05026 
05027        // make sure there aren't any totally empty processors
05028        int bx = ( dimx + nrp - 1 ) / nrp;
05029        int nrpbx = ( dimx + bx - 1 ) / bx;
05030        int by = ( dimy + nrp - 1 ) / nrp;
05031        int nrpby = ( dimy + by - 1 ) / by;
05032        nrp = ( nrpby > nrpbx ? nrpby : nrpbx );
05033        if ( bx != ( dimx + nrp - 1 ) / nrp )
05034          NAMD_bug("Error in selecting number of PME processors.");
05035        if ( by != ( dimy + nrp - 1 ) / nrp )
05036          NAMD_bug("Error in selecting number of PME processors.");
05037 
05038        // numGridPes = nrpbx;
05039        // numTransPes = nrpby;
05040        // numRecipPes = nrp;
05041        int block2 = (PMEGridSizeY + nrp - 1) / nrp;
05042        int block2_min = PMEGridSizeY % block2;
05043        if ( ! block2_min ) block2_min = block2;
05044        int dim3 = 2 * (PMEGridSizeZ/2 + 1);
05045 
05046        int n[3]; n[0] = PMEGridSizeX; n[1] = PMEGridSizeY; n[2] = PMEGridSizeZ;
05047        fftw_complex *work = new fftw_complex[n[0]];
05048        float *grid1 = (float *) fftwf_malloc(sizeof(float) *n[1]*dim3);
05049        float *grid2 = (float *) fftwf_malloc(sizeof(float) *n[0]*block2*dim3*2);
05050        iout << iINFO << "Optimizing 6 FFT steps.  1..." << endi;
05051 #ifdef NAMD_FFTW_3
05052        int fftwFlags = FFTWPatient ? FFTW_PATIENT  : FFTWEstimate ? FFTW_ESTIMATE  : FFTW_MEASURE ;
05053        int planLineSizes[1];
05054        planLineSizes[0]=n[2];
05055        int nx= n[0];
05056        int ny=block2;
05057        int sizeLines=nx*ny;
05058        int zdim = dim3;
05059        int nz=zdim;
05060        int xStride=block2 * dim3 / 2;
05061        fftwf_destroy_plan(
05062                          fftwf_plan_many_dft_r2c(1, planLineSizes, sizeLines,
05063                                                  (float *) grid2, NULL, 1, 
05064                                                  dim3,
05065                                                  (fftwf_complex *) grid2, 
05066                                                  NULL, 1,
05067                                                  dim3/2,
05068                                                  fftwFlags));
05069 
05070        iout << " 2..." << endi;
05071        fftwf_destroy_plan(
05072                          fftwf_plan_many_dft_c2r(1, planLineSizes, sizeLines,
05073                                                  (fftwf_complex *) grid2,
05074                                                  NULL, 1, 
05075                                                  dim3/2,
05076                                                  (float *) grid2, NULL, 1,
05077                                                  dim3,
05078                                                  fftwFlags));
05079        iout << " 3..." << endi;
05080        sizeLines=nz;
05081        planLineSizes[0]=block2;
05082        fftwf_destroy_plan(fftwf_plan_many_dft(1, planLineSizes, sizeLines, 
05083                                               (fftwf_complex *) grid2, NULL, 
05084                                               sizeLines, 1,
05085                                               (fftwf_complex *) grid2, NULL, 
05086                                               sizeLines, 1,
05087                                               FFTW_FORWARD, 
05088                                               fftwFlags));
05089        iout << " 4..." << endi;
05090        fftwf_destroy_plan(fftwf_plan_many_dft(1, planLineSizes, sizeLines, 
05091                                               (fftwf_complex *) grid2, NULL, 
05092                                               sizeLines, 1,
05093                                               (fftwf_complex *) grid2, NULL, 
05094                                               sizeLines, 1,
05095                                               FFTW_FORWARD, 
05096                                               fftwFlags));
05097        iout << " 5..." << endi;
05098        sizeLines=ny*nz;
05099        planLineSizes[0]=n[0];
05100        fftwf_destroy_plan(fftwf_plan_many_dft(1, planLineSizes, sizeLines,
05101                                               (fftwf_complex *) grid2, NULL,
05102                                               sizeLines, 1,
05103                                               (fftwf_complex *) grid2, NULL, 
05104                                               sizeLines, 1,
05105                                               FFTW_FORWARD,
05106                                               fftwFlags));
05107        iout << " 6..." << endi;
05108        fftwf_destroy_plan(fftwf_plan_many_dft(1, planLineSizes, sizeLines,
05109                                               (fftwf_complex *) grid2, NULL, 
05110                                               sizeLines, 1,
05111                                               (fftwf_complex *) grid2, NULL, 
05112                                               sizeLines, 1,
05113                                               FFTW_BACKWARD,
05114                                               fftwFlags));
05115 
05116 #else
05117        rfftwnd_destroy_plan( rfftwnd_create_plan_specific(
05118          2, n+1, FFTW_REAL_TO_COMPLEX,
05119          ( FFTWEstimate ? FFTW_ESTIMATE : FFTW_MEASURE )
05120          | FFTW_IN_PLACE | FFTW_USE_WISDOM, grid1, 1, 0, 0) );
05121        iout << " 2..." << endi;
05122        fftw_destroy_plan( fftw_create_plan_specific(n[0], FFTW_REAL_TO_COMPLEX,
05123          ( FFTWEstimate ? FFTW_ESTIMATE : FFTW_MEASURE )
05124          | FFTW_IN_PLACE | FFTW_USE_WISDOM, (fftw_complex *) grid2,
05125          block2*dim3/2, work, 1) );
05126        iout << " 3..." << endi;
05127        fftw_destroy_plan( fftw_create_plan_specific(n[0], FFTW_REAL_TO_COMPLEX,
05128          ( FFTWEstimate ? FFTW_ESTIMATE : FFTW_MEASURE )
05129          | FFTW_IN_PLACE | FFTW_USE_WISDOM, (fftw_complex *) grid2,
05130          block2_min*dim3/2, work, 1) );
05131        iout << " 4..." << endi;
05132        fftw_destroy_plan( fftw_create_plan_specific(n[0], FFTW_COMPLEX_TO_REAL,
05133          ( FFTWEstimate ? FFTW_ESTIMATE : FFTW_MEASURE )
05134          | FFTW_IN_PLACE | FFTW_USE_WISDOM, (fftw_complex *) grid2,
05135          block2*dim3/2, work, 1) );
05136        iout << " 5..." << endi;
05137        fftw_destroy_plan( fftw_create_plan_specific(n[0], FFTW_COMPLEX_TO_REAL,
05138          ( FFTWEstimate ? FFTW_ESTIMATE : FFTW_MEASURE )
05139          | FFTW_IN_PLACE | FFTW_USE_WISDOM, (fftw_complex *) grid2,
05140          block2_min*dim3/2, work, 1) );
05141        iout << " 6..." << endi;
05142        rfftwnd_destroy_plan( rfftwnd_create_plan_specific(
05143          2, n+1, FFTW_COMPLEX_TO_REAL,
05144          ( FFTWEstimate ? FFTW_ESTIMATE : FFTW_MEASURE )
05145          | FFTW_IN_PLACE | FFTW_USE_WISDOM, grid1, 1, 0, 0) );
05146 #endif
05147        iout << "   Done.\n" << endi;
05148        delete [] work;
05149        fftwf_free(grid1);
05150        fftwf_free(grid2);
05151 
05152        iout << iINFO << "Writing FFTW data to "
05153                 << FFTWWisdomFile << "\n" << endi;
05154        wisdom_file = fopen(FFTWWisdomFile,"w");
05155        if ( wisdom_file ) {
05156 #ifdef NAMD_FFTW_3 
05157          fftwf_export_wisdom_to_file(wisdom_file);
05158 #else
05159          fftw_export_wisdom_to_file(wisdom_file);
05160 #endif
05161          fclose(wisdom_file);
05162        }
05163 
05164 #ifdef NAMD_FFTW_3 
05165        FFTWWisdomString = fftwf_export_wisdom_to_string();
05166 #else
05167        FFTWWisdomString = fftw_export_wisdom_to_string();
05168 #endif
05169      }
05170 #endif
05171      iout << endi;
05172    }
05173    if (fullDirectOn)
05174    {
05175      iout << iINFO << "DIRECT FULL ELECTROSTATIC CALCULATIONS ACTIVE\n";
05176      iout << endi;
05177    }
05178 
05179    // MSM configure
05180    if (MSMOn)
05181    {
05182      // check MSMQuality
05183      enum { LO=0, MED, MEDHI, HI };
05184 
05185      // MSMApprox
05186      enum { CUBIC=0, QUINTIC, QUINTIC2,
05187        SEPTIC, SEPTIC3, NONIC, NONIC4, C1HERMITE, NUM_APPROX };
05188 
05189      // MSMSplit
05190      enum { TAYLOR2=0, TAYLOR3, TAYLOR4,
05191        TAYLOR5, TAYLOR6, TAYLOR7, TAYLOR8, NUM_SPLIT };
05192 
05193      if (MSMApprox || MSMSplit) {  // take these definitions
05194        if (MSMApprox < 0 || MSMApprox >= NUM_APPROX) {
05195          NAMD_die("MSM: unknown approximation requested (MSMApprox)");
05196        }
05197        if (MSMSplit < 0 || MSMSplit >= NUM_SPLIT) {
05198          NAMD_die("MSM: unknown splitting requested (MSMSplit)");
05199        }
05200      }
05201      else {  // otherwise use MSMQuality to set MSMApprox and MSMSplit
05202        switch (MSMQuality) {
05203          case LO:
05204            MSMApprox = CUBIC;
05205            MSMSplit = TAYLOR2;
05206            break;
05207          case MED:
05208            MSMApprox = QUINTIC;
05209            MSMSplit = TAYLOR3;
05210            break;
05211          case MEDHI:
05212            MSMApprox = SEPTIC;
05213            MSMSplit = TAYLOR4;
05214            break;
05215          case HI:
05216            MSMApprox = NONIC;
05217            MSMSplit = TAYLOR5;
05218            break;
05219          default:
05220            NAMD_die("MSM: unknown quality requested (MSMQuality)");
05221        }
05222      }
05223 
05224      iout << iINFO
05225        << "MULTILEVEL SUMMATION METHOD (MSM) FOR ELECTROSTATICS ACTIVE\n";
05226      if (MsmSerialOn) {
05227        iout << iINFO
05228          << "PERFORMING SERIAL MSM CALCULATION FOR LONG-RANGE PART\n";
05229      }
05230      char *approx_str, *split_str;
05231      switch (MSMApprox) {
05232        case CUBIC:    approx_str = "C1 CUBIC";   break;
05233        case QUINTIC:  approx_str = "C1 QUINTIC"; break;
05234        case QUINTIC2: approx_str = "C2 QUINTIC"; break;
05235        case SEPTIC:   approx_str = "C1 SEPTIC";  break;
05236        case SEPTIC3:  approx_str = "C3 SEPTIC";  break;
05237        case NONIC:    approx_str = "C1 NONIC";   break;
05238        case NONIC4:   approx_str = "C4 NONIC";   break;
05239        default:       approx_str = "UNKNOWN";    break;
05240      }
05241      switch (MSMSplit) {
05242        case TAYLOR2:  split_str = "C2 TAYLOR";   break;
05243        case TAYLOR3:  split_str = "C3 TAYLOR";   break;
05244        case TAYLOR4:  split_str = "C4 TAYLOR";   break;
05245        case TAYLOR5:  split_str = "C5 TAYLOR";   break;
05246        case TAYLOR6:  split_str = "C6 TAYLOR";   break;
05247        case TAYLOR7:  split_str = "C7 TAYLOR";   break;
05248        case TAYLOR8:  split_str = "C8 TAYLOR";   break;
05249        default:       split_str = "UNKNOWN";     break;
05250      }
05251      iout << iINFO
05252        << "MSM WITH " << approx_str << " INTERPOLATION "
05253        << "AND " << split_str << " SPLITTING\n"
05254        << endi;
05255 
05256    } // end MSM configure
05257    if (FMMOn)
05258    {
05259 #ifdef FMM_SOLVER
05260      iout << iINFO << "FAST MULTIPOLE METHOD (FMM) FOR ELECTROSTATICS ACTIVE\n";
05261      iout << iINFO << "PERFORMING SERIAL FMM CALCULATION\n";
05262      iout << iINFO << "FMM LEVELS = " << FMMLevels << "\n";
05263      iout << iINFO << "FMM PADDING = " << FMMPadding << " ANGSTROMS\n";
05264      iout << endi;
05265 #else
05266      NAMD_die("Must link to FMM library to use FMM\n");
05267 #endif
05268    }
05269 
05270    if ( FMAOn || PMEOn || MSMOn || fullDirectOn || GBISOn || FMMOn )
05271    {
05272      iout << iINFO << "FULL ELECTROSTATIC EVALUATION FREQUENCY      "
05273         << fullElectFrequency << "\n";
05274      iout << endi;
05275 
05276      if ( ( outputEnergies % fullElectFrequency ) &&
05277           ( fullElectFrequency % outputEnergies ) )
05278         NAMD_die("Either outputEnergies must be a multiple of fullElectFrequency or vice versa.\n");
05279    }
05280 
05281   if (MTSAlgorithm == NAIVE)
05282   {
05283     iout << iINFO << "USING NAIVE (CONSTANT FORCE) MTS SCHEME.\n" << endi;
05284   }
05285   if (MTSAlgorithm == VERLETI )
05286   {
05287     iout << iINFO << "USING VERLET I (r-RESPA) MTS SCHEME.\n" << endi;
05288   }
05289 
05290    if (longSplitting == SHARP)
05291   iout << iINFO << "SHARP SPLITTING OF LONG RANGE ELECTROSTATICS\n";
05292    else if (longSplitting == XPLOR)
05293   iout << iINFO << "XPLOR SPLITTING OF LONG RANGE ELECTROSTATICS\n";
05294    else if (longSplitting == C1)
05295   iout << iINFO << "C1 SPLITTING OF LONG RANGE ELECTROSTATICS\n";
05296    else if (longSplitting == C2)
05297   iout << iINFO << "C2 SPLITTING OF LONG RANGE ELECTROSTATICS\n";
05298 
05299    if (splitPatch == SPLIT_PATCH_POSITION)
05300   iout << iINFO << "PLACING ATOMS IN PATCHES BY POSITION\n";
05301    else if (splitPatch == SPLIT_PATCH_HYDROGEN)
05302   iout << iINFO << "PLACING ATOMS IN PATCHES BY HYDROGEN GROUPS\n";
05303 
05304    iout << endi;
05305 
05306    if (mollyOn)
05307    {
05308      iout << iINFO << "SLOW FORCE MOLLIFICATION : \n";
05309      iout << iINFO << "         ERROR TOLERANCE : " << mollyTol << "\n";
05310      iout << iINFO << "          MAX ITERATIONS : " << mollyIter << "\n";
05311      iout << endi;
05312    }
05313 
05314    if (rigidBonds != RIGID_NONE)
05315    {
05316      iout << iINFO << "RIGID BONDS TO HYDROGEN : ";
05317      if (rigidBonds == RIGID_ALL)    iout << "ALL\n";
05318      if (rigidBonds == RIGID_WATER)  iout << "WATER\n";
05319      iout << iINFO << "        ERROR TOLERANCE : " << rigidTol << "\n";
05320      iout << iINFO << "         MAX ITERATIONS : " << rigidIter << "\n";
05321      if (useSettle) iout << iINFO << "RIGID WATER USING SETTLE ALGORITHM\n";
05322      iout << endi;
05323    }
05324    
05325 
05326    if (nonbondedFrequency != 1)
05327    {
05328      iout << iINFO << "NONBONDED FORCES EVALUATED EVERY " << nonbondedFrequency << " STEPS\n";
05329    }
05330 
05331    iout << iINFO << "RANDOM NUMBER SEED     "
05332       << randomSeed << "\n";
05333 
05334    iout << endi;
05335 
05336    iout << iINFO << "USE HYDROGEN BONDS?    ";
05337    if (HydrogenBonds)
05338    {
05339   iout << "YES\n" << endi;
05340   iout << iINFO << "USE ANTECEDENT ATOMS?  ";
05341   iout << (useAntecedent ? "YES" : "NO");
05342         iout << "\nHB DIST CUT, ON, OFF   ";
05343   iout << daCutoffDist << " , " << daOnDist << " , " << daOffDist;
05344         iout << "\nHB ANGLE CUT, ON, OFF  ";
05345   iout << dhaCutoffAngle << " , " << dhaOnAngle << " , ";
05346   iout << dhaOffAngle;
05347         iout << "\nHB ATT, REP exponents  ";
05348   iout << distAttExp << " , " << distRepExp;
05349         iout << "\nHB AA, HA exponents    ";
05350   iout << aaAngleExp << " , " << haAngleExp;
05351   iout << "\n" << endi;
05352    }
05353    else
05354    {
05355   iout << "NO\n" << endi;
05356    }
05357 
05358 // If this is AMBER, then print AMBER options
05359 
05360    if (amberOn)
05361    { iout << iINFO << "Using AMBER format force field!\n";
05362      current = config->find("parmfile");
05363      iout << iINFO << "AMBER PARM FILE        " << current->data << '\n';
05364      if (opts.defined("coordinates"))
05365      { current = config->find("coordinates");
05366        iout << iINFO << "COORDINATE PDB         " << current->data << '\n';
05367      }
05368      else
05369      { current = config->find("ambercoor");
05370        iout << iINFO << "AMBER COORDINATE FILE  " << current->data << '\n';
05371      }
05372      if (readExclusions)
05373        iout << iINFO << "Exclusions will be read from PARM file!\n";
05374      else
05375        iout << iINFO << "Exclusions in PARM file will be ignored!\n";
05376      iout << iINFO << "SCNB (VDW SCALING)     " << vdwscale14 << "\n" << endi;
05377    }
05378    else if(gromacsOn)
05379    {
05380      iout << iINFO << "Using GROMACS format force field!\n";
05381 
05382      current = config->find("grotopfile");
05383      // it should be defined, but, just in case...
05384      if (current == NULL)
05385        NAMD_die("no GROMACS topology file defined!?");
05386      iout << iINFO << "GROMACS TOPO FILE        " << current->data << '\n';
05387 
05388      // XXX handle the two types of coordinates more gracefully
05389      current = config->find("grocoorfile");
05390      if (current == NULL) {
05391        current = config->find("coordinates");
05392        if (current == NULL) {
05393          NAMD_die("no coordinate file defined!?");
05394        }
05395      }
05396      iout << iINFO << "GROMACS COOR FILE        " << current->data << '\n' 
05397           << endi;
05398 
05399    }
05400    else {
05401      if ( !usePluginIO ) {
05402        if ( current = config->find("coordinates") )
05403        iout << iINFO << "COORDINATE PDB         " << current->data << '\n' << endi;
05404      }
05405 
05406      current = config->find("structure");
05407 
05408      iout << iINFO << "STRUCTURE FILE         " 
05409         << current->data << "\n" << endi;
05410 
05411      if (cosAngles) 
05412      {
05413        iout << iINFO << "COSANGLES ON. SOME ANGLES WILL BE COSINE-BASED\n" << endi;
05414      }
05415 
05416      //****** BEGIN CHARMM/XPLOR type changes
05417      if (paraTypeXplorOn)
05418      {
05419        iout << iINFO << "PARAMETER file: XPLOR format! (default) \n" << endi;
05420      }
05421      else if (paraTypeCharmmOn)
05422      {
05423        iout << iINFO << "PARAMETER file: CHARMM format! \n" << endi;
05424      }
05425      //****** END CHARMM/XPLOR type changes
05426 
05427      current = config->find("parameters");
05428 
05429      while (current != NULL)
05430      {
05431        iout << iINFO << "PARAMETERS             " 
05432           << current->data << "\n" << endi;
05433        current = current->next;
05434      }
05435    }
05436 
05437      iout << iINFO << "USING " <<
05438         ( vdwGeometricSigma ? "GEOMETRIC" : "ARITHMETIC" ) <<
05439         " MEAN TO COMBINE L-J SIGMA PARAMETERS\n" << endi;
05440 
05441    if (opts.defined("bincoordinates"))
05442    {
05443      current = config->find("bincoordinates");
05444 
05445      iout << iINFO << "BINARY COORDINATES     " 
05446               << current->data << "\n";
05447    }
05448 
05449 #ifdef MEM_OPT_VERSION
05450    if (opts.defined("binrefcoords"))
05451    {
05452      current = config->find("binrefcoords");
05453 
05454      iout << iINFO << "BINARY REF COORDS      " 
05455               << current->data << "\n";
05456    }
05457 #endif
05458 
05459    if (firstTimestep)
05460    {
05461   iout << iINFO << "FIRST TIMESTEP         "
05462      << firstTimestep << "\n" << endi;
05463    }
05464 }
05465 /*    END OF FUNCTION initialize_config_data    */
05466    
05467    
05468 /****************************************************************/
05469 /*                                                              */
05470 /*      FUNCTION parse_mgrid_params                             */
05471 /*                                                              */
05472 /*                                                              */
05473 /****************************************************************/
05474 void SimParameters::parse_mgrid_params(ConfigList *config)
05475 {
05476   StringList *current;
05477 
05478   mgridforcelist.clear();
05479   char *key = new char[81];
05480   char *valstr = new char[256];
05481   // If the old gridforce commands are still in use, parse them too.
05482   if (gridforceOn) {
05483     mgridforceOn = TRUE;
05484     char *default_key = MGRIDFORCEPARAMS_DEFAULTKEY;
05485     MGridforceParams* mgfp = NULL;
05486     mgfp = mgridforcelist.find_key(default_key);
05487     if (mgfp != NULL) {
05488       iout << iINFO << "MGRIDFORCEPOTFILE key " 
05489         << key << " redefined for file " << valstr << "\n" << endi;
05490     } else {
05491       mgfp = mgridforcelist.add(default_key);
05492     }
05493     mgfp->gridforceVolts = gridforceVolts;
05494     mgfp->gridforceScale = gridforceScale;
05495 
05496     parse_mgrid_string_param(config,"gridforcefile",&(mgfp->gridforceFile));
05497     parse_mgrid_string_param(config,"gridforcecol",&(mgfp->gridforceCol));
05498     parse_mgrid_string_param(config,"gridforcechargecol",&(mgfp->gridforceQcol));
05499     parse_mgrid_string_param(config,"gridforcepotfile",&(mgfp->gridforceVfile));
05500     
05501     mgfp->gridforceCont[0] = gridforceContA1;
05502     mgfp->gridforceCont[1] = gridforceContA2;
05503     mgfp->gridforceCont[2] = gridforceContA3;
05504     mgfp->gridforceVOffset = gridforceVOffset;
05505     
05506     mgfp->gridforceLite = gridforceLite;
05507   }
05508   
05509   // Create multigrid parameter structures
05510   current = config->find("mgridforcepotfile");
05511   while (current != NULL) {
05512     int curlen = strlen(current->data);
05513     //    iout << iINFO << "MGRIDFORCEPOTFILE " << current->data 
05514     //         << " " << curlen << "\n"  << endi;
05515     sscanf(current->data,"%80s%255s",key,valstr);
05516     
05517     MGridforceParams* mgfp = NULL;
05518     mgfp = mgridforcelist.find_key(key);
05519     if ( mgfp != NULL) {
05520       iout << iINFO << "MGRIDFORCEPOTFILE key " 
05521         << key << " redefined for file " << valstr << "\n" << endi;
05522     } else {
05523       mgfp = mgridforcelist.add(key);
05524     }
05525     int fnamelen = strlen(valstr);
05526     mgfp->gridforceVfile = new char[fnamelen+1];
05527     strncpy(mgfp->gridforceVfile,valstr,fnamelen+1);
05528     mgfp->gridforceScale.x = 
05529       mgfp->gridforceScale.y = 
05530         mgfp->gridforceScale.z = 1.;
05531     mgfp->gridforceVOffset.x = 
05532       mgfp->gridforceVOffset.y = 
05533         mgfp->gridforceVOffset.z = 0.;
05534     
05535     current = current->next;
05536   }
05537 
05538   current = config->find("mgridforcefile");
05539   while (current != NULL) {
05540     int curlen = strlen(current->data);
05541     //    iout << iINFO << "MGRIDFORCEFILE " << current->data          
05542     //         << " " << curlen << "\n"  << endi;
05543     sscanf(current->data,"%80s%255s",key,valstr);
05544     
05545     MGridforceParams* mgfp = NULL;
05546     mgfp = mgridforcelist.find_key(key);
05547     if ( mgfp == NULL) {
05548       iout << iINFO << "MGRIDFORCEFILE no key " 
05549       << key << " defined for file " << valstr << "\n" << endi;
05550     } else {
05551       int fnamelen = strlen(valstr);
05552       if (mgfp->gridforceFile != NULL) {
05553         delete [] mgfp->gridforceFile;
05554       }
05555       mgfp->gridforceFile = new char[fnamelen+1];
05556       strncpy(mgfp->gridforceFile,valstr,fnamelen+1);
05557     }
05558 
05559     current = current->next;
05560   }
05561   
05562   current = config->find("mgridforcevolts");
05563   while (current != NULL) {
05564     //    iout << iINFO << "MGRIDFORCEVOLTS " << current->data << "\n"  
05565     //         << endi;
05566     int curlen = strlen(current->data);
05567     sscanf(current->data,"%80s%255s",key,valstr);
05568     
05569     MGridforceParams* mgfp = NULL;
05570     mgfp = mgridforcelist.find_key(key);
05571     if ( mgfp == NULL) {
05572       iout << iINFO << "MGRIDFORCEVOLTS no key " 
05573       << key << " defined for file " << valstr << "\n" << endi;
05574     } else {
05575       int boolval = MGridforceParamsList::atoBool(valstr);
05576       if (boolval == -1) {
05577         iout << iINFO << "MGRIDFORCEVOLTS  key " 
05578           << key << " boolval " << valstr << " badly defined" << endi;
05579       } else {
05580         mgfp->gridforceVolts = (boolval == 1);
05581       }
05582     }
05583 
05584     current = current->next;
05585   }
05586   
05587   current = config->find("mgridforcescale");
05588   while (current != NULL) {
05589     //    iout << iINFO << "MGRIDFORCESCALE " << current->data 
05590     //         << "\n"  << endi;
05591     int curlen = strlen(current->data);
05592     int nread;
05593     sscanf(current->data,"%80s%n",key,&nread);
05594     char *val = current->data + nread + 1;
05595     
05596     MGridforceParams* mgfp = NULL;
05597     mgfp = mgridforcelist.find_key(key);
05598     if ( mgfp == NULL) {
05599       iout << iINFO << "MGRIDFORCESCALE no key " 
05600       << key << " defined for vector " << val << "\n" << endi;
05601     } else {
05602       mgfp->gridforceScale.set(val);
05603     }
05604     
05605     current = current->next;
05606   }
05607   
05608   current = config->find("mgridforcevoff");
05609   while (current != NULL) {
05610     //    iout << iINFO << "MGRIDFORCEVOFF " << current->data 
05611     //         << "\n"  << endi;
05612     int curlen = strlen(current->data);
05613     int nread;
05614     sscanf(current->data,"%80s%n",key,&nread);
05615     char *val = current->data + nread + 1;
05616     
05617     MGridforceParams* mgfp = NULL;
05618     mgfp = mgridforcelist.find_key(key);
05619     if ( mgfp == NULL) {
05620       iout << iINFO << "MGRIDFORCEVOFF no key " 
05621       << key << " defined for vector " << val << "\n" << endi;
05622     } else {
05623       mgfp->gridforceVOffset.set(val);
05624     }
05625     
05626     current = current->next;
05627   }
05628   
05629   current = config->find("mgridforcecol");
05630   while (current != NULL) {
05631     //    iout << iINFO << "MGRIDFORCECOL " << current->data 
05632     //         << "\n"  << endi;
05633     int curlen = strlen(current->data);
05634     sscanf(current->data,"%80s%255s",key,valstr);
05635     
05636     MGridforceParams* mgfp = NULL;
05637     mgfp = mgridforcelist.find_key(key);
05638     if ( mgfp == NULL) {
05639       iout << iINFO << "MGRIDFORCECOL no key " 
05640       << key << " defined for file " << valstr << "\n" << endi;
05641     } else {
05642       int collen = strlen(valstr);
05643       if (mgfp->gridforceCol != NULL) {
05644         delete [] mgfp->gridforceCol;
05645       }
05646       mgfp->gridforceCol = new char[collen+1];
05647       strncpy(mgfp->gridforceCol,valstr,collen+1);
05648      }
05649     
05650     current = current->next;
05651   }
05652   
05653   current = config->find("mgridforcechargecol");
05654   while (current != NULL) {
05655     //    iout << iINFO << "MGRIDFORCECHARGECOL " << current->data << "\n"  
05656     //         << endi;
05657     int curlen = strlen(current->data);
05658     sscanf(current->data,"%80s%255s",key,valstr);
05659     
05660     MGridforceParams* mgfp = NULL;
05661     mgfp = mgridforcelist.find_key(key);
05662     if ( mgfp == NULL) {
05663       iout << iINFO << "MGRIDFORCECHARGECOL no key " 
05664       << key << " defined for file " << valstr << "\n" << endi;
05665     } else {
05666       int collen = strlen(valstr);
05667       if (mgfp->gridforceQcol != NULL) {
05668         delete [] mgfp->gridforceQcol;
05669       }
05670       mgfp->gridforceQcol = new char[collen+1];
05671       strncpy(mgfp->gridforceQcol,valstr,collen+1);
05672     }
05673     
05674     current = current->next;
05675   }
05676   
05677   current = config->find("mgridforcecont1");
05678   while (current != NULL) {
05679     //    iout << iINFO << "MGRIDFORCECONT1 " << current->data 
05680     //         << "\n"  << endi;
05681     int curlen = strlen(current->data);
05682     sscanf(current->data,"%80s%255s",key,valstr);
05683     
05684     MGridforceParams* mgfp = NULL;
05685     mgfp = mgridforcelist.find_key(key);
05686     if ( mgfp == NULL) {
05687       iout << iINFO << "MGRIDFORCECONT1 no key " 
05688       << key << " defined for file " << valstr << "\n" << endi;
05689     } else {
05690       int boolval = MGridforceParamsList::atoBool(valstr);
05691       if (boolval == -1) {
05692         iout << iINFO << "MGRIDFORCECONT1  key " 
05693         << key << " boolval " << valstr << " badly defined" << endi;
05694       } else {
05695         mgfp->gridforceCont[0] = (boolval == 1);
05696       }
05697     }
05698     
05699     current = current->next;
05700   }
05701   
05702   current = config->find("mgridforcecont2");
05703   while (current != NULL) {
05704     //    iout << iINFO << "MGRIDFORCECONT2 " << current->data 
05705     //         << "\n"  << endi;
05706     int curlen = strlen(current->data);
05707     sscanf(current->data,"%80s%255s",key,valstr);
05708     
05709     MGridforceParams* mgfp = NULL;
05710     mgfp = mgridforcelist.find_key(key);
05711     if ( mgfp == NULL) {
05712       iout << iINFO << "MGRIDFORCECONT2 no key " 
05713       << key << " defined for file " << valstr << "\n" << endi;
05714     } else {
05715       int boolval = MGridforceParamsList::atoBool(valstr);
05716       if (boolval == -1) {
05717         iout << iINFO << "MGRIDFORCECONT2  key " 
05718         << key << " boolval " << valstr << " badly defined" << endi;
05719       } else {
05720         mgfp->gridforceCont[1] = (boolval == 1);
05721       }
05722     }
05723     
05724     current = current->next;
05725   }
05726   current = config->find("mgridforcecont3");
05727   while (current != NULL) {
05728     //    iout << iINFO << "MGRIDFORCECONT3 " << current->data 
05729     //         << "\n"  << endi;
05730     int curlen = strlen(current->data);
05731     sscanf(current->data,"%80s%255s",key,valstr);
05732     
05733     MGridforceParams* mgfp = NULL;
05734     mgfp = mgridforcelist.find_key(key);
05735     if ( mgfp == NULL) {
05736       iout << iINFO << "MGRIDFORCECONT3 no key " 
05737       << key << " defined for file " << valstr << "\n" << endi;
05738       NAMD_die("MGRIDFORCE error");
05739     } else {
05740       int boolval = MGridforceParamsList::atoBool(valstr);
05741       if (boolval == -1) {
05742         iout << iINFO << "MGRIDFORCECONT3  key " 
05743         << key << " boolval " << valstr << " badly defined" << endi;
05744       } else {
05745         mgfp->gridforceCont[2] = (boolval == 1);
05746       }
05747     }
05748     
05749     current = current->next;
05750   }
05751   
05752   current = config->find("mgridforcelite");
05753   while (current != NULL) {
05754     //    iout << iINFO << "MGRIDFORCELITE " << current->data << "\n"  
05755     //         << endi;
05756     int curlen = strlen(current->data);
05757     sscanf(current->data,"%80s%255s",key,valstr);
05758     
05759     MGridforceParams* mgfp = NULL;
05760     mgfp = mgridforcelist.find_key(key);
05761     if ( mgfp == NULL) {
05762       iout << iINFO << "MGRIDFORCELITE no key " 
05763       << key << " defined for file " << valstr << "\n" << endi;
05764     } else {
05765       int boolval = MGridforceParamsList::atoBool(valstr);
05766       if (boolval == -1) {
05767         iout << iINFO << "MGRIDFORCELITE  key " 
05768           << key << " boolval " << valstr << " badly defined" << endi;
05769       } else {
05770         mgfp->gridforceLite = (boolval == 1);
05771       }
05772     }
05773 
05774     current = current->next;
05775   }
05776   
05777   delete [] valstr;
05778   delete [] key;
05779 
05780   // Fill in default values for optional items
05781   
05782   MGridforceParams* params = mgridforcelist.get_first();
05783   
05784   while (params != NULL) {
05785     if (params->gridforceFile == NULL) {
05786       char errmsg[255];
05787       sprintf(errmsg,"Value undefined for gridforceFile for key %s\n",
05788               params->gridforceKey);
05789          NAMD_die(errmsg);
05790     }
05791     if (params->gridforceCol == NULL) {
05792       char errmsg[255];
05793       sprintf(errmsg,"Value undefined for gridforceCol for key %s\n",
05794               params->gridforceKey);
05795          NAMD_die(errmsg);
05796     }
05797     if (params->gridforceLite) {
05798         for (int i = 0; i < 3; i++) {
05799             if (params->gridforceCont[i]) {
05800                 char errmsg[255];
05801                 sprintf(errmsg, "Cannot use gridforceCont and gridforceLite together for key %s\n",
05802                         params->gridforceKey);
05803                 NAMD_die(errmsg);
05804             }
05805         }
05806     }
05807     params = params->next;
05808   }
05809   
05810 }
05811 
05812 void SimParameters::parse_mgrid_string_param(ConfigList *cl,
05813                                              const char *fieldname,
05814                                              char **dest) 
05815 {
05816   StringList *vallist = cl->find(fieldname);
05817   char *val = NULL;
05818   
05819   if (vallist != NULL) {
05820     val = vallist->data;
05821   } else {
05822     return;
05823   }
05824   
05825   int len = 0;
05826   if (val == NULL) {
05827     *dest = NULL;
05828   } else {
05829     len = strlen(val);
05830     if (len == 0) {
05831       *dest = NULL;
05832     } else {
05833       *dest = new char[len+1];
05834       strncpy(*dest,val,len+1);
05835     }
05836   }
05837 }
05838    
05839 //This function is used to create directories when outputing into
05840 //multiple files, i.e. used for Parallel IO. -Chao Mei
05841 void SimParameters::create_output_directories(const char *dirname){
05842         //output files organization:
05843         //$outputFilename/$dirname/$outputproc_rank
05844         
05845         //Step 1: create $outputFilename if necessary
05846         int baselen = strlen(outputFilename);
05847         char *filename = new char[baselen+32];
05848         memset(filename, 0, baselen+32);
05849         strcpy(filename, outputFilename);
05850         if(access(filename, F_OK)!=0) {
05851                 int ret = MKDIR(filename);
05852                 if(ret!=0) {
05853                         char errmsg[512];
05854                         sprintf(errmsg, "Error in creating top-level directory %s!", filename);
05855                         NAMD_die(errmsg);
05856                 }
05857         }
05858 
05859         //Step 2: create $dirname if necessary
05860         strcat(filename, PATHSEPSTR);
05861         strcat(filename, dirname);
05862         //check if the directory exists or not  
05863         if(access(filename, F_OK)!=0) {
05864                 int ret = MKDIR(filename);
05865                 if(ret!=0) {
05866                         char errmsg[512];
05867                         sprintf(errmsg, "Error in creating middle-level directory %s!", filename);
05868                         NAMD_die(errmsg);
05869                 }
05870         }
05871 
05872         //step 3: create $outputproc_rank if necessary
05873         char tmpstr[256];
05874         for(int i=0; i<numoutputprocs; i++) {
05875                 memset(tmpstr, 0, 256);
05876                 sprintf(tmpstr, "%s%s%d", filename, PATHSEPSTR, i);
05877                 if(access(tmpstr, F_OK)!=0) {
05878                         int ret = MKDIR(tmpstr);
05879                         if(ret!=0) {
05880                                 char errmsg[512];
05881                                 sprintf(errmsg, "Error in creating last-level directory %s!", tmpstr);
05882                                 NAMD_die(errmsg);
05883                         }
05884                 }
05885         }
05886 }
05887    
05888 /****************************************************************/
05889 /*                                                              */
05890 /*      FUNCTION print_mgrid_params                             */
05891 /*                                                              */
05892 /*                                                              */
05893 /****************************************************************/
05894 #define BoolToString(b) ((b) ? "TRUE" : "FALSE")
05895   
05896 void SimParameters::print_mgrid_params()
05897 {
05898   const MGridforceParams* params = mgridforcelist.get_first();
05899   
05900   while (params != NULL) {
05901     iout << iINFO << "MGRIDFORCE key " << params->gridforceKey << "\n" << endi;
05902     iout << iINFO << "           Potfile " << params->gridforceVfile 
05903       << "\n" << endi;
05904     iout << iINFO << "           Scale " << params->gridforceScale 
05905       << "\n" << endi;
05906     iout << iINFO << "           File " << params->gridforceFile 
05907       << "\n" << endi;
05908     iout << iINFO << "           Col " << params->gridforceCol 
05909       << "\n" << endi;
05910     
05911     char *qcol_msg = "Use atom charge";
05912     if (params->gridforceQcol != NULL) {
05913       qcol_msg = params->gridforceQcol;
05914     }
05915     iout << iINFO << "           ChargeCol " << qcol_msg
05916       << "\n" << endi;
05917     iout << iINFO << "           VOffset " << params->gridforceVOffset 
05918       << "\n" << endi;
05919     iout << iINFO << "           Continuous K1 " 
05920       << BoolToString(params->gridforceCont[0]) 
05921       << "\n" << endi;
05922     iout << iINFO << "           Continuous K2 " 
05923       << BoolToString(params->gridforceCont[1]) 
05924     << "\n" << endi;
05925     iout << iINFO << "           Continuous K3 " 
05926       << BoolToString(params->gridforceCont[2]) 
05927       << "\n" << endi;
05928     iout << iINFO << "           Volts " 
05929       << BoolToString(params->gridforceVolts)
05930       << "\n" << endi;
05931     iout << iINFO << "           Gridforce-Lite " 
05932       << BoolToString(params->gridforceLite)
05933       << "\n" << endi;
05934     params = params->next;
05935   }
05936 }
05937    
05938    
05939 /****************************************************************/
05940 /*                */
05941 /*    FUNCTION send_SimParameters      */
05942 /*                */
05943 /*  This function is used by the master process to broadcast*/
05944 /*  the parameter data to all the other nodes.  It just builds  */
05945 /*  a message with all the relevant data and broadcasts it to   */
05946 /*  the other nodes.  The routine receive_SimParameters is used */
05947 /*  by all the other nodes to receive this message.    */
05948 /*                */
05949 /****************************************************************/
05950 
05951 void SimParameters::send_SimParameters(MOStream *msg)
05952 
05953 {
05954   /*MOStream *msg = com_obj->newOutputStream(ALLBUTME, SIMPARAMSTAG, BUFSIZE);
05955   if ( msg == NULL )
05956   {
05957     NAMD_die("memory allocation failed in SimParameters::send_SimParameters");
05958   }*/
05959 
05960   msg->put(sizeof(SimParameters),(char*)this);
05961   if ( FFTWWisdomString ) {
05962     int fftwlen = strlen(FFTWWisdomString) + 1;
05963     msg->put(fftwlen);
05964     msg->put(fftwlen,FFTWWisdomString);
05965   }
05966   if ( tclBCScript ) {
05967     int tcllen = strlen(tclBCScript) + 1;
05968     msg->put(tcllen);
05969     msg->put(tcllen,tclBCScript);
05970   }
05971 
05972 #ifdef MEM_OPT_VERSION
05973   int filelen = strlen(binAtomFile)+1;
05974   msg->put(filelen);
05975   msg->put(filelen, binAtomFile);
05976 
05977   filelen = strlen(binCoorFile)+1;
05978   msg->put(filelen);
05979   msg->put(filelen, binCoorFile);
05980 
05981   if(binVelFile) {
05982     filelen = strlen(binVelFile)+1;
05983     msg->put(filelen);
05984     msg->put(filelen, binVelFile);
05985   }
05986 
05987   if(binRefFile) {
05988     filelen = strlen(binRefFile)+1;
05989     msg->put(filelen);
05990     msg->put(filelen, binRefFile);
05991   }
05992 #endif
05993 
05994   mgridforcelist.pack_data(msg);
05995   
05996   msg->end();
05997 }
05998 /*    END OF FUNCITON send_SimParameters    */
05999 
06000 /****************************************************************/
06001 /*                */
06002 /*      FUNCTION receive_SimParameters    */
06003 /*                */
06004 /*  This function is used by all the child nodes to   */
06005 /*  receive the simulation parameters from the master node.  */
06006 /*                */
06007 /****************************************************************/
06008 
06009 void SimParameters::receive_SimParameters(MIStream *msg)
06010 
06011 {
06012   msg->get(sizeof(SimParameters),(char*)this);
06013   if ( FFTWWisdomString ) {
06014     int fftwlen;
06015     msg->get(fftwlen);
06016     FFTWWisdomString = new char[fftwlen];
06017     msg->get(fftwlen,FFTWWisdomString);
06018 #ifdef NAMD_FFTW
06019 #ifdef NAMD_FFTW_3
06020     fftwf_import_wisdom_from_string(FFTWWisdomString);
06021 #else
06022     fftw_import_wisdom_from_string(FFTWWisdomString);
06023 #endif
06024 #endif
06025   }
06026   if ( tclBCScript ) {
06027     int tcllen;
06028     msg->get(tcllen);
06029     tclBCScript = new char[tcllen];
06030     msg->get(tcllen,tclBCScript);
06031   }
06032 
06033 #ifdef MEM_OPT_VERSION
06034   int filelen;
06035   msg->get(filelen);
06036   binAtomFile = new char[filelen];
06037   msg->get(filelen, binAtomFile);
06038   
06039   msg->get(filelen);
06040   binCoorFile = new char[filelen];
06041   msg->get(filelen, binCoorFile);
06042 
06043   if(binVelFile) {    
06044     msg->get(filelen);
06045     binVelFile = new char[filelen];
06046     msg->get(filelen, binVelFile);
06047   }
06048 
06049   if(binRefFile) {    
06050     msg->get(filelen);
06051     binRefFile = new char[filelen];
06052     msg->get(filelen, binRefFile);
06053   }
06054 #endif
06055 
06056 
06057   // The simParameters bit copy above put illegal values in the list pointers
06058   // So this resets everything so that unpacking will work.
06059   mgridforcelist.clear();
06060   mgridforcelist.unpack_data(msg);
06061   
06062   delete msg;
06063 }
06064 /*      END OF FUNCTION receive_SimParameters  */
06065 

Generated on Thu May 23 04:07:18 2013 for NAMD by  doxygen 1.3.9.1