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

How to obtain energy components in OpenForceField Toolkit #1752

Closed Guackkk closed 8 months ago

Guackkk commented 8 months ago

Hi!

I've been following the tutorial provided in the OpenForceField Toolkit documentation (https://docs.openforcefield.org/en/latest/examples/openforcefield/openff-toolkit/conformer_energies/conformer_energies.html) to compute conformer energies. It's been a great resource, but I have a question regarding the extraction of energy components.

I'd like to know if it's possible to obtain individual energy components such as bond energy, angle energy, and dihedral energy for each conformer.

Any insights or guidance on this matter would be greatly appreciated!

Yoshanuikabundi commented 8 months ago

Hi! Good question. The notebook in question uses OpenMM to calculate energies, so if you take the simulation object you should be able to use the last cell or two from this notebook in their cookbook: https://openmm.github.io/openmm-cookbook/latest/notebooks/cookbook/Analyzing%20Energy%20Contributions.html

mattwthompson commented 8 months ago

Accessing the Simulation object directly is probably the easiest way to get potential energy components in that notebook.

Another route, which for OpenMM calls the same code, is to use these functions in the openff.interchange.drivers module. I haven't written proper documentation for that module, but here's a snippet and the functions have docstrings:

In [1]: from openff.toolkit import ForceField, Molecule

In [2]: from openff.interchange.drivers import get_openmm_energies

In [3]: molecule = Molecule.from_smiles("c1cc(Cl)ccc1C(=O)CS[C]1=CO[C](F)(F)CC1")

In [4]: molecule.generate_conformers()

In [5]: for index in range(molecule.n_conformers):
   ...:     molecule_copy = Molecule(molecule)
   ...:     molecule_copy._conformers = [molecule.conformers[index]]
   ...:     interchange = ForceField("2.1.0.offxml").create_interchange(
   ...:         molecule_copy.to_topology()
   ...:     )
   ...:     print(get_openmm_energies(interchange))
/Users/mattthompson/mambaforge/envs/openff-interchange-dev/lib/python3.11/site-packages/smirnoff99frosst/smirnoff99frosst.py:11: UserWarning: Module openff was already imported from None, but /Users/mattthompson/software/openff-interchange is being added to sys.path
  from pkg_resources import resource_filename
Energies:

Bond:               5.436333226473517 kilojoule / mole
Angle:              16.643336347971697 kilojoule / mole
Torsion:            17.72449198253925 kilojoule / mole
RBTorsion:          None
Nonbonded:          -10.911396912621075 kilojoule / mole
vdW:                None
Electrostatics:     None

Energies:

Bond:               5.436307019040268 kilojoule / mole
Angle:              16.64333666596585 kilojoule / mole
Torsion:            17.72449114558837 kilojoule / mole
RBTorsion:          None
Nonbonded:          -8.02872088810233 kilojoule / mole
vdW:                None
Electrostatics:     None

etc., etc.

There's also get_openmm_energies and a few others in case that's useful.

Guackkk commented 8 months ago

Thanks for your kind reply! It works well.