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
319 stars 93 forks source link

Openff cutoff distances for Sage? #1936

Closed dkchalmers closed 2 months ago

dkchalmers commented 2 months ago

Hi All,

An openff question - I have been looking for this information but have not been able to pin it down.

What are the recommended cutoff values for Sage when running simulations with Gromacs?

I have this at the moment.

cutoff-scheme = Verlet nstlist = 20 rlist = 1.2 vdwtype = Cut-off vdw-modifier = Force-switch rvdw_switch = 1.0 rvdw = 1.2 coulombtype = PME rcoulomb = 1.2

Cheers

David

mattwthompson commented 2 months ago

Short answer: 0.9 nm cutoff for all vdW interactions, "switching" starting at 0.8 nm, and compatible vanilla settings for PME electrostatics.


Sage (including all current versions) has vdW and electrostatics treatments specified in these lines:

Everything is meant to follow the SMIRNOFF specification, which should be considered correct over anything I say if we disagree. (I'm not always correct and, well, it's the spec.)

Here's what our infrastructure thinks should be used for non-bonded settings for a simple GROMACS simulation:

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

In [2]: molecule = Molecule.from_smiles("CCO")

In [3]: molecule.generate_conformers()

In [4]: topology = molecule.to_topology()

In [5]: topology.box_vectors = Quantity([5, 5, 5], "nanometer")

In [6]: interchange = ForceField("openff-2.2.0.offxml").create_interchange(topology)

In [7]: from openff.interchange.components.mdconfig import MDConfig

In [8]: MDConfig.from_interchange(interchange).write_mdp_file("1936.mdp")

In [9]: !cat 1936.mdp

nsteps                   = 0
nstenergy                = 1000
continuation             = yes
cutoff-scheme            = verlet

DispCorr                 = Ener
pbc = xyz
constraints = h-bonds
coulombtype = PME
rcoulomb = 0.9
coulomb-modifier = None
fourier-spacing = 0.12
ewald-rtol = 1e-4
vdwtype = cutoff
rvdw = 0.9
vdw-modifier = Potential-switch
rvdw-switch = 0.8

Interchange 0.4.0 will have a Interchange.to_gromacs method that produces this file alongside the topology and coordinate files, but that's still in a beta release.

dkchalmers commented 2 months ago

Thanks Matt, that is really helpful.

I had come across a page referencing the Interchange.to_gromacs method - but it is clearly not present in the version of Interchange that I have currently.

Cheers

David

On 12 Sep 2024, at 11:23 PM, Matt Thompson @.***> wrote:

Sage (including all current versions) has vdW and electrostatics treatments specified in these lines:

vdW: https://github.com/openforcefield/openff-forcefields/blob/main/openforcefields/offxml/openff-2.2.0.offxml#L339 Electrostatics: https://github.com/openforcefield/openff-forcefields/blob/main/openforcefields/offxml/openff-2.2.0.offxml#L379 Everything is meant to follow the SMIRNOFF specification https://openforcefield.github.io/standards/standards/smirnoff/, which should be considered correct over anything I say if we disagree. (I'm not always correct and, well, it's the spec.)

Here's what our infrastructure thinks should be used for non-bonded settings for a simple GROMACS simulation:

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

In [2]: molecule = Molecule.from_smiles("CCO")

In [3]: molecule.generate_conformers()

In [4]: topology = molecule.to_topology()

In [5]: topology.box_vectors = Quantity([5, 5, 5], "nanometer")

In [6]: interchange = ForceField("openff-2.2.0.offxml").create_interchange(topology)

In [7]: from openff.interchange.components.mdconfig import MDConfig

In [8]: MDConfig.from_interchange(interchange).write_mdp_file("1936.mdp")

In [9]: !cat 1936.mdp

nsteps = 0 nstenergy = 1000 continuation = yes cutoff-scheme = verlet

DispCorr = Ener pbc = xyz constraints = h-bonds coulombtype = PME rcoulomb = 0.9 coulomb-modifier = None fourier-spacing = 0.12 ewald-rtol = 1e-4 vdwtype = cutoff rvdw = 0.9 vdw-modifier = Potential-switch rvdw-switch = 0.8 Interchange 0.4.0 will have a Interchange.to_gromacs method that produces this file alongside the topology and coordinate files, but that's still in a beta release.

Sage was trained with a 0.9 nanometer cutoff for all condensed-phase fitting procedures (no cutoff for gas-phase properties, which GROMACS doesn't really support) so using a longer cutoff will be more computationally expensive and likely worse physical properties.

— Reply to this email directly, view it on GitHub https://github.com/openforcefield/openff-toolkit/issues/1936#issuecomment-2346272417, or unsubscribe https://github.com/notifications/unsubscribe-auth/AKY7UKYBTWY65VQ4URQPA33ZWGITXAVCNFSM6AAAAABOCFWFNOVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDGNBWGI3TENBRG4. You are receiving this because you authored the thread.


David Chalmers PhD BSc(Hons) Associate Professor

Monash University Medicinal Chemistry. Monash Institute of Pharmaceutical Sciences (MIPS) Faculty of Pharmacy and Pharmaceutical Sciences

Room 119. Level 1, Manning Building, Parkville Campus 381 Royal Parade, Parkville, VIC 3052. Australia

T: +61 3 9903 9110 E: @.*** Twitter: @davidkchalmers

mattwthompson commented 2 months ago

Yep - keep an eye out for 0.4.0 being released, hope to have it out in a few weeks depending on feedback in the beta/RC period. Until then, you can get by with calling each of the methods it wraps (Interchange.to_top, Interchange.to_gro, and the lines above to get the .mdp file).

mattwthompson commented 2 months ago

I'm closing this as I think the original question has been answered - please feel free to re-open if there are more details that are not discoverable in the documentation