mdtraj
bond angle computation results in 'nan'
I have a molecule 'mol.gro' ''' Computing the bond angle results in 'nan' 3 1MOL C 1 0.475 -1.452 -0.467 1MOL C 2 0.459 -1.348 -0.418 1MOL C 3 0.443 -1.244 -0.369 9.0 9.0 9.0 0.0 0.0 0.0 0.0 0.0 0.0 ''' Do this to compute the bond angle ''' import mdtraj as md a=md.load('mol.gro') b=md.compute_angles(a, [[0,1,2]]) ''' That result in b=[[nan]]. I think this is because the bond angle is close to 180 and numerical errors make something like 'arccos -1.00000000001'
This is a pretty common issue in this space, I'm not sure if it's properly been dealt with here
In general you might not want to be using low-precision formats like GRO files for trajectory analysis, but I don't know the context here.
https://github.com/mdtraj/mdtraj/issues/1499
This has definitely come up often before. In some cases I've found it maybe useful to flip around the bond angles because then you might get zero instead (although this would require extra bookeeping).
Aside from using higher precision formats for storing coordinates, I'm not entirely too sure how to fix this so I'm open to suggestions!
Unfortunately there's some sort of issue with your mol.gro file that you pasted in as I'm having trouble loading it so I can't replicate this
I can reproduce this error with a version of the .gro file with the formatting fixed (see https://manual.gromacs.org/archive/5.0.3/online/gro.html). I agree with @Enzymesf that it looks like the result of numerical instabilities as the bond angle gets close towards 180 degrees. Adding code to angle.py and anglekernels.h to clip the calculated cosine into the range (-1.0, 1.0) fixes the problem, at least for this example. I can provide a PR if wanted.
it looks like the result of numerical instabilities as the bond angle gets close towards 180 degrees. Adding code to angle.py and anglekernels.h to clip the calculated cosine into the range (-1.0, 1.0) fixes the problem
Ahhhhhh that makes sense!
Adding code to angle.py and anglekernels.h to clip the calculated cosine into the range (-1.0, 1.0) fixes the problem, at least for this example. I can provide a PR if wanted.
I would welcome a PR to address this long-standing issue! We'll probably also need to include some tests simultaneously to ensure angle computation still works as expected elsewhere within a reasonable degree of precision.