REST2 implementation in NAMD 88/3588/10
authorWei Jiang <weijiang731@gmail.com>
Mon, 29 Jan 2018 18:26:28 +0000 (12:26 -0600)
committerDavid Hardy <dhardy@ks.uiuc.edu>
Wed, 2 May 2018 05:50:19 +0000 (00:50 -0500)
Contributed by Wei Jiang, Jonathan Thirman and Sunhwan Jo. REST2 denotes
second generation Replica Exchange Solute Tempering, where all heated
atoms get scaled charge and van der Waals parameters, as well as dihedral
and, optionally, bond and angle parameters.

The config file parameters are:
  soluteScaling - enable REST2, off by default
  soluteScalingFactor - scaling factor for solute interactions
  soluteScalingFactorCharge - optional factor for electrostatics
  soluteScalingFactorVdw - optional factor for van der Waals
  soluteScalingAll - apply scaling to bonds and angles, on by default
  ssFile - optional PDB for marking solute atoms, defaults to main PDB
  ssCol - optional column from ssFile, defaults to Occupancy column

Parameter soluteScalingFactor is the primary scaling factor for all
solute interactions.  Electrostatics and van der Waals interactions can
also be independently scaled using soluteScalingFactorCharge and
soluteScalingFactorVdw, respectively.

Rather than scaling individual energy and force terms, the rescaling
coefficients can be directly applied to the force field parameters
during the parameter reading and setup phase.  An extended LJ table is
tabulated for the van der Waals interactions, and the charges are scaled
before integration.  Each is updated at the Tcl scripting level, so that
there is no impact to performance during time stepping.

Bonded interaction rescaling is also applied to the force field
parameters rather than the individual interactions.  The scaling is
performed in ComputeHome(Self)Tuples.h, which is then offloaded
correctly to GPU.

However, the van der Waals parameter rescaling for now works correctly
only for CPU.  Testing shows that the CUDA kernels are not receiving the
updated LJ table parameters.

The charge rescaling is working correctly for the short-range non-bonded
CUDA kernels.  However, the electrostatic potential is not rescaling
correctly unless "usePMECUDA off" is specified (testing on one node).

Change-Id: Ic51dc0036310e85a6d6172fa475034f467ca18ab

14 files changed:
src/Broadcasts.h
src/ComputeHomeTuples.h
src/ComputeSelfTuples.h
src/Controller.C
src/LJTable.C [changed mode: 0644->0755]
src/Molecule.C
src/Molecule.h
src/NamdState.C
src/ParseOptions.h
src/ScriptTcl.C
src/Sequencer.C
src/Sequencer.h
src/SimParameters.C
src/SimParameters.h

index cc644d1..3f37e08 100644 (file)
@@ -21,6 +21,7 @@ enum {
   SCRIPT_REINITVELS,
   SCRIPT_RESCALEVELS,
   SCRIPT_RELOADCHARGES,
+  SCRIPT_RESCALESOLUTECHARGES,
   SCRIPT_CHECKPOINT,
   SCRIPT_REVERT,
   SCRIPT_CHECKPOINT_STORE,
index 181115e..c508429 100644 (file)
@@ -203,9 +203,12 @@ template <class T, class S, class P> class HomeTuples : public Tuples {
       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);
@@ -252,13 +255,16 @@ template <class T, class S, class P> class HomeTuples : public Tuples {
             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);
             }
@@ -338,9 +344,12 @@ template <class T, class S, class P> class ComputeHomeTuples : public Compute {
       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);
@@ -374,13 +383,16 @@ template <class T, class S, class P> class ComputeHomeTuples : public Compute {
              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);
              }
index beeb7dd..c950dfd 100644 (file)
@@ -53,9 +53,12 @@ template <class T, class S, class P> class SelfTuples : public HomeTuples<T, S,
       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!
