Add RPMD features

src/data/npt.py

@@ -312,7 +312,7 @@ def main():
     # Number of threads, multiple timestep integrator, anisotropic box etc.
     threads = TgtOptions.get('md_threads', 1)
     mts = TgtOptions.get('mts_integrator', 0)
-    rpmd_beads = TgtOptions.get('rpmd_beads', 0)
+    rpmd_beads = TgtOptions.get('rpmd_beads', [])
     force_cuda = TgtOptions.get('force_cuda', 0)
     anisotropic = TgtOptions.get('anisotropic_box', 0)
     minimize = TgtOptions.get('minimize_energy', 1)
@@ -366,7 +366,7 @@ def main():
         EngOpts["gas"]["gmx_top"] = os.path.splitext(gas_fnm)[0] + ".top"
         EngOpts["gas"]["gmx_mdp"] = os.path.splitext(gas_fnm)[0] + ".mdp"
         if force_cuda: logger.warn("force_cuda option has no effect on Gromacs engine.")
-        if rpmd_beads > 0: raise RuntimeError("Gromacs cannot handle RPMD.")
+        if len(rpmd_beads) > 0: raise RuntimeError("Gromacs cannot handle RPMD.")
         if mts: logger.warn("Gromacs not configured for multiple timestep integrator.")
         if anisotropic: logger.warn("Gromacs not configured for anisotropic box scaling.")
     elif engname == "tinker":
@@ -375,7 +375,7 @@ def main():
         EngOpts["liquid"]["tinker_key"] = os.path.splitext(liquid_fnm)[0] + ".key"
         EngOpts["gas"]["tinker_key"] = os.path.splitext(gas_fnm)[0] + ".key"
         if force_cuda: logger.warn("force_cuda option has no effect on Tinker engine.")
-        if rpmd_beads > 0: raise RuntimeError("TINKER cannot handle RPMD.")
+        if len(rpmd_beads) > 0: raise RuntimeError("TINKER cannot handle RPMD.")
         if mts: logger.warn("Tinker not configured for multiple timestep integrator.")
     EngOpts["liquid"].update(GenOpts)
     EngOpts["gas"].update(GenOpts)

src/forcefield.py

@@ -223,7 +223,7 @@ def __init__(self, options, verbose=True):
         ## AMOEBA mutual dipole convergence tolerance.
         self.set_option(options, 'amoeba_eps')
         ## Switch for rigid water molecules
-        self.set_option(options, 'rigid_water')
+        self.set_option(options, 'rigid_water', forceprint=True)
         ## Bypass the transformation and use physical parameters directly
         self.set_option(options, 'use_pvals')
 
@@ -1115,7 +1115,8 @@ def build_qtrans2(tq, qid, qmap):
     def list_map(self):
         """ Create the plist, which is like a reversed version of the parameter map.  More convenient for printing. """
         if len(self.map) == 0:
-            warn_press_key('The parameter map has no elements (Okay if we are not actually tuning any parameters.)')
+            #warn_press_key('The parameter map has no elements (Okay if we are not actually tuning any parameters.)')
+            logger.warning('The parameter map has no elements (Okay if we are not actually tuning any parameters.)\n')
         else:
             self.plist = [[] for j in range(max([self.map[i] for i in self.map])+1)]
             for i in self.map:

src/openmmio.py

@@ -209,6 +209,16 @@ def CopyNonbondedParameters(src, dest):
     for i in range(src.getNumExceptions()):
         dest.setExceptionParameters(i,*src.getExceptionParameters(i))
 
+def CopyCustomBondParameters(src, dest):
+    dest.setEnergyFunction(src.getEnergyFunction())
+    for i in range(src.getNumGlobalParameters()):
+        dest.setGlobalParameterName(i,src.GlobalParameterName(i))
+    for i in range(src.getNumGlobalParameters()):
+        dest.setGlobalParameterDefaultValue(i,*src.getGlobalParameterDefaultValue(i))
+    for i in range(src.getNumPerBondParameters()):
+        dest.setPerBondParameterName(i,src.getPerBondParameterName(i))
+    for i in range(src.getNumBonds()):
+        dest.setBondParameters(i,*src.getBondParameters(i))
 def do_nothing(src, dest):
     return
 
@@ -226,6 +236,7 @@ def CopySystemParameters(src,dest):
                'HarmonicAngleForce':CopyHarmonicAngleParameters,
                'PeriodicTorsionForce':CopyPeriodicTorsionParameters,
                'NonbondedForce':CopyNonbondedParameters,
+               'CustomBondForce':CopyCustomBondParameters,
                'CMMotionRemover':do_nothing}
     for i in range(src.getNumForces()):
         nm = src.getForce(i).__class__.__name__
