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'),