Monitor energies and RMSD on the fly to decide when to stop the MD simulation.

[Thomas Evangelidis]

Greetings,

I conduct high through-put MD simulations to equilibrate multiple protein-ligand complexes. However, the required simulation time for equilibration is different for each complex. Is it possible to monitor properties like the protein's and ligand's potential energies, the interaction energy and the ligand RMSD from the previous frame, on the fly and set up some criteria to stop the MD simulation when these properties pass bellow some threshold value?

[Peter Eastman]

There's no simple way to do that, and some of those quantities aren't that well defined. Does the "ligand's potential energy" include its interactions with solvent? And by its presence it alters how the protein interacts with solvent. How do you compute that, and what do you count it as? Also, with PME it's hard to break up the system into components like that. It computes the interaction of each atom with the total electric field at its location, without distinguishing what that field came from.

Anyway, the best approach is probably to do these calculations only occasionally. Run for 1000 steps, compute whatever metrics you want, then decide whether to continue for another 1000 steps.

[Thomas Evangelidis]

Thank you. I don't use PME for the long-range electrostatics but the reaction field method.
At least regarding the RMSD of the ligand atoms wrt previous frames, can I do this calculation using the Python API or do I have to call an external code like ptraj?

[Peter Eastman]

Using an external code to compute the RMSD is probably simplest.

There might be reasonably simple ways to calculate the energies you want, depending on exactly how you define them. For example, you could create a separate System and Context containing only the ligand and nothing else. Whenever you want to compute the ligand's potential energy, just copy over the coordinates and evaluate the energy of that Context.