@@ -103,13 +106,16 @@ template <class T, class S, class P> class SelfTuples : public HomeTuples<T, S,
             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++) {
@@ -174,9 +180,12 @@ template <class T, class S, class P> class ComputeSelfTuples :
       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!
@@ -212,13 +221,16 @@ template <class T, class S, class P> class ComputeSelfTuples :
              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++) {
index 8ef5d55..5d045c9 100644 (file)
@@ -302,6 +302,10 @@ void Controller::algorithm(void)
         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;
old mode 100644 (file)
new mode 100755 (executable)
index aa2fc36..6dab772
 #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++;
@@ -50,9 +57,13 @@ void LJTable::compute_vdw_params(int i, int j,
   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
@@ -99,11 +110,20 @@ void LJTable::compute_vdw_params(int i, int j,
     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 =
@@ -125,6 +145,27 @@ void LJTable::compute_vdw_params(int i, int j,
     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");
   }
index 755ff04..79739a4 100644 (file)
@@ -299,7 +299,11 @@ void Molecule::initialize(SimParameters *simParams, Parameters *param)
 //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);
@@ -681,6 +685,15 @@ Molecule::~Molecule()
        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;
   
@@ -5537,6 +5550,13 @@ void Molecule::send_Molecule(MOStream *msg){
   }
   //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 ) {
@@ -6003,6 +6023,20 @@ void Molecule::receive_Molecule(MIStream *msg){
       }
 //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) {
@@ -8710,7 +8744,183 @@ void Molecule::build_fep_flags(StringList *alchfile, StringList *alchcol,
   }
 }
 // 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
index df7e421..dfed468 100644 (file)
@@ -342,6 +342,11 @@ private:
        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;
@@ -444,6 +449,12 @@ private:
        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
@@ -930,6 +941,22 @@ public:
         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)
@@ -1304,6 +1331,13 @@ public:
     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 
index ef7f913..6e7233d 100644 (file)
@@ -112,6 +112,8 @@ void NamdState::checkMemOptCompatibility(){
         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");
     }
@@ -504,6 +506,10 @@ int NamdState::loadStructure(const char *molFilename, const char *pdbFilename, i
            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");
@@ -682,6 +688,10 @@ int NamdState::loadStructure(const char *molFilename, const char *pdbFilename, i
            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 "
index d876548..79fd569 100644 (file)
@@ -172,6 +172,11 @@ class ParseOptions {
    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,
index 76acd76..bcc7da9 100644 (file)
@@ -963,6 +963,12 @@ int ScriptTcl::Tcl_param(ClientData clientData,
   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;
 }
 
index a37b609..67f0619 100644 (file)
@@ -71,6 +71,14 @@ Sequencer::Sequencer(HomePatch *p) :
     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;
@@ -137,6 +145,9 @@ void Sequencer::algorithm(void)
       case SCRIPT_RESCALEVELS:
        rescaleVelocitiesByFactor(simParams->scriptArg1);
        break;
+      case SCRIPT_RESCALESOLUTECHARGES:
+        rescaleSoluteCharges(simParams->soluteScalingFactorCharge);
+        break;
       case SCRIPT_RELOADCHARGES:
        reloadCharges();
        break;
@@ -1654,6 +1665,22 @@ void Sequencer::reloadCharges()
   }
 }
 
+// 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 )
index f7b86fe..9973ac4 100644 (file)
@@ -77,6 +77,7 @@ protected:
     void minimizationQuenchVelocity(void);
 
     void reloadCharges();
+    void rescaleSoluteCharges(BigReal);
 
     BigReal adaptTempT;         // adaptive tempering temperature
     void adaptTempUpdate(int); // adaptive tempering temperature update
index 8e08dfa..5ee0de6 100644 (file)
@@ -19,6 +19,7 @@
 
 #include "InfoStream.h"
 #include "ComputeNonbondedUtil.h"
+#include "LJTable.h"
 #include "ComputePme.h"
 #include "ConfigList.h"
 #include "SimParameters.h"
@@ -370,12 +371,45 @@ void SimParameters::scriptSet(const char *param, const char *value) {
     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);
@@ -1137,6 +1171,32 @@ void SimParameters::config_parser_methods(ParseOptions &opts) {
    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);
@@ -3401,6 +3461,19 @@ void SimParameters::check_config(ParseOptions &opts, ConfigList *config, char *&
       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;
@@ -5096,6 +5169,26 @@ if ( openatomOn )
      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";
index c1036c9..1a6b611 100644 (file)
@@ -428,6 +428,20 @@ public:
        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];