Changing molecule force field during simulation?

[George Pantelopulos]

Hi all,

Is it at all possible to, for a particular molecule, to scale its force field parameters between two sets of force field parameters? This is like doing an alchemical transformation of one molecule to another molecule, but without growing in any particles or changing the coordinates of the molecule.

Of course someone could rebuild the system with a new force field but this would obviously be very slow, and I would like to be able to do such a transformation very frequently.

Thank you for any help,
George

[Peter Eastman]

Yes, there are two ways you could do this. The first is to modify the parameters directly. For example, call setParticleParameters() on a NonbondedForce to modify the parameters for one or more particles. After doing that, call updateParametersInContext() so the changes will be propagated to the existing context.

The other approach is to replace the standard forces with custom ones, for example a CustomNonbondedForce instead of NonbondedForce. Make that directly depend on a global parameter which you then vary between 0 and 1. This can be a little more efficient in some situations, but for your case, the first approach is probably the easier one. And unless you're modifying the parameters almost every time step, the difference in efficiency will probably be negligible.

If you're doing alchemical simulations, you might also want to look at Yank (http://getyank.org/latest), which is based on OpenMM.

Peter

[George Pantelopulos]

Hi Peter,

Thank you again for the helpful information. I am not necessarily trying to do FEP, but the idea is very similar, of course. I am trying to do a sort of lambda dynamics in which the force fields of a selection of molecules may be mutated to another force field via linear combinations of their parameters with a penalty function which restrains the average value of the lambdas for the selected molecules.

I do understand much better how to modify every particle parameter at every time step and update the context parameters, but that would require re-writing the script every time a system containing new force fields to perform this calculation on, right? It would also, as you say, cause major slowdowns.

Might it be possible to write a CustomNonbondedForce which not only scales the forces per particle by a lambda which is unique to that particle's residue but also employs a non-physical penalty function which restrains the average value of the lambdas for all selected molecules? And would it be possible to linearlly combine the parameters that would be on a particle between two sets of force fields in such a force?

[Peter Eastman]

but that would require re-writing the script every time a system containing new force fields to perform this calculation on, right?

Why do you say that? It should be possible to implement it in a generic way. Obviously you need to specify in some way which parameters are changing, and what values they're changing between. But given that information in some reasonable format, everything else should be generic.

It would also, as you say, cause major slowdowns.

Actually I said it wouldn't cause major slowdowns. Especially if you're only updating the parameters periodically, say every 1000 time steps.

Might it be possible to write a CustomNonbondedForce which not only scales the forces per particle by a lambda which is unique to that particle's residue but also employs a non-physical penalty function which restrains the average value of the lambdas for all selected molecules?

If you want to do lambda dynamics, that's a bit different. It requires not only having the lambdas scale the forces, but also that their values get updated by the integrator, which in turn requires computing the forces acting on them (that is, derivatives of the energy with respect to the lambdas). Is that what you need? Or is it sufficient to just modify the lambdas directly?

[George Pantelopulos]

Hi Peter,

Why do you say that? It should be possible to implement it in a generic way. Obviously you need to specify in some way which parameters are changing, and what values they're changing between. But given that information in some reasonable format, everything else should be generic.

Oh, really? Is it possible to do something like load both force fields and scale the subsets of forces that belong to each of them based on their index?

Actually I said it wouldn't cause major slowdowns. Especially if you're only updating the parameters periodically, say every 1000 time steps.

Nice! Well, a MC idea of changing parameters would be pretty frequent, like in the tens of steps, I think would work well for the problem I'm trying to solve. I remember re-setting the context to get new masses wasn't too painful with that older problem you helped with! I should post that somewhere -- it was based on your ST code, actually!

If you want to do lambda dynamics, that's a bit different. It requires not only having the lambdas scale the forces, but also that their values get updated by the integrator, which in turn requires computing the forces acting on them (that is, derivatives of the energy with respect to the lambdas). Is that what you need? Or is it sufficient to just modify the lambdas directly?

Right, it is necessary that the lambda values gets updated by the integrator, and that the lambda values scale the forces between sets of forces which would be input by two different force fields, on each molecule.

[Peter Eastman]

Oh, really? Is it possible to do something like load both force fields and scale the subsets of forces that belong to each of them based on their index?

You aren't interpolating between force fields; you're interpolating between values for specific parameters. You can determine those parameters however you want. For example, you might construct systems from two force fields, then loop over atoms and see which ones have different parameters. Or you might do it in some other way. But the only job of a force field is to construct a system. Once it has done that, it's no longer involved in the simulation. You can then modify the system however you want.

I remember re-setting the context to get new masses wasn't too painful with that older problem you helped with!

Updating parameters is much less expensive than reinitializing the context. Still not free of course, but not terrible.

[George Pantelopulos]

Hi Peter,

I believe I have a pretty good idea of how to try this now -- thank you for the help!