Refactor wrapper and add Gromacs support

[jeff231li]

This PR refactors the Amber simulation wrapper module paprika/amber.py and adds a new wrapper for Gromacs. Simulation wrappers are now stored in paprika/simulate.py.

from paprika.simulate import Amber
simulation = Amber()

from paprika.simulate import Gromacs
simulation = Gromacs()

Also added a small function in tleap.System to convert Amber prmtop/rst7 to Gromacs top/gro format.

Question: do you think it's better to have all caps for the simulation wrapper class name?

EDIT: I've made simulation class names all caps.

TODO

• Compare Amber and Gromacs results for CB6-BUT-TIP3P system.

Status

• Ready to go

[codecov-commenter]

Codecov Report

Merging #140 into master will decrease coverage by 2.92%.
The diff coverage is 39.92%.

[slochower]

Do we have to change the openff-evaluator paprika protocol in light of this?

[slochower]

I didn't run through the tutorial, so I trust that it works.

However, it struck me that perhaps we should create a simplified workflow because not every step is obvious. This could live in the documentation. For example, I noticed after converting the System to Gromacs, format you do:

    # Add dummy atom restraints - load in the Gromacs files instead 
    # of the amber files because the coordinates a shifted.
    structure = pmd.load_file(
        os.path.join(work_dir, window, f"{base_name}-sol.top"),
        xyz=os.path.join(work_dir, window, f"{base_name}-sol.gro"), 
        structure=True
    )
    plumed.add_dummy_atoms_to_file(structure,window=window)

...and I was thinking to myself, how would someone know to do this? It seems like the dummy atoms are in the structure already, because tleap is adding them, but I guess they still need to be added? (Also, should add_dummy_atoms_to_file be renamed to add_dummy_atoms_to_structure)?

Maybe the workflow could look something like:

1. Setup the complex for the calculation by creating a bound pose and aligning it.
2. Determine the number of windows for the umbrella sampling (attach, pull, and release phases of the calculation).
3. Setup the restraints.
4. Create the windows.
5. Write the restraints to a file in each window (either AMBER-format or PLUMED-format).
6. Create coordinates for all windows.
7. Solvate the complex in each window and add dummy atoms. (Is this how you'd phrase things?)
8. If using PLUMED, add dummy atoms to the restraint file (??)
   [...etc...]

What do you think?

[jeff231li]

Do we have to change the openff-evaluator paprika protocol in light of this?

I don't think we need to change the openff-evaluator protocol, this protocol is self-contained. I think if we change analysis.py and restraints/restraints.py then we will have to adapt the openff-evaluator protocol.

[jeff231li]

(Also, should add_dummy_atoms_to_file be renamed to add_dummy_atoms_to_structure)?

Hmm, maybe add_dummy_atom_restraints would be more accurate as we are adding positional restraints on the dummy atoms.

However, it struck me that perhaps we should create a simplified workflow because not every step is obvious
...and I was thinking to myself, how would someone know to do this?

Yeah, I agree that the workflow/protocol of doing APR calculations with pAPRika is not clear. I am planning on adding more info to the documentation and maybe have a flowchart diagram to demonstrate this.

7. Solvate the complex in each window and add dummy atoms. (Is this how you'd phrase things?)

The dummy atoms are already added with tleap it's just that the positional restraints on the dummy atoms require absolute coordinates, unlike all the other restraints. When tleap solvates the host-guest-dum system in vacuum, the program shifts the coordinates. This is not a problem when we run Amber simulations because we take care of the dummy atom restraints in the simulation config file.

[slochower]

Hmm, maybe add_dummy_atom_restraints would be more accurate as we are adding positional restraints on the dummy atoms.

I think that is a bit more explicit.

Yeah, I agree that the workflow/protocol of doing APR calculations with pAPRika is not clear. I am planning on adding more info to the documentation and maybe have a flowchart diagram to demonstrate this.

+1

The dummy atoms are already added with tleap it's just that the positional restraints on the dummy atoms require absolute coordinates, unlike all the other restraints. When tleap solvates the host-guest-dum system in vacuum, the program shifts the coordinates. This is not a problem when we run Amber simulations because we take care of the dummy atom restraints in the simulation config file.

This makes sense. I think I didn't totally appreciate the distinction before.

[slochower]

The binding free-energy for CB6-BUT-TIP3P:

@slochower LGTM?

I'll give it another once-over tomorrow.

[slochower]

Do we need to upload the tutorial/complex files? That's going to make this repo heavyweight. Can we not just load the relevant files from ../data and not save tutorial-generated files?

[slochower]

Looks pretty good. I added a bunch of comments about minor wording changes, mostly in the documentation.

[jeff231li]

Do we need to upload the tutorial/complex files? That's going to make this repo heavyweight. Can we not just load the relevant files from ../data and not save tutorial-generated files?

Yeah, I can remove the complex folder from the repository and the tutorials can use the files from paprika/data

[codecov-io]

Codecov Report

Merging #140 into master will decrease coverage by 2.92%.
The diff coverage is 39.92%.