5
$\begingroup$

I am trying to get the potential energy of a system at each frame of a molecular dynamics trajectory after running the simulation. My current approach is as follows:

  1. Use MDAnalysis to iterate through the trajectory frames and write out the atom coordinates to separate topology files
  2. Iterate through each of the files and create a simulation using OpenMM and then use their getPotentialEnergy() method to extract the energy from the simulation state

The issue is that when I test this using the same forcefield that was used for the simulation, I get potential energy values that differ on the order of 10% from the potential energy reported during the simulation. I have also noticed that the first topology file created by step one has atom coordinates that differ significantly from the input topology file, but I do not know why.

What could the issues with my implementation be, or is there an alternate way to accomplish this accurately?

Improve this question
$\endgroup$
6
  • $\begingroup$ Well if the atomic coordinates "differ significantly" it shouldn't be too surprising that the energy is different. I suspect debug step 1) and step 2) will solve itself. As to why step 1) is wrong, well we can't tell with what is provided. The few initial thoughts: Are there periodic boundary conditions? If so are they being implemented correctly? Have you got the lattice vectors correct? If it is a constant pressure ensemble are the changing lattice vectors being updated correctly? $\endgroup$
    – Ian Bush
    Feb 16 at 22:50
  • $\begingroup$ This seems a related question modulo DFT/MD engine: mattermodeling.stackexchange.com/q/10102/4515 ? Try a tighter tolerance on both the MD and the energy calculation? $\endgroup$
    – Andrey Poletayev
    Feb 17 at 1:55
  • $\begingroup$ Are you using NAMD for the simulations? If yes, the issues are fairly easy to solve $\endgroup$
    – Hemanth Haridas
    Feb 17 at 3:52
  • 2
    $\begingroup$ There's not enough information here -- it would be tremendously helpful if you could tell us how you ran the simulation, what system, what force field, etc. Trying to replicate an energy calculation in particular is tricky because there are so many parameters (cutoff distance and handling, long range electrostatics, handling of 1-2 and 1-3 and 1-4 neighbours, etc). $\endgroup$
    – Shern Ren Tee
    Feb 17 at 8:26
  • 1
    $\begingroup$ I think that you are double-working. If you already have a trajectory file, this means that you already run a MD and the energy of each frame is known. About mentioned before, obtaining different energy from two runs is normal (each run will start with a different set of velocities, so, getting a different final state). $\endgroup$
    – Camps ♦
    Feb 17 at 11:43

0

Reset to default

You must log in to answer this question.