openforcefield / openff-toolkit

The Open Forcefield Toolkit provides implementations of the SMIRNOFF format, parameterization engine, and other tools. Documentation available at http://open-forcefield-toolkit.readthedocs.io
http://openforcefield.org
MIT License
302 stars 88 forks source link

swap_existing_ligand_parameters with OPC water #1769

Open emainas opened 7 months ago

emainas commented 7 months ago

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 commented 7 months ago

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

emainas commented 7 months ago

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 commented 7 months ago

parmed_ligand.zip

All the files are here

mattwthompson commented 7 months ago

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:

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 [^1] 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

[^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 commented 7 months ago

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 commented 7 months ago

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: https://github.com/openforcefield/openff-interchange/issues/783).