SCRIPT_REINITVELS,
SCRIPT_RESCALEVELS,
SCRIPT_RELOADCHARGES,
+ SCRIPT_RESCALESOLUTECHARGES,
SCRIPT_CHECKPOINT,
SCRIPT_REVERT,
SCRIPT_CHECKPOINT_STORE,
LocalID aid[T::size];
const int lesOn = node->simParameters->lesOn;
+ const int soluteScalingOn = node->simParameters->soluteScalingOn;
Real invLesFactor = lesOn ?
1.0/node->simParameters->lesFactor :
1.0;
+ const Real soluteScalingFactor = node->simParameters->soluteScalingFactor;
+ const Bool soluteScalingAll = node->simParameters->soluteScalingAll;
// cycle through each patch and gather all tuples
TuplePatchListIter ai(tuplePatchList);
int homepatch = aid[0].pid;
int samepatch = 1;
int has_les = lesOn && node->molecule->get_fep_type(t.atomID[0]);
+ int has_ss = soluteScalingOn && node->molecule->get_ss_type(t.atomID[0]);
for (i=1; i < T::size; i++) {
aid[i] = atomMap->localID(t.atomID[i]);
samepatch = samepatch && ( homepatch == aid[i].pid );
has_les |= lesOn && node->molecule->get_fep_type(t.atomID[i]);
+ has_ss |= soluteScalingOn && node->molecule->get_ss_type(t.atomID[i]);
}
+ if (T::size < 4 && !soluteScalingAll) has_ss = false;
if ( samepatch ) continue;
- t.scale = has_les ? invLesFactor : 1;
+ t.scale = (!has_les && !has_ss) ? 1.0 : ( has_les ? invLesFactor : soluteScalingFactor );
for (i=1; i < T::size; i++) {
homepatch = patchMap->downstream(homepatch,aid[i].pid);
}
LocalID aid[T::size];
const int lesOn = node->simParameters->lesOn;
+ const int soluteScalingOn = node->simParameters->soluteScalingOn;
Real invLesFactor = lesOn ?
1.0/node->simParameters->lesFactor :
1.0;
+ const Real soluteScalingFactor = node->simParameters->soluteScalingFactor;
+ const Bool soluteScalingAll = node->simParameters->soluteScalingAll;
// cycle through each patch and gather all tuples
TuplePatchListIter ai(tuplePatchList);
int homepatch = aid[0].pid;
int samepatch = 1;
int has_les = lesOn && node->molecule->get_fep_type(t.atomID[0]);
+ int has_ss = soluteScalingOn && node->molecule->get_ss_type(t.atomID[0]);
for (i=1; i < T::size; i++) {
aid[i] = atomMap->localID(t.atomID[i]);
samepatch = samepatch && ( homepatch == aid[i].pid );
has_les |= lesOn && node->molecule->get_fep_type(t.atomID[i]);
+ has_ss |= soluteScalingOn && node->molecule->get_ss_type(t.atomID[i]);
}
+ if (T::size < 4 && !soluteScalingAll) has_ss = false;
if ( samepatch ) continue;
- t.scale = has_les ? invLesFactor : 1;
+ t.scale = (!has_les && !has_ss) ? 1.0 : ( has_les ? invLesFactor : soluteScalingFactor );
for (i=1; i < T::size; i++) {
homepatch = patchMap->downstream(homepatch,aid[i].pid);
}
LocalID aid[T::size];
const int lesOn = node->simParameters->lesOn;
+ const int soluteScalingOn = node->simParameters->soluteScalingOn;
Real invLesFactor = lesOn ?
1.0/node->simParameters->lesFactor :
1.0;
+ const Real soluteScalingFactor = node->simParameters->soluteScalingFactor;
+ const Bool soluteScalingAll = node->simParameters->soluteScalingAll;
// cycle through each patch and gather all tuples
// There should be only one!
int homepatch = aid[0].pid;
int samepatch = 1;
int has_les = lesOn && node->molecule->get_fep_type(t.atomID[0]);
+ int has_ss = soluteScalingOn && node->molecule->get_ss_type(t.atomID[0]);
for (i=1; i < T::size; i++) {
aid[i] = atomMap->localID(t.atomID[i]);
samepatch = samepatch && ( homepatch == aid[i].pid );
has_les |= lesOn && node->molecule->get_fep_type(t.atomID[i]);
+ has_ss |= soluteScalingOn && node->molecule->get_ss_type(t.atomID[i]);
}
+ if (T::size < 4 && !soluteScalingAll) has_ss = false;
if ( samepatch ) {
- t.scale = has_les ? invLesFactor : 1;
+ t.scale = (!has_les && !has_ss) ? 1.0 : ( has_les ? invLesFactor : soluteScalingFactor );
TuplePatchElem *p;
p = tuplePatchList.find(TuplePatchElem(homepatch));
for(i=0; i < T::size; i++) {
LocalID aid[T::size];
const int lesOn = node->simParameters->lesOn;
+ const int soluteScalingOn = node->simParameters->soluteScalingOn;
Real invLesFactor = lesOn ?
1.0/node->simParameters->lesFactor :
1.0;
+ const Real soluteScalingFactor = node->simParameters->soluteScalingFactor;
+ const Bool soluteScalingAll = node->simParameters->soluteScalingAll;
// cycle through each patch and gather all tuples
// There should be only one!
int homepatch = aid[0].pid;
int samepatch = 1;
int has_les = lesOn && node->molecule->get_fep_type(t.atomID[0]);
+ int has_ss = soluteScalingOn && node->molecule->get_ss_type(t.atomID[0]);
for (i=1; i < T::size; i++) {
aid[i] = this->atomMap->localID(t.atomID[i]);
samepatch = samepatch && ( homepatch == aid[i].pid );
has_les |= lesOn && node->molecule->get_fep_type(t.atomID[i]);
+ has_ss |= soluteScalingOn && node->molecule->get_ss_type(t.atomID[i]);
}
+ if (T::size < 4 && !soluteScalingAll) has_ss = false;
if ( samepatch ) {
- t.scale = has_les ? invLesFactor : 1;
+ t.scale = (!has_les && !has_ss) ? 1.0 : ( has_les ? invLesFactor : soluteScalingFactor );
TuplePatchElem *p;
p = this->tuplePatchList.find(TuplePatchElem(homepatch));
for(i=0; i < T::size; i++) {
iout << "RESCALING VELOCITIES AT STEP " << simParams->firstTimestep
<< " BY " << simParams->scriptArg1 << "\n" << endi;
break;
+ case SCRIPT_RESCALESOLUTECHARGES:
+ // Parameter setting already reported in ScriptTcl
+ // Nothing to do!
+ break;
case SCRIPT_CHECKPOINT:
iout << "CHECKPOINTING AT STEP " << simParams->firstTimestep
<< "\n" << endi;
#include "Parameters.h"
// #define DEBUGM
#include "Debug.h"
-
+#include "Molecule.h"
//----------------------------------------------------------------------
LJTable::LJTable()
{
+ Bool soluteScalingOn = Node::Object()->simParameters->soluteScalingOn;
+
+ if (!soluteScalingOn) {
table_dim = Node::Object()->parameters->get_num_vdw_params();
+ } else {
+ int ss_dim = Node::Object()->molecule->ss_num_vdw_params;
+ table_dim = ss_dim + Node::Object()->parameters->get_num_vdw_params();
+ }
table_alloc = new char[2*table_dim*table_dim*sizeof(TableEntry) + 31];
char *table_align = table_alloc;
while ( (long)table_align % 32 ) table_align++;
SimParameters *simParams = Node::Object()->simParameters;
int useGeom = simParams->vdwGeometricSigma;
Bool tabulatedEnergies = simParams->tabulatedEnergies;
-
+ Bool soluteScalingOn = simParams->soluteScalingOn;
+ BigReal soluteScalingFactor = simParams->soluteScalingFactorVdw;
+ unsigned int table_dim_org = params->get_num_vdw_params();
+ int ss_dim = Node::Object()->molecule->ss_num_vdw_params;
Real A, B, A14, B14;
int K = -1;
+ int *ss_vdw_type = Node::Object()->molecule->ss_vdw_type;
// BigReal sigma_max;
// We need the A and B parameters for the Van der Waals. These can
// be explicitly be specified for this pair or calculated from the
Real sigma_i, sigma_i14, epsilon_i, epsilon_i14;
Real sigma_j, sigma_j14, epsilon_j, epsilon_j14;
+ if (!soluteScalingOn) {
+
params->get_vdw_params(&sigma_i, &epsilon_i, &sigma_i14,
&epsilon_i14,i);
params->get_vdw_params(&sigma_j, &epsilon_j, &sigma_j14,
&epsilon_j14,j);
-
+ } else {
+ int i_type = (i >= table_dim_org)? ss_vdw_type[i-table_dim_org]:i;
+ int j_type = (j >= table_dim_org)? ss_vdw_type[j-table_dim_org]:j;
+ params->get_vdw_params(&sigma_i, &epsilon_i, &sigma_i14,
+ &epsilon_i14,i_type);
+ params->get_vdw_params(&sigma_j, &epsilon_j, &sigma_j14,
+ &epsilon_j14,j_type);
+ }
BigReal sigma_ij =
useGeom ? sqrt(sigma_i*sigma_j) : 0.5*(sigma_i+sigma_j);
BigReal sigma_ij14 =
cur_scaled->B = 4.0 * sigma_ij14 * epsilon_ij14;
cur_scaled->A = cur_scaled->B * sigma_ij14;
+ if (soluteScalingOn) {
+ if (i >= table_dim_org && i < (table_dim_org+ss_dim) && j < table_dim_org) {
+ cur->A *= sqrt(soluteScalingFactor);
+ cur->B *= sqrt(soluteScalingFactor);
+ cur_scaled->A *= sqrt(soluteScalingFactor);
+ cur_scaled->B *= sqrt(soluteScalingFactor);
+ }
+ if (i < table_dim_org && j >= table_dim_org && j < (table_dim_org+ss_dim)) {
+ cur->A *= sqrt(soluteScalingFactor);
+ cur->B *= sqrt(soluteScalingFactor);
+ cur_scaled->A *= sqrt(soluteScalingFactor);
+ cur_scaled->B *= sqrt(soluteScalingFactor);
+ }
+ if (i >=table_dim_org && i < (table_dim_org+ss_dim) && j >= table_dim_org && j < (table_dim_org+ss_dim)) {
+ cur->A *= soluteScalingFactor;
+ cur->B *= soluteScalingFactor;
+ cur_scaled->A *= soluteScalingFactor;
+ cur_scaled->B *= soluteScalingFactor;
+ }
+ }
+
if ( tabulatedEnergies && ( cur->A < 0 || cur_scaled->A < 0 ) )
NAMD_die("LJ A is negative with tabulatedEnergies enabled");
}
//fepb
fepAtomFlags=NULL;
//fepe
-
+//soluteScaling
+ ssAtomFlags=NULL;
+ ss_vdw_type=NULL;
+ ss_index=NULL;
+//soluteScaling
nameArena = new ObjectArena<char>;
// nameArena->setAlignment(8);
// arena->setAlignment(32);
delete [] fepAtomFlags;
//fepe
+//soluteScaling
+ if (ssAtomFlags != NULL)
+ delete [] ssAtomFlags;
+ if (ss_vdw_type != NULL)
+ delete [] ss_vdw_type;
+ if (ss_index != NULL)
+ delete [] ss_index;
+//soluteScaling
+
if (qmAtomGroup != NULL)
delete [] qmAtomGroup;
}
//fepe
+ if (simParams->soluteScalingOn) {
+ msg->put(numAtoms*sizeof(char), (char*)ssAtomFlags);
+ msg->put(ss_num_vdw_params);
+ msg->put(params->get_num_vdw_params()*sizeof(int), (char*)ss_vdw_type);
+ msg->put(numAtoms*sizeof(int), (char*)ss_index);
+ }
+
#ifdef OPENATOM_VERSION
// needs to be refactored into its own openatom version
if (simParams->openatomOn ) {
}
//fepe
+//soluteScaling
+ if (simParams->soluteScalingOn) {
+ delete [] ssAtomFlags;
+ delete [] ss_vdw_type;
+ delete [] ss_index;
+ ssAtomFlags = new unsigned char[numAtoms];
+ ss_vdw_type = new int [params->get_num_vdw_params()];
+ ss_index = new int [numAtoms];
+ msg->get(numAtoms*sizeof(unsigned char), (char*)ssAtomFlags);
+ msg->get(ss_num_vdw_params);
+ msg->get(params->get_num_vdw_params()*sizeof(int), (char*)ss_vdw_type);
+ msg->get(numAtoms*sizeof(int), (char*)ss_index);
+ }
+//soluteScaling
#ifdef OPENATOM_VERSION
// This needs to be refactored into its own version
if (simParams->openatomOn) {
}
}
// End of function build_fep_flags
+
+// XXX Passing in cwd is useless, since the only caller (NamdState) always
+// sends NULL - note that several other routines have this same form,
+// which probably dates back to much earlier NAMD
+// XXX Should not be necessary to pass PDB pointer as nonconst when
+// we just want to read from it.
+//
+void Molecule::build_ss_flags(
+ const StringList *ssfile,
+ const StringList *sscol,
+ PDB *initial_pdb,
+ const char *cwd
+ ) {
+ PDB *bPDB;
+ int bcol = 4;
+ Real bval = 0;
+ int i, j;
+ char filename[129];
+
+ if (ssfile == NULL) {
+ if ( ! initial_pdb ) {
+ NAMD_die("Initial PDB file unavailable, ssFile required.");
+ }
+ bPDB = initial_pdb;
+ strcpy(filename, "coordinate PDB file (default)");
+ }
+ else {
+ if (ssfile->next != NULL) {
+ NAMD_die("Multiple definitions of ssFile in configuration file");
+ }
+
+ if ((cwd == NULL) || (ssfile->data[0] == '/')) {
+ strcpy(filename, ssfile->data);
+ }
+ else {
+ strcpy(filename, cwd);
+ strcat(filename, ssfile->data);
+ }
+
+ bPDB = new PDB(filename);
+ if (bPDB == NULL) {
+ NAMD_die("Memory allocation failed in Molecule::build_ss_flags");
+ }
+
+ if (bPDB->num_atoms() != numAtoms) {
+ NAMD_die("Number of atoms in ssFile PDB does not match coordinate PDB");
+ }
+ }
+
+ if (sscol == NULL) {
+ bcol = 4;
+ }
+ else {
+ if (sscol->next != NULL) {
+ NAMD_die("Multiple definitions of ssCol value in config file");
+ }
+
+ if (strcasecmp(sscol->data, "X") == 0) {
+ bcol = 1;
+ }
+ else if (strcasecmp(sscol->data, "Y") == 0) {
+ bcol = 2;
+ }
+ else if (strcasecmp(sscol->data, "Z") == 0) {
+ bcol = 3;
+ }
+ else if (strcasecmp(sscol->data, "O") == 0) {
+ bcol = 4;
+ }
+ else if (strcasecmp(sscol->data, "B") == 0) {
+ bcol = 5;
+ }
+ else {
+ NAMD_die("ssCol must have value of X, Y, Z, O or B");
+ }
+ }
+
+ iout << iINFO << "Reading solute scaling data from file: "
+ << filename << "\n" << endi;
+ iout << iINFO << "Reading solute scaling flags from column: "
+ << bcol << "\n" << endi;
+
+ ssAtomFlags = new unsigned char[numAtoms];
+ ss_index = new int[numAtoms];
+
+ if (ssAtomFlags == NULL || ss_index == NULL) {
+ NAMD_die("Memory allocation failed in Molecule::build_ss_params()");
+ }
+
+ num_ss = 0;
+ for (i = 0; i < numAtoms; i++) {
+ switch (bcol) {
+ case 1:
+ bval = (bPDB->atom(i))->xcoor();
+ break;
+ case 2:
+ bval = (bPDB->atom(i))->ycoor();
+ break;
+ case 3:
+ bval = (bPDB->atom(i))->zcoor();
+ break;
+ case 4:
+ bval = (bPDB->atom(i))->occupancy();
+ break;
+ case 5:
+ bval = (bPDB->atom(i))->temperaturefactor();
+ break;
+ }
+ if (simParams->soluteScalingOn) {
+ if (bval == 1.0) {
+ ssAtomFlags[i] = 1;
+ ss_index[num_ss] = i;
+ num_ss++;
+ }
+ else {
+ ssAtomFlags[i] = 0;
+ }
+ }
+ }
+
+ if (ssfile != NULL) {
+ delete bPDB;
+ }
+
+ // number of LJtypes read in from params files
+ int LJtypecount = params->get_num_vdw_params();
+
+ // generate a new array of LJtypecount elements.
+ // Each element stores number of REST2 atoms of that LJType.
+ int *numAtomsByLjType = new int[LJtypecount];
+
+ // array that stores LJTypes for REST2 atoms based on array numAtomsByLjType.
+ // The 'new' LJTypes will be used to construct extended LJTable later.
+ ss_vdw_type = new int[LJtypecount];
+
+ // zero number of REST2 atoms per LJType.
+ for (i = 0; i < LJtypecount; i++) {
+ numAtomsByLjType[i] = 0;
+ }
+
+ // count number of REST2 atoms (histogram) per LJType.
+ // The num_ss is the total number of REST2 atoms.
+ for (i = 0; i < num_ss; i++) {
+ numAtomsByLjType[atoms[ss_index[i]].vdw_type]++;
+ }
+
+ //zero number of vdw param types for REST2 atoms
+ ss_num_vdw_params = 0;
+ for (i = 0; i < LJtypecount; i++) { //loop all LJTypes.
+ // only process LJTypes that have nonzero REST2 atoms.
+ if (numAtomsByLjType[i] != 0) {
+ // Build a subset of vdw params for REST2 atoms.
+ // Each REST2 atom will have a new vdw type index
+ ss_vdw_type[ss_num_vdw_params] = i;
+ // once meets a LJType of nonzero REST2 atoms,
+ // number of vdw param types of REST2 increments.
+ ss_num_vdw_params++;
+ }
+ }
+
+ for (i = 0; i < num_ss; i++) { // loop over all REST2 atoms
+ // loop over all vdw param types of REST2 atoms
+ for (j = 0; j < ss_num_vdw_params; j++) {
+ // Extends number of LJTypes with REST2 atoms.
+ if (atoms[ss_index[i]].vdw_type == ss_vdw_type[j]) {
+ // The LJType of a REST2 atom now is equal to sum of original #LJTypes
+ // and its vdw type index within REST2 atoms (ss_vdw_type)
+ atoms[ss_index[i]].vdw_type = LJtypecount + j;
+ }
+ }
+ }
+
+ delete [] numAtomsByLjType;
+
+} // End of function build_ss_flags
+
//
//
// FUNCTION delete_alch_bonded
unsigned char *fepAtomFlags;
//fepe
+//soluteScaling
+ unsigned char *ssAtomFlags;
+ unsigned int num_ss;
+//soluteScaling
+
//occupancy and bfactor data from plugin-based IO implementation of loading structures
float *occupancy;
float *bfactor;
void read_parm(const GromacsTopFile *);
public:
+
+ // for solute scaling
+ int ss_num_vdw_params;
+ int *ss_vdw_type;
+ int *ss_index;
+
// DRUDE
int is_drude_psf; // flag for reading Drude PSF
int is_lonepairs_psf; // flag for reading lone pairs from PSF
void delete_alch_bonded(void);
//fepe
+ /**
+ * Build the flags needed for solute scaling.
+ *
+ * A PDB file is read, indicating which atoms are to be scaled,
+ * and an array is maintained marking which are to be included.
+ * Each marked atom then has its corresponding van der Waals
+ * type number reassigned to enable extending the LJTable with
+ * scaled interaction values.
+ */
+ void build_ss_flags(
+ const StringList *ssfile, ///< config "ssfile = my.pdb" for PDB filename
+ const StringList *sscol, ///< config "ssfile = O" indicating column
+ PDB *initial_pdb, ///< the initial PDB file
+ const char *cwd ///< current working directory
+ );
+
void build_exPressure_atoms(StringList *, StringList *, PDB *, char *);
// Determine which atoms are excluded from
// pressure (if any)
return(fepAtomFlags[anum]);
}
+//soluteScaling
+ unsigned char get_ss_type(int anum) const
+ {
+ return(ssAtomFlags[anum]);
+ }
+//soluteScaling
+
/* BKR - Get the FEP type (i.e. 0, 1, or 2) of a bonded _interaction_ based
on the atom indices of the atoms involved (internally converted to FEP
types) and the order of the bonded interaction (i.e. 2, 3, or 4). This
NAMD_die("MEMOPT: alchThermIntOn could not be enabled in memory optimized version");
if(simParameters->lesOn)
NAMD_die("MEMOPT: lesOn could not be enabled in memory optimized version");
+ if(simParameters->soluteScalingOn)
+ NAMD_die("MEMOPT: soluteScalingOn could not be enabled in memory optimized version");
if(simParameters->lonepairs) {
NAMD_die("MEMOPT: lonepairs could not be enabled in memory optimized version");
}
molecule->build_fep_flags(configList->find("lesfile"),
configList->find("lescol"), pdb, NULL, "les");
}
+ if (simParameters->soluteScalingOn) {
+ molecule->build_ss_flags(configList->find("ssfile"),
+ configList->find("sscol"), pdb, NULL);
+ }
if (simParameters->pairInteractionOn) {
molecule->build_fep_flags(configList->find("pairInteractionFile"),
configList->find("pairInteractionCol"), pdb, NULL, "pairInteraction");
iout << iINFO << molecule->numFepInitial <<
" LOCALLY ENHANCED ATOMS ENABLED\n";
}
+
+ if (simParameters->soluteScalingOn) {
+ iout << iINFO << " SOLUTE SCALING ENABLED\n";
+ }
if (simParameters->pairInteractionOn) {
iout << iINFO << "PAIR INTERACTION GROUP 1 CONTAINS "
int require(const char *newname, const char *parent, const char *msg,
char *ptr);
+ //
+ // XXX prototypes for optional reverse "newname" and "parent"
+ // (the argument names in a prototype don't matter
+ // but it makes the code more difficult to understand)
+ //
int optional(const char *newname, const char *parent, const char *msg,
BigReal *ptr, BigReal defalt);
int optional(const char *newname, const char *parent, const char *msg,
ScriptTcl *script = (ScriptTcl *)clientData;
script->setParameter(param,value);
+ // deal with some possible specifics
+ if ( ! strncasecmp(param,"soluteScalingFactor",MAX_SCRIPT_PARAM_SIZE) ||
+ ! strncasecmp(param,"soluteScalingFactorCharge",MAX_SCRIPT_PARAM_SIZE)) {
+ script->runController(SCRIPT_RESCALESOLUTECHARGES);
+ }
+
return TCL_OK;
}
random = new Random(simParams->randomSeed);
random->split(patch->getPatchID()+1,PatchMap::Object()->numPatches()+1);
+ // Is soluteScaling enabled?
+ if (simParams->soluteScalingOn) {
+ // If so, we must "manually" perform charge scaling on startup because
+ // Sequencer will not get a scripting task for initial charge scaling.
+ // Subsequent rescalings will take place through a scripting task.
+ rescaleSoluteCharges(simParams->soluteScalingFactorCharge);
+ }
+
rescaleVelocities_numTemps = 0;
stochRescale_count = 0;
berendsenPressure_count = 0;
case SCRIPT_RESCALEVELS:
rescaleVelocitiesByFactor(simParams->scriptArg1);
break;
+ case SCRIPT_RESCALESOLUTECHARGES:
+ rescaleSoluteCharges(simParams->soluteScalingFactorCharge);
+ break;
case SCRIPT_RELOADCHARGES:
reloadCharges();
break;
}
}
+// REST2 solute charge scaling
+void Sequencer::rescaleSoluteCharges(BigReal factor)
+{
+ FullAtom *a = patch->atom.begin();
+ int numAtoms = patch->numAtoms;
+ Molecule *molecule = Node::Object()->molecule;
+ BigReal sqrt_factor = sqrt(factor);
+ // apply scaling to the original charge (stored in molecule)
+ // of just the marked solute atoms
+ for ( int i = 0; i < numAtoms; ++i ) {
+ if (molecule->get_ss_type(a[i].id)) {
+ a[i].charge = sqrt_factor * molecule->atomcharge(a[i].id);
+ }
+ }
+}
+
void Sequencer::tcoupleVelocities(BigReal dt_fs, int step)
{
if ( simParams->tCoupleOn )
void minimizationQuenchVelocity(void);
void reloadCharges();
+ void rescaleSoluteCharges(BigReal);
BigReal adaptTempT; // adaptive tempering temperature
void adaptTempUpdate(int); // adaptive tempering temperature update
#include "InfoStream.h"
#include "ComputeNonbondedUtil.h"
+#include "LJTable.h"
#include "ComputePme.h"
#include "ConfigList.h"
#include "SimParameters.h"
ComputeNonbondedUtil::select();
return;
}
+
if ( ! strncasecmp(param,"commOnly",MAX_SCRIPT_PARAM_SIZE) ) {
commOnly = atobool(value);
ComputeNonbondedUtil::select();
return;
}
+ // REST2 - We don't have to make any changes to ComputeNonbondedUtil
+ // other than recalculating the LJTable. Skip doing the other stuff.
+ if ( ! strncasecmp(param,"soluteScalingFactor",MAX_SCRIPT_PARAM_SIZE)) {
+ soluteScalingFactor = atof(value);
+ if (soluteScalingFactor < 0.0) {
+ NAMD_die("Solute scaling factor should be non-negative\n");
+ }
+ soluteScalingFactorCharge = soluteScalingFactor;
+ soluteScalingFactorVdw = soluteScalingFactor;
+ if ( ComputeNonbondedUtil::ljTable ) delete ComputeNonbondedUtil::ljTable;
+ ComputeNonbondedUtil::ljTable = new LJTable;
+ return;
+ }
+ if ( ! strncasecmp(param,"soluteScalingFactorVdw",MAX_SCRIPT_PARAM_SIZE)) {
+ soluteScalingFactorVdw = atof(value);
+ if (soluteScalingFactorVdw < 0.0) {
+ NAMD_die("Solute scaling factor for van der Waals "
+ "should be non-negative\n");
+ }
+ if ( ComputeNonbondedUtil::ljTable ) delete ComputeNonbondedUtil::ljTable;
+ ComputeNonbondedUtil::ljTable = new LJTable;
+ return;
+ }
+ if ( ! strncasecmp(param,"soluteScalingFactorCharge",MAX_SCRIPT_PARAM_SIZE)) {
+ soluteScalingFactorCharge = atof(value);
+ if (soluteScalingFactorCharge < 0.0) {
+ NAMD_die("Solute scaling factor for electrostatics "
+ "should be non-negative\n");
+ }
+ return;
+ }
+
char *error = new char[2 * MAX_SCRIPT_PARAM_SIZE + 100];
sprintf(error,"Setting parameter %s from script failed!\n",param);
NAMD_die(error);
opts.optionalB("les", "lesReduceMass", "Reduce enhanced atom mass?",
&lesReduceMass, FALSE);
+ // REST2 (replica exchange solute tempering) parameters
+ opts.optionalB("main", "soluteScaling",
+ "Is replica exchange solute tempering enabled?",
+ &soluteScalingOn, FALSE);
+ opts.require("soluteScaling", "soluteScalingFactor",
+ "Solute scaling factor",
+ &soluteScalingFactor);
+ opts.range("soluteScalingFactor", NOT_NEGATIVE);
+ opts.optional("soluteScaling", "soluteScalingFactorCharge",
+ "Solute scaling factor for electrostatic interactions",
+ &soluteScalingFactorCharge);
+ opts.range("soluteScalingFactorVdw", NOT_NEGATIVE);
+ opts.optional("soluteScaling", "soluteScalingFactorVdw",
+ "Solute scaling factor for van der Waals interactions",
+ &soluteScalingFactorVdw);
+ opts.range("soluteScalingFactorVdw", NOT_NEGATIVE);
+ opts.optional("soluteScaling", "ssFile",
+ "PDB file with scaling flags; if undefined, defaults to main PDB file",
+ PARSE_STRING);
+ opts.optional("soluteScaling", "ssCol",
+ "Column in the ssFile providing the scaling flag",
+ PARSE_STRING);
+ opts.optionalB("main", "soluteScalingAll",
+ "Apply scaling also to bond and angle interactions?",
+ &soluteScalingAll, TRUE);
+
// Drude oscillators
opts.optionalB("main", "drude", "Perform integration of Drude oscillators?",
&drudeOn, FALSE);
randomSeed = (unsigned int) time(NULL) + 31530001 * CmiMyPartition();
}
+ // REST2
+ if (opts.defined("soluteScaling")) {
+ // Parameters soluteScalingFactorCharge and soluteScalingFactorVdw
+ // allow independent scaling of electrostatics and van der Waals.
+ // Initialize with soluteScalingFactor if either is not already set.
+ if ( ! opts.defined("soluteScalingFactorCharge") ) {
+ soluteScalingFactorCharge = soluteScalingFactor;
+ }
+ if ( ! opts.defined("soluteScalingFactorVdw") ) {
+ soluteScalingFactorVdw = soluteScalingFactor;
+ }
+ }
+
//fepb
alchFepOnAtStartup = alchFepOn = FALSE;
alchThermIntOnAtStartup = alchThermIntOn = FALSE;
if ( lesReduceMass ) iout << iINFO
<< "SCALING ENHANCED ATOM MASS BY 1/" << lesFactor << "\n";
}
+
+ // REST2
+ if ( soluteScalingOn ) {
+ iout << iINFO << "SOLUTE SCALING IS ACTIVE\n";
+ if (soluteScalingFactorCharge != soluteScalingFactorVdw) {
+ iout << iINFO << "SCALING FOR ELECTROSTATIC INTERACTIONS IS "
+ << soluteScalingFactorCharge << "\n";
+ iout << iINFO << "SCALING FOR VAN DER WAALS INTERACTIONS IS "
+ << soluteScalingFactorVdw << "\n";
+ iout << iINFO << "SCALING FOR BONDED INTERACTIONS IS "
+ << soluteScalingFactor << "\n";
+ }
+ else {
+ iout << iINFO << "SOLUTE SCALING FACTOR IS "
+ << soluteScalingFactor << "\n";
+ }
+ if ( ! soluteScalingAll ) {
+ iout << iINFO << "SOLUTE SCALING DISABLED FOR BONDS AND ANGLES\n";
+ }
+ }
if ( pairInteractionOn ) {
iout << iINFO << "PAIR INTERACTION CALCULATIONS ACTIVE\n";
Bool lesReduceTemp; // Reduce enhanced atom temperature?
Bool lesReduceMass; // Reduce enhanced atom mass?
+ // REST2
+ Bool soluteScalingOn;
+ /**< REST2 (replica exchange solute tempering) on? */
+ BigReal soluteScalingFactor;
+ /**< scaling factor for solute interactions */
+ BigReal soluteScalingFactorCharge;
+ /**< optional independent control over scaling factor for
+ * electrostatic interactions of solute */
+ BigReal soluteScalingFactorVdw;
+ /**< optional independent control over scaling factor for
+ * van der Waals interactions of solute */
+ Bool soluteScalingAll;
+ /**< enables scaling for bond and angle terms (default is off) */
+
Bool extForcesOn; // Are ext command forces present?
char extForcesCommand[256];
char extCoordFilename[128];
Code Review / namd.git / commitdiff