openff-toolkit

openff-toolkit copied to clipboard

I am trying to use ParmEd v3.4.3 to convert a system (with a "sage" ligand) solvated in tip4pew water from openMM v. 7.7.0 to GROMACS topology and configuration files (gro). Python lines doing the conversion are reported below (I can provide the full script with generation of openMM system + sage ligand + solvation if needed). Both openMM "export_system" and "pmd_complex_struct" seem to correctly include the tip4pew virtual sites. However, the generated GROMACS topology file does not contain the virtual site and the generation of the "gro" file is interrupted with this error

"RuntimeError: Could not find <Atom O [0]; In HOH 0>"

and no virtual sites are present in the "gro" file written before the error. Any tips?

ParmEd's GROMACS exporter can't handle constraints from openmm, so we need a variant for export without them

export_system = omm_forcefield.createSystem( modeller.topology, nonbondedMethod=openmm.app.PME, constraints=None, rigidWater=False, )

Combine the topology, system and positions into a ParmEd Structure

pmd_complex_struct = pmd.openmm.load_topology(topology, export_system, positions)

Export GROMACS files

pmd_complex_struct.save("system.top", overwrite=True) pmd_complex_struct.save("system.gro", overwrite=True)

Hi @maxbonomi, unfortunately we don't maintain ParmEd so there's not much we can do. I haven't personally used virtual sites with ParmEd, so I can't be of much more help there - maybe somebody else in the community has and will comment here or on the issue you raised over there. I see you've already started hacking into the ParmEd source code to get this working, though, and if you want to stick with using GROMACS I don't think there is any path forward besides that one.

We are working on a tool that fits into our infrastructure in the ways that ParmEd historically has. However, exports to GROMACS with virtual sites are not sufficiently supported and might not be for a couple of months.

Something that caught my eye, however, is that you're calling omm_forcefield.createSystem(), something that the OpenFF Toolkit's ForceField objeect doesn't support ... are you using openmmforcefields elsewhere in your pipeline?

Yes of course. I am using sage for my ligand, plus I have protein and sometimes RNA. I need to create a new system with "omm_forcefield.createSystem" without constraints and rigid water to pass to ParmEd. This is the procedure I found on ParmEd website (if I remember correctly) and works nicely with 3p waters.

@mattwthompson "However, exports to GROMACS with virtual sites are not sufficiently supported and might not be for a couple of months." Happy to test an alpha version and provide feedback if you want.

Wonderful - I'll try to put something together in the next couple of days and ping you when it's ready for a test spin!

Super, thanks!