CHARMM Development Project Forums User Discussion & Questions Molecular Dynamics XYZ file output bug

I am using CHARMM to run simualtions of a pure water box with ~2000 TIP4 water. I found the force outputed in XYZ file is not correct. For example, here are the last several lines of the file:

OH -0.117556615442641E+02 -0.460737086741699E+01 -0.504966436489195E+01 -0.896962810587684E-01 -0.173079977399105E+00 -0.481165302110584E-01 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00
OM -0.118031743479324E+02 -0.451734875393761E+01 -0.502282548502928E+01 -0.871213924496163E-01 -0.144800493693594E+00 -0.138447897817385E+00 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00
H1 -0.113836050592983E+02 -0.372543913170627E+01 -0.505115624108677E+01 0.129661680997219E+00 -0.266193695140701E+00 -0.320327127050701E+00 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00
H2 -0.126565853167218E+02 -0.448726215730393E+01 -0.474942763475455E+01 -0.280393034970152E+00 0.234814024439068E+00 -0.781388701678074E+00 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00
OH -0.920377932247452E+01 -0.292489385283675E+01 0.171612180879107E+02 0.261133883253496E+00 0.641836712759974E-02 -0.224156577754866E+00 0.106645767748745E+01 -0.109439792523077E+01 -0.262052863952775E+01
OM -0.914211865733862E+01 -0.295232157244103E+01 0.172420107575801E+02 0.194001142856658E+00 -0.626893644332873E-01 -0.196464073114893E+00 -0.302238294894556E+01 0.950654595829137E+01 0.117020104788113E+02
H1 -0.829465994835501E+01 -0.266934800539263E+01 0.173175194750189E+02 0.280628971884755E+00 -0.510378296330748E+00 0.505309228793788E+00 -0.128182206452567E+01 -0.471861955630476E+01 -0.299333692575178E+01
H2 -0.942655088720442E+01 -0.348573895752148E+01 0.179042237666623E+02 -0.505618946457559E+00 -0.246026430692034E+00 -0.645375792911392E+00 0.240800395831427E+01 -0.208514960990031E+01 -0.348649064796158E+01


The forces (last three columns) on the second last water are all exact zeros, which is impossible. Actually, the force on them are always zeros throughout the XYZ file. Only the first ~60 water molecules seem to have correct force (I calculated based on leap-frog algorithm, as well as those OMs have zero force).

This is verly likely a bug of outputing force in parallel CHARMM MD:
1) I used 32 cores, and ~2000/32 ~= ~60 (number of molecules with correct forces).
2) When I used one core, everything is correct.


Some other points:

1) I found the direction of the force (in XYZ file) is reversed, i.e., it is actually the force acting on the system by the atom, rather than the force acting on the atom. This is even the case in COOR FORCE. Is this the default setup in CHARMM?
2) The coordinates of the XYZ file is shifted by one step, 1st step in XYZ stores coordinate of the 2nd step..., etc.

The 'dynamc' info file mentions the CHARMM convention-- "The force on a particle is equal to the negative gradient of the potential energy of the particle."

The remaining issues appear to be bug in the placement of the write statement for the XYZ file within the dynamics control loop done for each integration step. To me, it also suggest the MXYZ output feature is older code that has not been used or rigorously tested for years. Depending on your goals, post-processing stored dynamics coordinates to get the forces might be a better option.

Thanks, Rick!

I fount that the results from COOR FORCE are a little different from XYZ's, but the resultant force and troque on one molecule are very close.