search

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

topology for drug like molecule #1793

Closed Lili-Cao closed 6 months ago

Lili-Cao commented 6 months ago

I want to create a topology for a drug like molecule which is not in the 20 canonical amino acids

Expand the from_pdb method to accept drug like molecule

Create a new method to read a drug like molecule

**I tried to use pass a sdf file to Molecule from_files() and create a topology by from_molecules but get several complains, when convert it to gromacs files.

But I get several error message when validating the intercange file for example: get_gromacs_energies(interchange).energies UnsupportedExportError: Electrostatics method Coulomb not supported**

mattwthompson commented 6 months ago

Hi @Lili-Cao,

I'm going to move this to the repo where the issues are bubbling out from. I'll also try to give some more context:

from openff.units import unit, Quantity

# load up molecules, make a topology (Topology object)

...
toplogy.box_vectors = Quantity([4, 4, 4], unit.nanometer)
...

# write out to GROMACS files

If there's more issues you run into, I'm happy to help. I would need to see more code so I can run the same thing on my computer, though.

Lili-Cao commented 6 months ago

Hi Matt,

Thanks for your prompt response.

I tried with two things. Building a Topology object from a pdb file and from a sdf file.

below is the example for the sdf file, please let me know if you see any suspicious things Thank you! LC from openmm import app from openff.toolkit import Molecule, Topology from openff.toolkit.utils.utils import get_data_file_path

sdf_filepath = "xxxx/sdf/xxxx-N.sdf" unique_molecules = [Molecule.from_file(sdf_filepath)] topology = Topology.from_molecules(unique_molecules)

forcefield = ForceField("openff-2.1.0.offxml") omm_system = forcefield.create_openmm_system(topology)

from openff.interchange import Interchange

interchange = Interchange.from_smirnoff( force_field=forcefield, topology=topology, )

Export GROMACS files.

interchange.to_top("system.top") interchange.to_gro("system.gro") from openff.interchange.drivers import get_amber_energies, get_openmm_energies, get_gromacs_energies, get_lammps_energies gromacs_energies = get_gromacs_energies(interchange) gromacs_energies.energies 175 mdp.write(f"rcoulomb = {coul_cutoff}\n") 176 else: --> 177 raise UnsupportedExportError( 178 f"Electrostatics method {self.coul_method} not supported", 179 ) 181 if self.vdw_method == "cutoff": 182 mdp.write("vdwtype = cutoff\n") UnsupportedExportError: Electrostatics method Coulomb not supported

Den mån 18 dec. 2023 kl 14:07 skrev Matt Thompson @.***

:

Hi @Lili-Cao https://github.com/Lili-Cao,

I'm going to move this to the repo where the issues are bubbling out from. I'll also try to give some more context:

from openff.units import unit, Quantity

load up molecules, make a topology (Topology object)

...toplogy.box_vectors = Quantity([4, 4, 4], unit.nanometer) ...

write out to GROMACS files

If there's more issues you run into, I'm happy to help. I would need to see more code so I can run the same thing on my computer, though.

— Reply to this email directly, view it on GitHub https://github.com/openforcefield/openff-toolkit/issues/1793#issuecomment-1860454602, or unsubscribe https://github.com/notifications/unsubscribe-auth/ASHJ3CYI6OFHNBJWV2RFW43YKA52PAVCNFSM6AAAAABAZHZGT6VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTQNRQGQ2TINRQGI . You are receiving this because you were mentioned.Message ID: @.***>

Lili-Cao commented 6 months ago

Here is the sdf file I was referring to.

BR LC

Den mån 18 dec. 2023 kl 14:19 skrev Lili Cao @.***>:

Hi Matt,

Thanks for your prompt response.

I tried with two things. Building a Topology object from a pdb file and from a sdf file.

below is the example for the sdf file, please let me know if you see any suspicious things Thank you! LC from openmm import app from openff.toolkit import Molecule, Topology from openff.toolkit.utils.utils import get_data_file_path

sdf_filepath = "xxxx/sdf/xxxx-N.sdf" unique_molecules = [Molecule.from_file(sdf_filepath)] topology = Topology.from_molecules(unique_molecules)

forcefield = ForceField("openff-2.1.0.offxml") omm_system = forcefield.create_openmm_system(topology)

from openff.interchange import Interchange

interchange = Interchange.from_smirnoff( force_field=forcefield, topology=topology, )

Export GROMACS files.

