{"action":"created","issue":{"url":"https://api.github.com/repos/pandegroup/openmm/issues/1278","labels_url":"https://api.github.com/repos/pandegroup/openmm/issues/1278/labels{/name}","comments_url":"https://api.github.com/repos/pandegroup/openmm/issues/1278/comments","events_url":"https://api.github.com/repos/pandegroup/openmm/issues/1278/events","html_url":"https://github.com/pandegroup/openmm/issues/1278","id":119818696,"number":1278,"title":"Proposed API for constant-pH simulations","user":{"login":"jchodera","id":3656088,"avatar_url":"https://avatars.githubusercontent.com/u/3656088?v=3","gravatar_id":"","url":"https://api.github.com/users/jchodera","html_url":"https://github.com/jchodera","followers_url":"https://api.github.com/users/jchodera/followers","following_url":"https://api.github.com/users/jchodera/following{/other_user}","gists_url":"https://api.github.com/users/jchodera/gists{/gist_id}","starred_url":"https://api.github.com/users/jchodera/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/jchodera/subscriptions","organizations_url":"https://api.github.com/users/jchodera/orgs","repos_url":"https://api.github.com/users/jchodera/repos","events_url":"https://api.github.com/users/jchodera/events{/privacy}","received_events_url":"https://api.github.com/users/jchodera/received_events","type":"User","site_admin":false},"labels":[],"state":"open","locked":false,"assignee":null,"milestone":null,"comments":0,"created_at":"2015-12-01T22:05:58Z","updated_at":"2015-12-01T23:50:59Z","closed_at":null,"body":"I'm discussing with @bas-rustenburg (the student heading up the constant-pH work) what setting up a constant-pH simulation should look like to the OpenMM user using the `app` layer.  We'd like your input in refining these ideas.\r\n\r\nCurrently, setting up an implicit solvent simulation with the `app` layer looks like this:\r\n```python\r\npdb = PDBFile('input.pdb')\r\nforcefield = ForceField('amber99sbildn.xml', 'amber99_obc.xml')\r\nmodeller = Modeller(pdb.topology, pdb.positions)\r\nmodeller.addHydrogens(forcefield, pH=7.0)\r\nsystem = forcefield.createSystem(modeller.topology, nonbondedMethod=NoCutoff)\r\nintegrator = LangevinIntegrator(300*kelvin, 5/picoseconds, 2*femtoseconds)\r\nsimulation = Simulation(modeller.topology, system, integrator)\r\nsimulation.minimizeEnergy()\r\nsimulation.step(10000)\r\n```\r\nIdeally, setting up a constant-pH simulation would change this process very little.\r\n\r\nWe need to specify the following important additional pieces of information:\r\n* the fact that we want to run a constant-pH (rather than fixed protonation state) simulation\r\n* the pH\r\n* which residues/components are allowed to sample multiple protonation states, perhaps making the assumption that all allowed residues/components can do so unless otherwise specified when a constant-pH simulation is selected\r\n* which constant-pH capable forcefield (or scheme for computing protomeric populations, for small molecules) is to be used\r\n* optionally, the user should be able to tune some parameters for the constant-pH simulation (update frequencies, probabilities of joint moves, NCMC switching trajectory lengths, etc.), but it should do the right thing (e.g. automatically calibrate reference energies) by default\r\n\r\nHere are a few possibilities for how this could work:\r\n\r\n### Adding simulation \"updaters\" in analogy to \"reporters\"\r\n\r\nWe could easily fit what we have in [`openmm-constph`](https://github.com/choderalab/openmm-constph) if we extend `Simulation` to have `updaters` (in analogy to `reporters`) that are periodically called:\r\n```python\r\npdb = PDBFile('input.pdb')\r\nforcefield = ForceField('amber99sbildn.xml', 'amber99_obc.xml', 'amber99sbildn-constph.xml') # add additional titratable residue definitions\r\nmodeller = Modeller(pdb.topology, pdb.positions)\r\nmodeller.addHydrogens(forcefield, pH=7.0, titratable=True) # replace titratable residues with fully protonated form\r\nsystem = forcefield.createSystem(modeller.topology, nonbondedMethod=NoCutoff)\r\nintegrator = LangevinIntegrator(300*kelvin, 5/picoseconds, 2*femtoseconds)\r\nconstph = TitrationUpdater(system, forcefield, modeller.topology, pH=7.0, temperature=300*kelvin, updateInterval=500, filename='constph.out')\r\nsimulation = Simulation(modeller.topology, system, integrator)\r\nsimulation.updaters.append(constph) # similar to simulation.reporter\r\nsimulation.minimizeEnergy()\r\nconstph.calibrate() # calibrate reference energies\r\nsimulation.step(10000)\r\n```\r\nDescribing a few changes:\r\n*  `amber99sbildn-constph.xml` is an `ffxml` file that uses a few extensions to the `ffxml` schema (#1183) that specifies titratable variants for titratable residues (in addition to some other info described later)\r\n* The flag `titratable` defaults to `False` for standard behavior, but when set to `True`,\r\n```python\r\nmodeller.addHydrogens(forcefield, pH=7.0, titratable=True)\r\n```\r\ncauses the `variants` residue slots for titratable residues to be populated with their fully protonated titratable variants [instead of just `None`](https://github.com/pandegroup/openmm/blob/master/wrappers/python/simtk/openmm/app/modeller.py#L662-L666). Users can prevent certain residues from being titratable by specifying their variants manually via the `variants` argument, as usual.\r\n* The `pHTitrator` is a type of `SimulationUpdater` that can modify the `Context` being simulated and is called with a given update interval.  This class uses `updateParametersInContext` to carry out constant pH simulations and track statistics of the constant-pH simulation.  In principle, this could be a `reporter` since it might also be expected to output data about the protonation state change statistics.\r\n* `ForceFIeld` will be extended to load some additional information about the titratable sites in the `ffxml` extension\r\n* `constph.calibrate()` needs to run some short simulations on reference compounds using the same nonbonded and GB parameters as are used in `system` in order to optimize the reference energies of different titration states\r\n\r\nBenefits:\r\n* Minimal impact to current API\r\n* Easy for users to specify constant pH\r\n\r\nDrawbacks:\r\n* We are relying on `updateParametersInContext` to carry out constant-pH simulations, which may be inefficient; possibly fine for implicit but maybe impractically slow for explicit solvent\r\n* May require API changes if we switch to more efficient models\r\n\r\n### Transform the system to use `Custom*Forces` and a `CustomIntegrator` for constant-pH\r\n\r\nWe could instead transform the `System` object into one that uses context parameters to allow protonation states to be rapidly switched:\r\n```python\r\npdb = PDBFile('input.pdb')\r\nforcefield = ForceField('amber99sbildn.xml', 'amber99_obc.xml', 'amber99sbildn-constph.xml') # add additional titratable residue definitions\r\nmodeller = Modeller(pdb.topology, pdb.positions)\r\nmodeller.addHydrogens(forcefield, pH=7.0, titratable=True) # replace titratable residues with fully protonated form\r\nsystem = forcefield.createSystem(modeller.topology, nonbondedMethod=NoCutoff)\r\nintegrator = LangevinIntegrator(300*kelvin, 5/picoseconds, 2*femtoseconds)\r\n# We might be able to hide some of this complexity inside of `Simulation`\r\nconstph_integrator = TitrationIntegrator(system, forcefield, temperature=300*kelvin, pH=7, ncmc_steps=50)\r\ncompound_integrator = CompoundIntegrator()\r\ncompound_integrator.addIntegrator(integrator)\r\nintegrator_index = compound_integrator.addIntegrator(constph_integrator)\r\nconstph = TitrationUpdater(system, forcefield, modeller.topology, compound_integrator, integrator_index, updateInterval=50, filename='constph.out')\r\nsimulation = Simulation(modeller.topology, system, compound_integrator)\r\nsimulation.updaters.append(constph) # would need to call `constph_integrator` to do constant pH updates and write out statistics\r\nsimulation.minimizeEnergy()\r\nconstph.calibrate() # calibrate constant-pH\r\nsimulation.step(10000)\r\n```\r\n* `ForceField` would create a `System` object with `NonbondedForce` terms for the titratable residues converted to `CustomNonbondedForce` terms with context parameters that select the protonation state\r\n* Much of the complexity of creating the modified system and creating a `CompoundIntegrator` could be hidden inside of `Simulation` if some extra parameters were added to `Simulation()`. \r\n* The `constph` here both *updates* the protonation state *and* prints out statistics, but could be split into an `updater` and a `reporter` if desired\r\n* `constph.calibrate()` needs to run some short simulations on reference compounds using the same nonbonded and GB parameters as are used in `system` in order to optimize the reference energies of different titration states\r\n* Some of this might be simplified by extending `CompoundIntegrator` to be able to define cycles of `nsteps1` with `integrator1` (e.g. Langevin) followed by `nsteps2` with `integrator2` (constant pH titration).\r\n\r\nPros:\r\n* Allows us to more efficiently switch protonation states with context parameters and `CustomIntegrators`\r\n\r\nCons:\r\n* More complicated since it juggles multiple integrators\r\n* Modifies system to use Custom*Force terms instead of original forces; full support for PME may still require use of `updateParametersInContext` to handle reciprocal-space terms\r\n* Some `CustomForce` implementations may need to make duplicates of atoms in residues representing different titration states and append them to `System`, which will then not match the number of atoms in the `Topology`.  To fix this, we would need to create modified trajectory reporters that are aware of this potential difference.\r\n\r\n### Open questions: \r\n* Do we still want to keep support for Amber `cpin` files for users bringing in systems through Amber `inpcrd/prmtop` file readers?\r\n"},"comment":{"url":"https://api.github.com/repos/pandegroup/openmm/issues/comments/161134433","html_url":"https://github.com/pandegroup/openmm/issues/1278#issuecomment-161134433","issue_url":"https://api.github.com/repos/pandegroup/openmm/issues/1278","id":161134433,"user":{"login":"peastman","id":4379786,"avatar_url":"https://avatars.githubusercontent.com/u/4379786?v=3","gravatar_id":"","url":"https://api.github.com/users/peastman","html_url":"https://github.com/peastman","followers_url":"https://api.github.com/users/peastman/followers","following_url":"https://api.github.com/users/peastman/following{/other_user}","gists_url":"https://api.github.com/users/peastman/gists{/gist_id}","starred_url":"https://api.github.com/users/peastman/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/peastman/subscriptions","organizations_url":"https://api.github.com/users/peastman/orgs","repos_url":"https://api.github.com/users/peastman/repos","events_url":"https://api.github.com/users/peastman/events{/privacy}","received_events_url":"https://api.github.com/users/peastman/received_events","type":"User","site_admin":false},"created_at":"2015-12-01T23:50:59Z","updated_at":"2015-12-01T23:50:59Z","body":"```XML\r\nforcefield = ForceField('amber99sbildn.xml', 'amber99_obc.xml', 'amber99sbildn-constph.xml') # add additional titratable residue definitions\r\n```\r\n\r\nWhat residue definitions would be in `amber99sbildn-constph.xml`?  Doesn't `amber99sbildn.xml` already include the fully protonated forms?  If so, those templates will already match.\r\n\r\n```XML\r\nmodeller.addHydrogens(forcefield, pH=7.0, titratable=True) # replace titratable residues with fully protonated form\r\n```\r\n\r\nIs that essentially the same as saying `pH=0`?  It means anything that can be protonated should be, right?\r\n\r\n> We are relying on updateParametersInContext to carry out constant-pH simulations, which may be inefficient; possibly fine for implicit but maybe impractically slow for explicit solvent\r\n\r\nThat's something we should test before we start.  How expensive is it to call `updateParametersInContext()`?  We should measure it and determine whether it's going to be a problem or not.  It's also possible it could be further optimized.\r\n\r\n> More complicated since it juggles multiple integrators\r\n\r\nAren't you going to be doing that either way?  You said before that each Monte Carlo step might require doing an integration.  That seems like a separate question from whether you use global or per-particle parameters to control protonation states."}}