Create qube_to_PrmRst.py

[SofiaBariami]

Reads qubekit output files (pdb/xml) and converts them to amber format (prm7/rst7) to use them with Sire

usage: ~/sire.app/bin/qube_to_PrmRst -p PDB_FILE -x XML_FILE

[lohedges]

Hi there,

I'll take a closer look tomorrow, but since this a standalone script I can't foresee any issues. Just a couple of quick comments:

• Do you have any tests? Even a simple single point energy calculation with Sire or sander for a known system might be useful to ensure that nothing gets broken going forward.
• The script name is a bit of tongue twister and would be a pain to type on the command-line (without autocompletion). How about qube_convert? Are there specifics that mean that it can only be converted to AMBER format, or could it also be converted to, e.g. GROMACS (even if you go via AMBER files as an intermediate). Perhaps a more general name might be useful if you do choose to expand functionality later.

Cheers.

[jmichel80]

hi @SofiaBariami

picking up on Lester's comments

1. can you also commit a unit test showing a consistent single point energy calculations between a prm7 created by your script from a pdb/xml input and a hardcoded reference energy value calculated previously for that system ? If that's not too hard you can add a dozen small molecules from your hydration free energy study. I would also like to see examples that mix vacuum and solvated systems to test that the space/periodicity treatment is correct.

2. yes needs a better name

3. Regarding output file format, that's a very good point. Would this work to create gromacs input ? If so it would be neat to give this option. Later we might think about your script into a general topology conversion utility accessible directly from BioSimSpace.

[SofiaBariami]

@lohedges , sorry I didn't see the notification earlier, good points!
@jmichel80 yes, no problem, I will add a unit test. So far we have only converted molecules in vacuum, so it would be good to see how it behaves with solvated systems. I will also see if we can generate GROMACS topologies. Will let you know later today, or on Monday.

[SofiaBariami]

I wrote a unit test that generates amber files, calculates the energies using the amber and qube files separately and then uses assert_almost_equal() to compare the energies. You can find it here.

Regarding the solvated system, the xml files that are available contain parameters just for one molecule, so I'm not sure we can test this case. QUBEKit generates parameters for individual molecules either way.

I was also looking at the documentation of Sire to find a way to convert to gro/top. Sire.IO.GroTop is the way to do it, but we're missing a property: The forcefield described by 'MMDetail::null' does not have a property called 'vdwstyle'.
The Sire MM Library should be modified to include the oplsaa forcefield. Would you like me to go ahead and find a way to do it using Sire? We can probably do it with BioSimSpace more easily.

[lohedges]

Hi @SofiaBariami, I should hopefully be able to take a proper look at this tomorrow or Friday.

I wrote a unit test that generates amber files, calculates the energies using the amber and qube files separately and then uses assert_almost_equal() to compare the energies. You can find it here.

Ideally the test would run your script as a subprocess, rather than re-implementing everything in the test. That way you only ever need to update things in one place, i.e. the original script, and it also serves as a test of the way the user would interact with it from the command-line.

I was also looking at the documentation of Sire to find a way to convert to gro/top. Sire.IO.GroTop is the way to do it, but we're missing a property: The forcefield described by 'MMDetail::null' does not have a property called 'vdwstyle'.

I assume this is because you are trying to convert from the Sire System you create directly to the GroTop, so you are missing some properties that are required by the parser. You could just manually create an appropriate MMDetail object and add that as a molecule property before passing the system to the GroTop constructor. Alternatively, you could see what happens when you create an AMBER system (like you currently do) then convert that system to GROMACS. Assuming the MMDetail created by the AMBER parser is compatible it should work okay. There are some issues with scale factors being lost on conversion with OPLS-style force fields, however. See this issue thread for a discussion and workaround.

[lohedges]

Hi @SofiaBariami. Just to note that I've moved your QUBE tests to a separate feature branch here and deleted them from devel. If you need to make edits, please do so on the feature branch then issue a pull request to the devel branch when you're happy. Since any commit to Sire's devel branch runs unit tests against the devel branch of SireUnitTests, they would fail if we happen to add more commits before accepting this pull request. (Since the QUBE conversion code wouldn't yet be included, although since you've reproduced the entire script in the test it might be okay.)

General practice is to create a matching feature branch on SireUnitTests when there are tests that accompany a pull request. That way I can check that they pass locally, then merge the SireUnitTests pull request, then finally merge the Sire pull request itself.

Cheers.