interchange.to_top("system.top") interchange.to_gro("system.gro") from openff.interchange.drivers import get_amber_energies, get_openmm_energies, get_gromacs_energies, get_lammps_energies gromacs_energies = get_gromacs_energies(interchange) gromacs_energies.energies 175 mdp.write(f"rcoulomb = {coul_cutoff}\n") 176 else: --> 177 raise UnsupportedExportError( 178 f"Electrostatics method {self.coul_method} not supported", 179 ) 181 if self.vdw_method == "cutoff": 182 mdp.write("vdwtype = cutoff\n") UnsupportedExportError: Electrostatics method Coulomb not supported

Den mån 18 dec. 2023 kl 14:07 skrev Matt Thompson < @.***>:

Hi @Lili-Cao https://github.com/Lili-Cao,

I'm going to move this to the repo where the issues are bubbling out from. I'll also try to give some more context:

from openff.units import unit, Quantity

load up molecules, make a topology (Topology object)

...toplogy.box_vectors = Quantity([4, 4, 4], unit.nanometer) ...

write out to GROMACS files

If there's more issues you run into, I'm happy to help. I would need to see more code so I can run the same thing on my computer, though.

— Reply to this email directly, view it on GitHub https://github.com/openforcefield/openff-toolkit/issues/1793#issuecomment-1860454602, or unsubscribe https://github.com/notifications/unsubscribe-auth/ASHJ3CYI6OFHNBJWV2RFW43YKA52PAVCNFSM6AAAAABAZHZGT6VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTQNRQGQ2TINRQGI . You are receiving this because you were mentioned.Message ID: @.***>

Lili-Cao commented 6 months ago

Hi,

Quick updates, after add a box with 4, 4, 4, I got a periodic Topology object which went through get_gromacs_enegies.

Do you know if it is possible to generate a position restraint file also? which will be posre_xx.itp file with the content looks like bellow: [image: image.png]

Den mån 18 dec. 2023 kl 14:07 skrev Matt Thompson @.***

:

Hi @Lili-Cao https://github.com/Lili-Cao,

I'm going to move this to the repo where the issues are bubbling out from. I'll also try to give some more context:

from openff.units import unit, Quantity

load up molecules, make a topology (Topology object)

...toplogy.box_vectors = Quantity([4, 4, 4], unit.nanometer) ...

write out to GROMACS files

If there's more issues you run into, I'm happy to help. I would need to see more code so I can run the same thing on my computer, though.

— Reply to this email directly, view it on GitHub https://github.com/openforcefield/openff-toolkit/issues/1793#issuecomment-1860454602, or unsubscribe https://github.com/notifications/unsubscribe-auth/ASHJ3CYI6OFHNBJWV2RFW43YKA52PAVCNFSM6AAAAABAZHZGT6VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTQNRQGQ2TINRQGI . You are receiving this because you were mentioned.Message ID: @.***>

mattwthompson commented 6 months ago

No, the toolkit nor Interchange do not store restraints

Lili-Cao commented 6 months ago

OK, May I ask in which way I could get it in openff?

Best Regards Lili

Den mån 18 dec. 2023 kl 14:58 skrev Matt Thompson @.***

:

No, the toolkit nor Interchange do not store restraints

— Reply to this email directly, view it on GitHub https://github.com/openforcefield/openff-toolkit/issues/1793#issuecomment-1860571934, or unsubscribe https://github.com/notifications/unsubscribe-auth/ASHJ3C5FWY6W7DHPGPTQBSLYKBDXTAVCNFSM6AAAAABAZHZGT6VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTQNRQGU3TCOJTGQ . You are receiving this because you were mentioned.Message ID: @.***>

mattwthompson commented 6 months ago

You'll have to generate that file externally, our models don't store that information

Lili-Cao commented 6 months ago

ok. Thanks.

Do you have an example the the generate .grp and .top file works for gromacs? In the converted top file, some info are quite different comparing to acpyp generated which I have confirmed works well.

Best Regards Lili

Den mån 18 dec. 2023 kl 15:21 skrev Matt Thompson @.***

:

You'll have to generate that file externally, our models don't store that information

— Reply to this email directly, view it on GitHub https://github.com/openforcefield/openff-toolkit/issues/1793#issuecomment-1860629479, or unsubscribe https://github.com/notifications/unsubscribe-auth/ASHJ3C7YELWGE5J6Y5A7FKDYKBGOBAVCNFSM6AAAAABAZHZGT6VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTQNRQGYZDSNBXHE . You are receiving this because you were mentioned.Message ID: @.***>

