michellab / Sire

Sire Molecular Simulations Framework
http://siremol.org
GNU General Public License v3.0
95 stars 26 forks source link

Velocities lost during read/write with Gro87 parser #387

Closed kexul closed 2 years ago

kexul commented 2 years ago

Dear BioSimSpace developers: I found that velocities in Gromacs .gro file would lose during reading and saving by BioSimSpace, is this a bug or feature?

Steps to reproduce:

import BioSimSpace as BSS
sys = BSS.IO.readMolecules(['gromacs.gro', 'gromacs.top'])
BSS.IO.saveMolecules('new', sys, ['gro87', 'grotop'])

Inputs: temp.zip

Diff: d5d69e82b417f9e77f6eb77d619d853

lohedges commented 2 years ago

I see mention of velocities in the Gro87 parser, so it looks like it should be possible to read/write them. It could be something to do with the numerical precision, in which case it's possible to specify the formatting using the property map.

I'll see if I can reproduce here and figure out why it's not working.

Cheers.

lohedges commented 2 years ago

I can't even read the example files that you provided:

import BioSimSpace as BSS
BSS.setVerbose(True)

sys = BSS.IO.readMolecules(['gromacs.gro', 'gromacs.top'])

gives...

...
/home/lester/sire.app/bin/python(+0x1be51d) [0x556cb724551d]
__EndTrace__
Exception 'SireError::io_error' thrown by the thread 'master'.
There is not enough information in this parser (Gro87( title() = GROtesk MACabre and Sinister, nAtoms() = 27716, nResidues() = 7996, nFrames() = 1, hasCoordinates() = 1, hasVelocities() = 1 )) to start the creation of a new System. You need to use a more detailed input file.
Thrown from FILE: /home/lester/Code/Sire/corelib/src/libs/SireIO/moleculeparser.cpp, LINE: 1466, FUNCTION: virtual SireSystem::System SireIO::MoleculeParser::startSystem(const QVector<QString>&, const SireBase::PropertyMap&) const

The above exception was the direct cause of the following exception:

OSError                                   Traceback (most recent call last)
<ipython-input-8-74aefe221fd5> in <module>
----> 1 sys = BSS.IO.readMolecules(['gromacs.gro', 'gromacs.top'])

~/Code/BioSimSpace/python/BioSimSpace/IO/_io.py in readMolecules(files, property_map)
    395                 msg = "Failed to read molecules from: %s" % files
    396                 if _isVerbose():
--> 397                     raise IOError(msg) from e
    398                 else:
    399                     raise IOError(msg) from None

OSError: Failed to read molecules from: ['gromacs.gro', 'gromacs.top']

Reading the files individually, it is the top file that fails.

from Sire.IO import GroTop

gro_top = GroTop("gromacs.top")

gives...


/home/lester/sire.app/bin/python(+0x1be51d) [0x556cb724551d]
__EndTrace__
Exception 'SireError::io_error' thrown by the thread 'master'.
Cannot find the file '/data/FEP/prepare/test_md/./npt2/bound/posre_0001.itp'. Please make sure that this file exists and is readable
Thrown from FILE: /home/lester/Code/Sire/corelib/src/libs/SireIO/grotop.cpp, LINE: 4989, FUNCTION: QString SireIO::GroTop::findIncludeFile(QString, QString)

It looks like you are including position restraint files that aren't included here.

I'll just try working with the Gro87 file, since (I believe) that's where the problem lies.

kexul commented 2 years ago

Sorry about that 😰, there are the restraint files: temp.zip

lohedges commented 2 years ago

This is strange. The velocities are written back to the file, but they are wrong/different (as you show) and aren't read back in at all. I'll try sorting the file by the velocities to see if they are the same, but being written out of order. Alternatively, there could be a unit conversion error. I'll also move this over to the Sire repository since it's a bug with the Sire.IO.Gro87 parser, not BioSimSpace.

lohedges commented 2 years ago

Okay, I've now fixed this. There were numerous incorrect unit conversions on reading / storing / extracting / writing the velocities. In some cases, (different) unit conversions were applied multiple times and the forward / reverse (read / write) procedure was also not symmetric.

Will close when I push to devel. Note that we still haven't resolved the issue with our macOS pipeline, so there won't be an updated conda package until we figure out a fix.

lohedges commented 2 years ago

Looks like some unit tests failed due to precision errors in comparisons. Will fix those when I get a chance and re-push. Ironically they passed before, despite the numbers being complete nonsense.

lohedges commented 2 years ago

Actually, there was an incorrect unit conversion in the test too, i.e. going from AMBER to GROMACS units. I think this was the original issue too, i.e. the GROMACS code had been copied from the AMBER parser, with various unit conversions not replaced.