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

Generated on Fri May 25 04:07:17 2012 for NAMD by  doxygen 1.3.9.1