mattwthompson commented 6 months ago

Unfortunately I don't know what a .grp file is; your example above shows how to generate the GROMACS topology file (Interchange.to_top, or Interchange.to_gromacs)

Lili-Cao commented 6 months ago

I meant .gro, sorry There is a jupyter notebook available for generating a gromacs .gro and .top file. Do you also have an example using this two files which works with Gromacs?

Den mån 18 dec. 2023 kl 16:40 skrev Matt Thompson @.***

:

Unfortunately I don't know what a .grp file is; your example above shows how to generate the GROMACS topology file (Interchange.to_top, or Interchange.to_gromacs)

— Reply to this email directly, view it on GitHub https://github.com/openforcefield/openff-toolkit/issues/1793#issuecomment-1860843920, or unsubscribe https://github.com/notifications/unsubscribe-auth/ASHJ3C4QPRWQNJAAVRCKIA3YKBPWFAVCNFSM6AAAAABAZHZGT6VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTQNRQHA2DGOJSGA . You are receiving this because you were mentioned.Message ID: @.***>

mattwthompson commented 6 months ago

Those two files do work with GROMACS, or at least serve as a starting point depending on what you want to do. OpenFF's objective is to parameterize chemical topologies with force fields and provide users the starting point of running simulations, not going all the way through production runs and trajectory analysis. GROMACS can do many, many things and most of them are outside the scope of what one tool could possibly handle. To continue with your workflow, you likely want to consult other resources like the official GROMACS documentation (and probably other tools like https://mdanalysis.org/).

Lili-Cao commented 6 months ago

To make it clear, I only want ff from openFF to parametrize a small molecule and convert the topology to whatever gromacs accept. Everything else I have code for that but don't fully understand the converted .top format, how it should be used by gromacs.

j-wags commented 6 months ago

Hi @Lili-Cao - I think you're asking for more detail based on some attached files, but it doesn't look like your email attachments are coming through. Could you post your attachments using the web browser interface at https://github.com/openforcefield/openff-toolkit/issues/1793?

Lili-Cao commented 6 months ago

Hi Jeff,

I hop this time the sdf file is uploaded successfully. (its with a txt suffix and please change it to sdf, otherwise I cannot upload it) CDDNTCUZCTVFFB-UYBDAZJANA-N.txt

from openmm import app from openff.toolkit import Molecule, Topology from openff.toolkit.utils.utils import get_data_file_path

sdf_filepath = "xxxx/sdf/xxxx-N.sdf" unique_molecules = [Molecule.from_file(sdf_filepath)] topology = Topology.from_molecules(unique_molecules) toplogy.box_vectors = Quantity([4, 4, 4], unit.nanometer)

forcefield = ForceField("openff-2.1.0.offxml") omm_system = forcefield.create_openmm_system(topology)

from openff.interchange import Interchange

interchange = Interchange.from_smirnoff( force_field=forcefield, topology=topology, )

interchange.to_top("system.top") interchaninterchange.to_gro("system.gro")

Till here, I got the system.top and system.gro. As Matt says, posre file is not saved in this module. But the .top file looks different to me, comparing to the .itp or top file generated from tleap, using Acpype, which can be directly used by Gromacs for small compounds. That's why I asking if there is any example workflow includes steps from opeff to gromacs.

I know how to combine those parameterization files for both protein and ligand, run md etc.

Best Regards Lili

j-wags commented 6 months ago

Ok, I've run the following to reproduce your case:

from openmm import app
from openff.toolkit import Molecule, Topology, ForceField
from openff.toolkit.utils.utils import get_data_file_path
from openff.units import unit, Quantity

sdf_filepath = "ligand.sdf"
unique_molecules = [Molecule.from_file(sdf_filepath)]
topology = Topology.from_molecules(unique_molecules)
topology.box_vectors = Quantity([4, 4, 4], unit.nanometer)

forcefield = ForceField("openff-2.1.0.offxml")
#omm_system = forcefield.create_openmm_system(topology)

from openff.interchange import Interchange

interchange = Interchange.from_smirnoff(
force_field=forcefield,
topology=topology,
)

