[BUG] BSS IO flipping the box angle

[xiki-tempula]

Describe the bug
I have an amber.crd file. If I use cpptraj to convert it to RST7 file, I got the angle
120.0000229 120.0000229 90.0000000
If I use BSS to read in the file and write it out as RST7 file, I got the angle
59.9999771 59.9999771 90.0000000

To Reproduce
To generate the rst7 from cpptraj

cpptraj -p amber.prm7
trajin amber.crd
trajout amber_out.rst7
run

To generate the rst7 from BSS

import BioSimSpace.Sandpit.Exscientia as BSS
top = BSS.IO.readMolecules(['amber.prm7', 'amber.crd'])
BSS.IO.saveMolecules('BSS_out', top, 'RST7')

Expected behavior
They are the same.

Input files
Archive.zip

[lohedges]

This is a result of the lattice reduction performed on the triclinic lattice vectors to make sure that they are consistent with all engines that we support. The two boxes are equivalent.

[xiki-tempula]

AMBER doesn't think they are equivalent. If you run a amber MPI run with replica exchange, amber would quit saying that not all replic has the same box shape/angle.

[lohedges]

Presumably because the AMBER code is simply checking the exact numerical values for the box dimensions and angles, not whether it is an equivalent triclinic space.

If possible, can you give some insight as to how you are setting up your replicas. Even if Sire/BioSimSpace are transforming things, I would assume that the transformation is the same for each replica, so the space would be the same. I can understand the result of a rounding issue here, since this could (potentially) differ between replicas.

[xiki-tempula]

So we have one lambda window where the simulation generates the output file as 120.0000229 120.0000229 90.0000000. When you do process.getSystem(), and then use Process to write a new RST7 file. The new RST7 file has 59.9999771 59.9999771 90.0000000.

For another lambda window, the simulation generates the output file as 120.0000229 120.0000229 90.0000000. When you do process.getSystem(), and then use Process to write a new RST7 file. the new RST7 file retains the input box dimension of 120.0000229 120.0000229 90.0000000, which is then read by AMBER and cause the error.

[lohedges]

Are you able to share the two RST7 files that were created by the two lambda windows? I want to check to see whether reading then writing them back produces different box angles. Assuming the same info in the RST7 files, the lattice reduction should be the same, so it's strange that you aren't getting the same output. Given the limited precision in the file, I'm wondering if some of the conditionals used in the reduction are giving different results due to numerical precision. I can also try reading/writing one repeatedly to make sure that the output is the same.

For your original input here I can read and write the file repeatedly and get the exact same output, i.e.:

import BioSimSpace as BSS

for x in range(0, 100):
    s = BSS.IO.readMolecules(["amber.prm7", "amber_out.rst7"])
    BSS.IO.saveMolecules(f"tmp_{x:03d}", s, "rst7")

Then:

tail -q -n 1 tmp_*rst7 > box.txt
uniq box.txt
  81.5793091  81.5793091  81.5792526  59.9999771  59.9999771  90.0000000

So, the exact same box information is written each time for the same input. (Note that I'm loading a different system each time, so it's not using the cached file.)

[xiki-tempula]

Archive.zip

mol_0 = BSS.IO.readMolecules(['lambda_0.0000/amber.prm7', 'lambda_0.0000/amber.crd'])
BSS.IO.saveMolecules('mol_0', mol_0, 'RST7')
mol_1 = BSS.IO.readMolecules(['lambda_1.0000/amber.prm7', 'lambda_1.0000/amber.crd'])
BSS.IO.saveMolecules('mol_1', mol_1, 'RST7')

You could see the first one is 59.9999771 59.9999771 90.0000000 and the second one is 120.0000229 120.0000229 90.0000000

One could also do the same conversion with amber

cpptraj -p lambda_0.0000/amber.prm7 -y lambda_0.0000/amber.crd -x lambda_0.rst7
cpptraj -p lambda_1.0000/amber.prm7 -y lambda_1.0000/amber.crd -x lambda_1.rst7

which gives the same box angle

[lohedges]

Okay, I see the issue. Although your box angles are the same in the input written for the lambda legs, the box dimensions are not, so the lattice reduction code is giving a different result on read due to numerical imprecision, even when using binary format for the coordinates.

In [4]: mol_0._sire_object.property("space")
Out[4]: TriclinicBox( ( 81.5793, 0, 0 ), ( 4.99529e-15, 81.5793, 0 ), ( 40.7897, 40.7897, 57.6852 ) )

In [5]: mol_1._sire_object.property("space")
Out[5]: TriclinicBox( ( 81.4219, 0, 0 ), ( 4.98565e-15, 81.4219, 0 ), ( -40.711, -40.711, 57.5739 ) )

I've modified the AMBER crd code to actually print out what's in the two files. The first line is the x, y, z box dimensions, the second line is the angles:

mol_0 = BSS.IO.readMolecules(['lambda_0.0000/amber.prm7', 'lambda_0.0000/amber.crd'])
81.5793 81.5793 81.5793
120 120 90

mol_1 = BSS.IO.readMolecules(['lambda_1.0000/amber.prm7', 'lambda_1.0000/amber.crd'])
81.4219 81.4219 81.4218
120 120 90

(It looks like the box precision isn't any higher in the DCD file.)

Looking at the bottom of the rst7 files gives the same thing:

tail -n 1 lambda_*/amber.rst7
==> lambda_0.0000/amber.rst7 <==
  81.5793091  81.5793091  81.5792527 120.0000229 120.0000229  90.0000000

==> lambda_1.0000/amber.rst7 <==
  81.4219024  81.4219024  81.4218459 120.0000229 120.0000229  90.0000000

Presumably these have been equilibrated independently using NPT.

[lohedges]

I'll close this since it's the same as OpenBioSim/sire#49, which I'll transfer to Sire.