BAT method Cartesian modifies input data

[JIMonroe]

Expected behavior

When calling the Cartesian method of a BAT class instance to convert BAT coordinates back to Cartesian, the input BAT data should not be modified.

Actual behavior

"Non-primary" torsions in the original input BAT data are modified due to issues with scope and pass-by-reference. The very easy fix is just to call copy.deepcopy() when defining the torsions from the input data here:

mdanalysis/package/MDAnalysis/analysis/bat.py

Line 502 in 4a6aee4

torsions = bat_frame[2 * n_torsions + 9:]

Code to reproduce the behavior

import copy
import numpy as np

import MDAnalysis as mda
from MDAnalysis.analysis import bat

#Just load in a simple test molecule - you can use any molecule as long as it has non-primary torsions
prmtopFile = 'alanine-dipeptide.prmtop'
strucFile = 'alanine-dipeptide.pdb'
uni = mda.Universe(prmtopFile, strucFile)

#Create BAT analysis and convert input coordinates to BAT, noting that we deep copy the data
bat_analysis = bat.BAT(uni.select_atoms('all'))
bat_analysis.run()
bat_coords = copy.deepcopy(bat_analysis.bat[0])
print('Original BAT coordinates:\n', bat_coords)

#Now run test to convert back to Cartesian
print('Original Cartesian coordinates:\n', uni.trajectory[0].positions)
print('First transformation back to Cartesian:\n', bat_analysis.Cartesian(bat_coords))

#It works the first time, but let's check bat_coords and try again...
print('Modified BAT coordinates after one transformation:\n', bat_coords)
print('Transforming again is now incorrect (obviously):\n', bat_analysis.Cartesian(bat_coords))
print('Non-primary torsions keep having primaries added to them:\n', bat_coords)

Dialanine.zip

Current version of MDAnalysis

2.0.0-dev0 (looking at https://github.com/MDAnalysis/mdanalysis/blob/develop/package/MDAnalysis/analysis/bat.py it looks like the current code will still have the same issue)
python 3.7
BigSur, Centos7

[JIMonroe]

By the way, I've found the BAT analysis class very helpful! I just noticed this bug recently when it created some issues for me, but it's a minor problem in a much-appreciated tool.

[AnirG]

@JIMonroe Hey there! The below code snippet doesn't print the same which I believe should be the case.

print('Original Cartesian coordinates:\n', uni.trajectory[0].positions)
print('First transformation back to Cartesian:\n', bat_analysis.Cartesian(bat_coords))

output.txt
Could you please verify? Or am I missing something?

[JIMonroe]

@AnirG Hey, thanks for checking this out! For me, downloading the files and then running the script I have outputs the same coordinates for the first transformation, within 1e-06 or so, as the original (see below). I just created a fresh conda environment and install of MDAnalysis and I get the same thing (though you do have to use bat_analysis.results.bat[0] to get the BAT coords in the latest MDAnalysis version). Based on that, I'm not sure what the issue is, though I don't think it has to do with the version.

Output:

Original Cartesian coordinates:
 [[ 2.     1.    -0.   ]
 [ 2.     2.09   0.   ]
 [ 1.486  2.454  0.89 ]
 [ 1.486  2.454 -0.89 ]
 [ 3.427  2.641 -0.   ]
 [ 4.391  1.877 -0.   ]
 [ 3.555  3.97  -0.   ]
 [ 2.733  4.556 -0.   ]
 [ 4.853  4.614 -0.   ]
 [ 5.408  4.316  0.89 ]
 [ 5.661  4.221 -1.232]
 [ 5.123  4.521 -2.131]
 [ 6.63   4.719 -1.206]
 [ 5.809  3.141 -1.241]
 [ 4.713  6.129  0.   ]
 [ 3.601  6.653  0.   ]
 [ 5.846  6.835  0.   ]
 [ 6.737  6.359 -0.   ]
 [ 5.846  8.284  0.   ]
 [ 4.819  8.648  0.   ]
 [ 6.36   8.648  0.89 ]
 [ 6.36   8.648 -0.89 ]]

First transformation back to Cartesian:
 [[ 1.99999980e+00  1.00000011e+00 -2.05560270e-07]
 [ 1.99999984e+00  2.09000003e+00 -1.57430843e-07]
 [ 1.48599982e+00  2.45400009e+00  8.89999858e-01]
 [ 1.48599977e+00  2.45400016e+00 -8.90000113e-01]
 [ 3.42699990e+00  2.64100010e+00 -1.72364129e-07]
 [ 4.39099963e+00  1.87700000e+00 -2.32622690e-07]
 [ 3.55499997e+00  3.97000009e+00 -1.17203382e-07]
 [ 2.73299996e+00  4.55600032e+00 -6.87115096e-08]
 [ 4.85300009e+00  4.61399986e+00 -1.24480879e-07]
 [ 5.40799993e+00  4.31599994e+00  8.89999833e-01]
 [ 5.66099965e+00  4.22100024e+00 -1.23200016e+00]
 [ 5.12300001e+00  4.52100001e+00 -2.13100018e+00]
 [ 6.63000001e+00  4.71899987e+00 -1.20600014e+00]
 [ 5.80899986e+00  3.14100007e+00 -1.24100027e+00]
 [ 4.71299980e+00  6.12900021e+00 -5.37333549e-08]
 [ 3.60100007e+00  6.65299988e+00  0.00000000e+00]
 [ 5.84600025e+00  6.83500005e+00 -5.37333572e-08]
 [ 6.73700002e+00  6.35900019e+00 -9.92665958e-08]
 [ 5.84600030e+00  8.28400041e+00  1.02478938e-08]
 [ 4.81899988e+00  8.64799981e+00  5.45776679e-08]
 [ 6.36000027e+00  8.64799972e+00  8.89999998e-01]
 [ 6.36000022e+00  8.64799980e+00 -8.89999974e-01]]

[AnirG]

@JIMonroe Looks like it. Apologies, trying it fresh gave a similar output as yours for the first transformation (probably jupyter messup). I will try to work on this issue and come up with a plausible fix.

[JIMonroe]

Great, thanks, @AnirG !

I think the simplest fix is what I described above, adding copy.deepcopy() when defining the torsions from the input data.

[Neel-Shah-29]

I would like to solve this issue, if no one is currently working on it.

[umak1106]

Hello @orbeckst , is @Neel-Shah-29 still working on this issue ? Or can I try coming up with a possible solution for this issue ?

[orbeckst]

PR #3563 has been closed so you're more than welcome to tackle this issue, @umak1106 .