interchange.to_top("system.top")
interchange.to_gro("system.gro")
from openff.interchange.drivers import get_amber_energies, get_openmm_energies, get_gromacs_energies, get_lammps_energies
gromacs_energies = get_gromacs_energies(interchange)
print(f"{gromacs_energies.energies=}")
openmm_energies = get_openmm_energies(interchange)
print(f"{openmm_energies.energies=}")

This prints:

gromacs_energies.energies={'Bond': <Quantity(21.5653725, 'kilojoule / mole')>, 'Angle': <Quantity(7.3002429, 'kilojoule / mole')>, 'Torsion': <Quantity(70.7501363, 'kilojoule / mole')>, 'RBTorsion': <Quantity(0.0, 'kilojoule / mole')>, 'vdW': <Quantity(206.390487, 'kilojoule / mole')>, 'Electrostatics': <Quantity(-432.642303, 'kilojoule / mole')>}
openmm_energies.energies={'Bond': <Quantity(21.5656497, 'kilojoule / mole')>, 'Angle': <Quantity(7.30049334, 'kilojoule / mole')>, 'Torsion': <Quantity(70.7501632, 'kilojoule / mole')>, 'Nonbonded': <Quantity(-226.183232, 'kilojoule / mole')>}

The output makes me think that the resulting energies are the same in both OpenMM and Gromacs. So I think the parameter assignment and energy calculation are happening correctly (but please correct me if I'm wrong).

I'm not sure I understand the comment about the .top and .itp files. I haven't used Gromacs in a little while, but I think itp files were a way to reference multiple .top files at once. Could you clarify how the .top files look different, and how we'd be able to tell whether they're right or wrong?

Also, thanks for your patience and willingness to work on this. We haven't had much feedback on our exporters to GROMACS, so it's helpful to know if there's something wrong, or if we're making files that the users don't expect :-)

Lili-Cao commented 6 months ago

example_files.zip Hi Jeff,

I attached folder which contains 4 example files (from acpype) which I could be used for Gromacs after small modification. The modification usually contains extracting relative part and integrating with those files of a protein.

Best Regards Lili

Yoshanuikabundi commented 6 months ago

Hi @Lili-Cao - I think you're expecting two ITP files for GROMACS so that you can incorporate it into an existing TOP file, one with force field parameters and the other with position restraints. This is the usual way parameters are combined in GROMACS - you import a whole bunch of ITP files in a TOP file, and then conditionally import the position restraints ITP when you want to equilibrate with them.

Unfortunately, supporting this pattern directly is out of scope for Interchange for the moment. Interchange is designed to produce entire systems that are identical across simulation engines. It doesn't produce position restraints at all, because other engines do this in very different ways (though this might be a useful feature @mattwthompson, there is a way to get opt-in, MDP-configurable position restraints into a single TOP file with the pre-processor that could probably be implemented in all GROMACS exports pretty easily - let me know if you'd like a write-up). And it doesn't break its GROMACS exports into individual ITP files - instead, it puts everything directly into the TOP file.

Fortunately @Lili-Cao, the strategy you've taken so far is very much on the right track to accomplishing the above. You just need to re-organise the TOP file into an ITP, and then write a separate ITP with nothing but a position restraint for every atom, and a conditional #IFDEF import. You can model the latter after an existing POSRES file, but unfortunately the former has to be done by hand at the moment, and requires you to understand GROMACS' topology format, which is described in the manual: https://manual.gromacs.org/current/reference-manual/topologies/topologies.html

One thing you'll notice is that to trim the TOP file down to a single [molecule] specification ITP file, you'll have to choose between the [defaults] section from your existing TOP file or the one from Interchange - which will change the force field, either of the small molecule or of the rest of the system. This is one of the many reasons it's a bad idea to combine different force fields - you'll be changing the LJ combination rules and a few other global settings. OpenFF force fields may be compatible with some Amber force fields, but if you're using CHARMM or GROMOS or something else then you won't be getting the parameters we intend, and the non-bonded parameters won't be appropriate for the rest of the system.

I'd recommend making life easier on yourself and preparing the whole system with the Toolkit, using the Amber ff19sb port.

Lili-Cao commented 6 months ago

Hi Josh,

Thanks for your great explanation. It will be great if you can tell me more about the "MDP-configurable position restraints ", put everything into the Top file is perfectly fine for me.

And great to hear that the strategy I've taken is on the right track. Yes, it is not impossible that I can make top and itp files working for Gromacs, just need some struggle : )

Best Regards Lili