@@ -566,35 +577,68 @@ def prepare(self, pbc=False, mmopts={}, **kwargs):
         self.AtomLists['ResidueNumber'] = [a.residue.index for a in Atoms]
         self.AtomMask = [a == 'A' for a in self.AtomLists['ParticleType']]
 
-    def create_simulation(self, timestep=1.0, faststep=0.25, temperature=None, pressure=None, anisotropic=False, mts=False, collision=1.0, nbarostat=25, rpmd_beads=0, **kwargs):
+    def create_simulation(self, timestep=1.0, faststep=0.25, temperature=None, pressure=None, anisotropic=False, mts=False, collision=1.0, nbarostat=25, rpmd_beads=[], **kwargs):
 
         """
         Create simulation object.  Note that this also takes in some
         options pertinent to system setup, including the type of MD
         integrator and type of pressure control.
         """
-
         # Divisor for the temperature (RPMD sets it to nonzero.)
         self.tdiv = 1
 
         ## Determine the integrator.
         if temperature:
             ## If temperature control is turned on, then run Langevin dynamics.
             if mts:
-                if rpmd_beads > 0:
+                if len(rpmd_beads) > 0:
                     raise RuntimeError("No multiple timestep integrator without temperature control.")
                 integrator = MTSVVVRIntegrator(temperature*kelvin, collision/picosecond,
                                                timestep*femtosecond, self.system, ninnersteps=int(timestep/faststep))
             else:
-                if rpmd_beads > 0:
-                    logger.info("Creating RPMD integrator with %i beads.\n" % rpmd_beads)
-                    self.tdiv = rpmd_beads
-                    integrator = RPMDIntegrator(rpmd_beads, temperature*kelvin, collision/picosecond, timestep*femtosecond)
+                if len(rpmd_beads) > 0:
+                    self.tdiv = int(rpmd_beads[0])
+                    if len(rpmd_beads) == 1:
+                        logger.info("Creating RPMD integrator with %i beads.\n" % int(rpmd_beads[0]))
+                        integrator = RPMDIntegrator(int(rpmd_beads[0]), temperature*kelvin, collision/picosecond, timestep*femtosecond)
+                    elif len(rpmd_beads) == 2:
+                        contract = False
+                        for frc in self.system.getForces():
+                            if any([isinstance(frc, fc) for fc in [NonbondedForce, AmoebaMultipoleForce, AmoebaVdwForce, CustomNonbondedForce]]):
+                                contract = True
+                                frc.setForceGroup(1)
+                        if contract:
+                            logger.info("Creating RPMD integrator with %i beads (NB forces contracted to %i).\n" % (int(rpmd_beads[0]), int(rpmd_beads[1])))
+                            integrator = RPMDIntegrator(int(rpmd_beads[0]), temperature*kelvin, collision/picosecond, timestep*femtosecond, {1:int(rpmd_beads[1])})
+                        else:
+                            logger.info("Creating RPMD integrator with %i beads (no NB forces to contract).\n" % (int(rpmd_beads[0])))
+                            integrator = RPMDIntegrator(int(rpmd_beads[0]), temperature*kelvin, collision/picosecond, timestep*femtosecond)
+                    elif len(rpmd_beads) == 3:
+                        contract = False
+                        contract_recip = False
+                        for frc in self.system.getForces():
+                            if any([isinstance(frc, fc) for fc in [NonbondedForce, AmoebaMultipoleForce, AmoebaVdwForce, CustomNonbondedForce]]):
+                                contract = True
+                                frc.setForceGroup(1)
+                                if isinstance(frc, NonbondedForce):
+                                    contract_recip = True
+                                    frc.setReciprocalSpaceForceGroup(2)
+                        if contract_recip:
+                            logger.info("Creating RPMD integrator with %i beads (NB/Recip forces contracted to %i/%i).\n" % (int(rpmd_beads[0]), int(rpmd_beads[1]), int(rpmd_beads[2])))
+                            integrator = RPMDIntegrator(int(rpmd_beads[0]), temperature*kelvin, collision/picosecond, timestep*femtosecond, {1:int(rpmd_beads[1]), 2:int(rpmd_beads[2])})
+                        elif contract:
+                            logger.info("Creating RPMD integrator with %i beads (NB forces contracted to %i, no Recip).\n" % (int(rpmd_beads[0]), int(rpmd_beads[1])))
+                            integrator = RPMDIntegrator(int(rpmd_beads[0]), temperature*kelvin, collision/picosecond, timestep*femtosecond, {1:int(rpmd_beads[1])})
+                        else:
+                            logger.info("Creating RPMD integrator with %i beads (no NB forces to contract).\n" % (int(rpmd_beads[0])))
+                            integrator = RPMDIntegrator(int(rpmd_beads[0]), temperature*kelvin, collision/picosecond, timestep*femtosecond)
+                    else:
+                        raise RuntimeError("Please provide a list of length 1, 2, or 3 to rpmd_beads")
                 else:
                     integrator = LangevinIntegrator(temperature*kelvin, collision/picosecond, timestep*femtosecond)
         else:
             ## If no temperature control, default to the Verlet integrator.
