{"action":"created","issue":{"url":"https://api.github.com/repos/choderalab/openmm-constph/issues/19","repository_url":"https://api.github.com/repos/choderalab/openmm-constph","labels_url":"https://api.github.com/repos/choderalab/openmm-constph/issues/19/labels{/name}","comments_url":"https://api.github.com/repos/choderalab/openmm-constph/issues/19/comments","events_url":"https://api.github.com/repos/choderalab/openmm-constph/issues/19/events","html_url":"https://github.com/choderalab/openmm-constph/issues/19","id":137082783,"number":19,"title":"Brainstorm simulation API proposal","user":{"login":"bas-rustenburg","id":3977731,"avatar_url":"https://avatars.githubusercontent.com/u/3977731?v=3","gravatar_id":"","url":"https://api.github.com/users/bas-rustenburg","html_url":"https://github.com/bas-rustenburg","followers_url":"https://api.github.com/users/bas-rustenburg/followers","following_url":"https://api.github.com/users/bas-rustenburg/following{/other_user}","gists_url":"https://api.github.com/users/bas-rustenburg/gists{/gist_id}","starred_url":"https://api.github.com/users/bas-rustenburg/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/bas-rustenburg/subscriptions","organizations_url":"https://api.github.com/users/bas-rustenburg/orgs","repos_url":"https://api.github.com/users/bas-rustenburg/repos","events_url":"https://api.github.com/users/bas-rustenburg/events{/privacy}","received_events_url":"https://api.github.com/users/bas-rustenburg/received_events","type":"User","site_admin":false},"labels":[],"state":"open","locked":false,"assignee":null,"milestone":null,"comments":20,"created_at":"2016-02-28T18:57:37Z","updated_at":"2016-03-06T01:00:19Z","closed_at":null,"body":"So we can discuss it here:  \r\n\r\n```py\r\n\r\nfrom __future__ import print_function\r\nfrom simtk.openmm import app\r\n# TODO  from constph import ConstpHForceField ?\r\nimport simtk.openmm as mm\r\nfrom simtk import unit\r\nfrom sys import stdout\r\n\r\npdb = app.PDBFile('input.pdb') # preprocessed with right residue names to indicate constph?\r\n\r\nforcefield = app.ConstpHForceField('amber99cph.xml', 'tip3p.xml', 'ligandcph.xml') # subclass of forcefield that supports custom format residues\r\n\r\nsystem = forcefield.createSystem(pdb.topology, \r\n nonbondedMethod=app.PME,nonbondedCutoff=1.0*unit.nanometers, constraints=app.HBonds, rigidWater=True, ewaldErrorTolerance=0.0005, cph_indices=None, ph=7.4)\r\n# cph_indices, if None, all that can be matched set up as constant ph, else, list of indices\r\nintegrator = mm.LangevinIntegrator(300*unit.kelvin, 1.0/unit.picoseconds, 2.0*unit.femtoseconds)\r\nintegrator.setConstraintTolerance(0.00001)\r\n\r\n# TODO define compound integrator here too, or leave for simulation?\r\nncmcintegrator =  VelocityVerletIntegrator(ncmc_timestep)\r\nsystem.addForce(mm.MonteCarloBarostat(1*unit.atmospheres, 300*unit.kelvin, 25))\r\n\r\nplatform = mm.Platform.getPlatformByName('CUDA')\r\nproperties = {'CudaPrecision': 'mixed'}\r\n# compound integrator under the hood, hide system update under the hood\r\nsimulation = app.Simulation(pdb.topology, system, {'md': integrator, 'ph': ncmcintegrator}, platform, properties)\r\nsimulation.context.setPositions(pdb.positions)\r\n\r\nprint('Minimizing...')\r\nsimulation.minimizeEnergy()\r\nsimulation.context.setVelocitiesToTemperature(300*unit.kelvin)\r\n\r\nsimulation.reporters.append(app.DCDReporter('trajectory.dcd', 1000)) # modify this to only write out active protons?\r\n# report protonation states\r\nsimulation.reporters.append(app.StateDataReporter(stdout, 1000, step=True,\r\n    potentialEnergy=True, temperature=True, progress=True, remainingTime=True,\r\n    speed=True, totalSteps=1000, protonationStates=True, separator='\\t'))\r\n\r\nprint('Running Production...')\r\nsimulation.calibrate(10000, ph_every=(100,1)) # 10000 cph calibration steps, updating ref energies\r\nsimulation.step(1000000, ph_every=(6000, 1)) # run mc step every 6000 md steps, 1 attempt\r\nprint('Done!')\r\n```\r\n"},"comment":{"url":"https://api.github.com/repos/choderalab/openmm-constph/issues/comments/192775380","html_url":"https://github.com/choderalab/openmm-constph/issues/19#issuecomment-192775380","issue_url":"https://api.github.com/repos/choderalab/openmm-constph/issues/19","id":192775380,"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},"created_at":"2016-03-06T01:00:19Z","updated_at":"2016-03-06T01:00:19Z","body":"Apologies for not responding to this one earlier---I hadn't seen it when I posted my last comment.\r\n\r\n> An important question is whether this should be implemented in Python or C++. If it's in Python, we wouldn't be able to use it in FAH, so if we wanted constant pH there, we would have to implement it separately as part of the core (but that could be a simplified version with fewer options than the more general OpenMM version). \r\n\r\nI think we need to find a way to get the fundamental capabilities into the C++ layer in order for this to be easily supported on FAH.\r\n\r\n> If it's in C++, it's impossible to handle parameter updating in any sort of general way. It would need hardcoded support for every Force subclass it could work with.\r\n\r\nI'm still trying to think of whether there is *some* way to make this general while not requiring extensions to *all* `Force` APIs.\r\n\r\n> If we decide to go with C++, what about implementing it as an Integrator rather than a Force? This would be less brittle and require fewer hacks. It could be similar to a CompoundIntegrator, in that you would create two other Integrators, one for it to use in ordinary integration, and one for it to use when attempting MC steps. Then you would just call step() on it as usual, and periodically it would throw in an MC step.\r\n\r\nThis is a really interesting idea.  Were you thinking that this integrator would be specific to constant pH, or that it would be somehow more general/programmable?  Could it be combined with `CompoundIntegrator`?  Why not just use `CompoundIntegrator` to combine this integrator (which could *just* do the constant-pH steps) with some other integrator?\r\n\r\nThe algorithm would roughly go something like this:\r\n* Use some scheme to select which titratable groups will attempt to change protonation state.  We're exploring this right now, but the probability one or more groups are selected will probably depend on the current protonation state to increase efficiency.\r\n* Select a new protonation state for each titratable group that will attempt to change protonation state.\r\n* If `nsteps = 0`, instantaneously change the parameters, compute the new potential energy, and accept/reject according to some criteria.\r\n* If `nsteps > 0`, store the initial total energy, and then take `nsteps+1` steps changing the parameters by linear interpolation in between each step.  Compute the final total energy and accept/reject according to some criteria.  If accepted, flip the sign of the velocity.\r\n\r\n> But even assuming you wanted to use only CustomNonbondedForce and CustomGBForce, what would you do? There would be possibly dozens of parameters, each controlling the state of one site. And then you would have to write a huge expression for the energy, involving all of those dozens of global parameters. It would be insanely inefficient.\r\n\r\nThere are currently more efficient ways to do this that involve duplicating the titratable residues and using a small set of global parameters to select the titration state of each one.  In the `Custom*Force` terms, there is simply a check to see if the current titration state is active or not, and the interactions are switched on and off according to a continuous parameter.  It's not incredibly elegant, but it doesn't require energy expressions with dozens of parameters.  But I wonder if there would be a way to set (or interpolate) parameters from a table using a global variable instead."}}