Add RPMD features

[frankmenglin]

Interface to RPMD contraction in OpenMM and kinetic energy estimator in openmmio.py .

[frankmenglin]

The last one "RPMD input should be list of integer not number" should be "RPMD input should be list of integer not string"

For now I just replace every rpmd_beads[#] by int(rpmd_beads[#]) where # is number from 0~2 in openmmio.py but I'm sure there is a better way doing it.

[leeping]

Hey :) It looks like the tests are passing. Is it ready to merge?

[leeping]

I looked over the code and it looks good.

[leeping]

Well, except you haven't implemented the estimator yet. :)

[frankmenglin]

I did not implement estimator as an extra function, just when RPMD integrator is called, energy is calculated differently(KE estimate+PE).

[leeping]

Are these comments intended to be the kinetic energy calculation?

[leeping]

From looking at the code, it appears ForceBalance isn't doing anything different to calculate the kinetic energy in RPMD. There is a commented-out section which appears to contain some code for the kinetic energy calculation - is that supposed to be there?

Do you mean that the KE estimate is implemented directly into OpenMM?

[frankmenglin]

Ouch! I forget to put it back. Yes those are supposed to be energy estimator.

Originally I comment it out to debug, then I find out rpmd bug is related to input format, but I forget to put it back

[frankmenglin]

KE estimate is not implemented in OpenMM yet.

[leeping]

Okay, I understand. Go ahead and fix the energy estimator. This is probably a good thing because we need to differentiate this quantity w/r.t. the force field parameters.

Note that you will need to pull the code from my repository before you push again, and there may be a merge conflict because we were both working on openmmio.py. If a merge conflict appears, you should edit the file manually so that the merge is resolved, commit again, and then push.

[leeping]

By the way, you may initialize pimdstate to an empty list [] and use the append() function to add elements to the list.

[frankmenglin]

A quick question, does your update involve some redefinition of objects like state,system,integrator,etc?

For some reasons I'll get error if using system.getNumParticles() to get # of particles(should use self.system.getNumParticles()).

Another issue is when I tried to call integrator.getState(...some arguments...) it'll create errors, but that used to work.

[frankmenglin]

Hm, that part actually does not do anything extra now. My original intention was to get a new way of calculating Multipole when RPMDintegrator is called, but then I find out I cannot avoid getting state information from context instead of RPMDintegrator when calculating the AMOEBA part.

Let me remove that part for now.

A slightly more general question is, I find out quite a few functions(energy_components,get_multipoles,UpdateSimulationParameter,evaluate_one_, etc) in openmmio.py are using context to get state information, and in RPMD simulation that will only get state information from a single copy of system. I'm not fully understand how each function works and how this should be adjusted in RPMD simulation.

[leeping]

Hi Frank,

UpdateSimulationParameters() is important because it updates the Simulation object with new force field parameters. In the case of changing Force object parameters, it is done using the method Force.updateParametersInContext(). In the case of changing virtual site parameters, the Simulation object is destroyed and created again. However, I don't think this function actually calls the context.getState() function.

Here are the places context.getState() is called:

evaluate_one_ is how we get the potential energy in the postprocessing step, so it is essential that we account for RPMD here. We might also need the quantum KE in the postprocessing step as well, since that is another way for the parameters to affect the Boltzmann factor (am I right?)

get_multipoles is used to get the simulation dipole moments, so it would be good to account for RPMD here.

energy_components is a way to break the total energy down into components, but it's not essential because we don't use it in the parameterization, so you can leave it alone.

Finally, optimize and energy_rmsd are not important because they are geometry optimizations and not MD simulations. You can leave them alone.

Since there are potentially multiple places where getState() needs to be modified for RPMD, I think we should not be adding if statements in several places - rather we should call another function that has the branch. For example:

def evaluate_potential(Sim):
    # This is how you tell whether the integrator is the RPMDIntegrator - better than using class name
    if isinstance(Sim.integrator, RPMDIntegrator):
        # Evaluate RPMD energy properly here.
        return rpmd_energy
    else:
        return Sim.context.getState(getEnergy=True).getPotentialEnergy()

I believe we can have this as a function rather than a method in Engine, because it shouldn't depend on anything else in our Engine object.

The other unresolved issue is how to set the virtual site positions in the RPMD images when we update the virtual site parameters. @peastman: Do you have advice for us on how to do this?

Thanks,

• Lee-Ping

[peastman]

Is it even valid to use virtual sites with RPMD? I know that it isn't valid to use constraints. Tom has said a lot of the equations would need to be reworked, and last I heard, he hadn't done that. I'd expect virtual sites to have similar issues.

[leeping]

Hi Peter,

Thanks. I am pretty sure that q-TIP4P/f (the model that Tom most often uses) does have a virtual site.

I think it's different from constraints, which are physically invalid in RPMD because they correspond to freezing out the same vibrational degrees of freedom that RPMD is intended to model. As long as the virtual site position is correctly computed from the current coordinates of the flexible molecule, it should make physical sense.

[tmarkland]

Hi Peter,

RPMD for models with virtual sites has worked since OpenMM 5.1 (you modified it then to fix the virtual site behavior). RPMD with constraints (e.g. bonds or angles) can be done but it is:

(a.) Not currently implemented in OpenMM.
(b.) Not really advisable since many of the interesting quantum effects come from quantizing the vibrations.

In practice the modification that is needed is to apply the constraints to each copy of the system independently but one has to be careful about integrating this with the free ring polymer evolution which does not preserve the constraints. That said it seems for this problem all that is needed is virtual sites which are already included and working with RPMD in OpenMM.

[peastman]

Ok, that's fine then. Anyway, if you want to change the virtual site parameters, you need to reinitialize or recreate your context. Currently there's no other way to do it.

[leeping]

Hi Peter,

That is how I change the virtual site parameters. My question is how we set the positions in the RPMD images in a trajectory loop over stored coordinates.

In the classical case, after I create a new Context and set the positions manually (from my stored positions), I use Context.computeVirtualSites() to create the virtual site positions for the currently set positions. I don't know how to do this for RPMD since the image positions are stored in the integrator.

Maybe I can do the following loop over stored image positions:

1. Call Context.setPositions() for the image.
2. Call Context.computeVirtualSites()
3. Transfer these positions to the appropriate image in the integrator

Thanks,

• Lee-Ping

[leeping]

Frank: I just talked to Peter and this is a good way to set the virtual site positions in the RPMD images. Basically you need to use self.simulation.context as an "processor" - you set the positions in the context, run context.computeVirtualSites(), then get the positions back out and set them in the integrator.

This seems complicated because normally the virtual site positions are set when you do any MD steps - however, in the postprocessing stage we don't do any MD steps ...

[frankmenglin]

OK, I'll do that.

Btw in the new build, I update the energy estimator. I think I'll open another issue to explain the detail.

[frankmenglin]

While I'm still updating the code, just curious, does it matter if the test in Jenkins shows red light before test finish but end up in success?

[leeping]

Hi Frank,

That doesn't matter. The reason for the red light is because the "Thermo2" pull request breaks the tests, and it is being built on the same "job" as this current pull request.

• Lee-Ping

[frankmenglin]

For some reasons "change nothing" on something fail actually success!