-            if rpmd_beads > 0:
+            if len(rpmd_beads) > 0:
                 raise RuntimeError("No RPMD integrator without temperature control.")
             if mts: warn_once("No multiple timestep integrator without temperature control.")
             integrator = VerletIntegrator(timestep*femtoseconds)
@@ -610,13 +654,13 @@ def create_simulation(self, timestep=1.0, faststep=0.25, temperature=None, press
         elif pressure != None: warn_once("Pressure is ignored because pbc is set to False.")
 
         ## Set up for energy component analysis.
-        if not mts:
+        if not mts and len(rpmd_beads) == 0:
             for i, j in enumerate(self.system.getForces()):
                 j.setForceGroup(i)
 
         ## If virtual particles are used with AMOEBA...
         SetAmoebaVirtualExclusions(self.system)
-        
+
         ## Finally create the simulation object.
         self.simulation = Simulation(self.mod.topology, self.system, integrator, self.platform)
 
@@ -951,8 +995,24 @@ def molecular_dynamics(self, nsteps, timestep, temperature=None, pressure=None,
             if iteration >= 0:
                 self.simulation.step(nsave)
             state = self.simulation.context.getState(getEnergy=True,getPositions=True,getVelocities=False,getForces=False)
+#####
+#            if self.tdiv>=2:
+#                hbar=0.0635078*kilojoule*picosecond/mole
+#                pimdstate=["?"]*self.tdiv
+#                kinetic=0.0*kilojoule/mole
+#                potential=0.0*kilojoule/mole
+#                for i in range(self.tdiv):
+#                    pimdstate[i]=integrator.getstate(i,getEnergy=True,getPositions=True,group=-1)
+#                    potential=potential+pimdstate[i].getPotentialEnergy()/self.tdiv
+#                    kinetic=kinetic+pimdstate[i].getKineticEnergy()/self.tdiv
+#                for i in range(self.tdiv):
+#                    ii=(i+1)%self.tdiv
+#                    deltabead=np.array(pimdstate[i].getPositions())-np.array(pimdstate[ii].getPositions())
+#                    kinetic=kinetic-np.sum(((deltabead*deltabead).sum(axis=1))*np.array([getParticleMass(j) for j in range (system.getNumParticles())]))*(kB**2*integrator.getTemperature()**2*self.tdiv)/(2.0*hbar**2)
+#            else:
             kinetic = state.getKineticEnergy()/self.tdiv
             potential = state.getPotentialEnergy()
+#####Energy Estimator
             if self.pbc:
                 box_vectors = state.getPeriodicBoxVectors()
                 volume = self.compute_volume(box_vectors)

src/parser.py

@@ -77,13 +77,13 @@
                  },
     'lists'   : {"forcefield"     : ([],  200, 'The names of force fields, corresponding to directory forcefields/file_name.(itp,xml,prm,frcmod,mol2)', 'All (important)'),
                  "scanindex_num"  : ([], -100, 'Numerical index of the parameter to scan over', 'Job types scan_mvals and scan_pvals'),
-                 "scanindex_name" : ([], -100, 'Parameter name to scan over (should convert to a numerical index)', 'Job types scan_mvals and scan_pvals')
+                 "scanindex_name" : ([], -100, 'Parameter name to scan over (should convert to a numerical index)', 'Job types scan_mvals and scan_pvals'),
+                 "rpmd_beads"     : ([], -160, 'Number of beads / NB contracted beads / Recip contracted beads in ring polymer MD', 'Condensed phase property targets (advanced usage)', 'liquid_openmm'),
                  },
     'ints'    : {"maxstep"      : (100, 50, 'Maximum number of steps in an optimization', 'Main Optimizer'),
                  "objective_history"  : (2, 20, 'Number of good optimization steps to average over when checking the objective convergence criterion', 'Main Optimizer (jobtype "newton")'),
                  "wq_port"   : (0, 0, 'The port number to use for Work Queue', 'Targets that use Work Queue (advanced usage)'),
                  "criteria"   : (1, 160, 'The number of convergence criteria that must be met for main optimizer to converge', 'Main Optimizer'),
-                 "rpmd_beads"       : (0, -160, 'Number of beads in ring polymer MD (zero to disable)', 'Condensed phase property targets (advanced usage)', 'liquid_openmm'),
                  },
     'bools'   : {"backup"           : (1,  10,  'Write temp directories to backup before wiping them'),
                  "writechk_step"    : (1, -50,  'Write the checkpoint file at every optimization step'),