swap_existing_ligand_parameters with OPC water

[emainas]

Has this been tested with OPC water. I keep getting

AmberWarning: Molecule atoms are not contiguous! Attempting to reorder to fix.
warn('Molecule atoms are not contiguous! Attempting to reorder to fix.', AmberWarning)

after combing my ligand structure with the opc structure.

[mattwthompson]

Could you describe, or share code for, how you're using OPC?

[emainas]

I created prmtop and rst7 using tleap:

source leaprc.protein.ff19SB
source leaprc.gaff
source leaprc.water.opc

BIL = loadmol2 biliverdin.mol2
loadamberparams biliverdin.frcmod
list
check BIL

solvatebox BIL OPCBOX 10.0
setbox BIL vdw
saveamberparm BIL gaff_opc_vdw.prmtop gaff_opc_vdw.rst7

then I used these two files as input into the notebook to make the parmed structure:

in_prmtop = "gaff_opc_vdw.prmtop" 
in_crd = "gaff_opc_vdw.rst7" 
orig_structure = parmed.amber.AmberParm(in_prmtop, in_crd)

pieces = orig_structure.split()

ligand_off_molecule = Molecule("biliverdin.sdf") 
ligand_pdbfile = PDBFile("biliverdin.pdb") 
ligand_off_topology = Topology.from_openmm(
    ligand_pdbfile.topology,
    unique_molecules=[ligand_off_molecule],
)
force_field = ForceField("openff_unconstrained-2.0.0.offxml")
ligand_system = force_field.create_openmm_system(ligand_off_topology)

new_ligand_structure = parmed.openmm.load_topology(
    ligand_off_topology.to_openmm(),
    ligand_system,
    xyz=pieces[0][0].positions,
)

just_water_structure = parmed.Structure()
just_water_structure += pieces[1][0]
just_water_structure *= len(pieces[1][1])

new_ligand_structure += parmed.amber.AmberParm.from_structure(just_water_structure)
new_ligand_structure.coordinates = orig_structure.coordinates
new_ligand_structure.box_vectors = orig_structure.box_vectors

new_ligand_structure.save("test.rst7", overwrite=True)
new_ligand_structure.save("test.prmtop", overwrite=True)

[emainas]

parmed_ligand.zip

All the files are here

[mattwthompson]

Thanks for the reproduction and files - I will try to dig into it more later, but a couple things stand out to me at surface-level, which is all I have time for immediately:

• I'm not sure if ParmEd well-supports ff19SB, because of CMAPs. It may now, but it did not a while back when I last checked.
• I think that warning is emitted by ParmEd itself, so we might not be able to fix things here
• I notice it's a warning, not an error, so I wonder if it's actually a problem or not

All that being said - is your goal to run a protein-ligand complex in Amber using Sage, OPC and ff19SB for the ligand, water, and protein, respectively? If you can substitute ff14SB instead ¹ you can do this all within OpenFF tools, and I'd point you to this example instead. I think the "swap existing parameters" example is there to show it's possible, not that it's necessarily the best way to prepare a system.

https://docs.openforcefield.org/en/latest/examples/openforcefield/openff-interchange/protein_ligand/protein_ligand.html

Footnotes

1. When Rosemary (OpenFF's 3.0 force field line) is released, there'll be no need to use separate force fields for the ligand and protein! ↩

[emainas]

Thanks for the quick response Matt. The endgoal it to add the protein but for now I just have a ligand in water. Sage for the ligand and OPC for the water. Indeed it is a warning, but it cannot be neglected because in the output parameter file it duplicates EP (extra points that the opc water molecules have). It does it right after:

new_ligand_structure.coordinates = orig_structure.coordinates
new_ligand_structure.box_vectors = orig_structure.box_vectors

but I do not understand why.

Also I forgot to mention, the callback comes from _amberparm.py in case you have a look at it. I did exactly the same procedure with TIP3P water (has no EPs) and it works just fine so I assume it is an EP issue.

I will look into Interchange, thank you for the suggestion.

[mattwthompson]

Sage with OPC ligands should work fine in OpenMM, and maybe in GROMACS but I forget right now. Interchange doesn't yet support writing OPC water to Amber, but I'm trying to figure out how to implement it (tracking: openforcefield/openff-interchange#783).