Tutorial 04 (Optimizing Force Fields)

[Shahriar-Keshvari]

Hello:

I'm following the tutorial number 4, and I have been able to run the program successfully with a customized forcefield. I did it for ethanol as described in the tutorial. However, I am interested in a different molecule. The program runs successfully with the new molecule and its customized forcefield but at some point it is crashing because of the kinetic energy is getting too large.

To solve the issue, I have run separately, openmm simulations and found that the problem comes from the time step and the damping coefficient of the Langevin coefficient and the attempt frequency of the MonteCarloBarostat (since I am interested in NPT ensemble).

I have two questions.
1- How can I set these parameters like time-step and the damping coefficient in tutorial 4?
2- How can I set the type of the integrator. Like If I am interested in MonteCarloBarostat or NoseHooverThermostat in tutorial 4?

Please let me know if you need more information
Thank
Shah

[lilyminium]

Hi Shah, glad you got the tutorial running with the custom molecule! I'll answer this in two parts, as the first question you asked is relatively easy and the second relatively hard.

Q1

To change the timestep (.timestep) and damping coefficient (.thermostat_friction), you need to modify the simulation schema generated in defining the estimation options.

from openff.evaluator.client import RequestOptions

# Create the options which evaluator should use.
evaluator_options = RequestOptions()
# Choose which calculation layers to make available.
evaluator_options.calculation_layers = ["SimulationLayer"]

# Reduce the default number of molecules
from openff.evaluator.properties import Density, EnthalpyOfVaporization
from openff.units import unit

density_schema = Density.default_simulation_schema(n_molecules=256)
h_vap_schema = EnthalpyOfVaporization.default_simulation_schema(n_molecules=256)

# now modify the equilibration and/or production simulation
# same code below applies to h_vap_schema
for protocol_schema in density_schema.workflow_schema.protocol_schemas:
    if "equilibration_simulation" in protocol_schema.id:
        protocol_schema.inputs[".timestep"] = 1 * unit.femtosecond
        protocol_schema.inputs[".thermostat_friction"] = 0.5 / unit.picosecond

    # this is for the production data collection simulation
    # conditional_groups are groups of protocols so a bit more complicated
    elif protocol_schema.id == "conditional_group":
        for protocol_name, conditional_protocol_schema in protocol_schema.protocol_schemas.items():
            if "production_simulation" in protocol_name:
                conditional_protocol_schema.inputs[".timestep"] = ...
                conditional_protocol_schema.inputs[".thermostat_friction"] = ...

evaluator_options.add_schema("SimulationLayer", "Density", density_schema)
evaluator_options.add_schema("SimulationLayer", "EnthalpyOfVaporization", h_vap_schema)

Q2

Unfortunately changing the type of integrator, barostat, etc. is relatively difficult as it's not natively supported in Evaluator.. You would need to do the below:

• fork Evaluator and work off your custom version
• modify this section of the code here, either by directly changing the OpenMMSimulation class or subclassing it: https://github.com/openforcefield/openff-evaluator/blob/98b50ba870d86c25f1fc9ba88319ca4d495bab9a/openff/evaluator/protocols/openmm.py#L644-L661
• If you subclass it, you will need to create a new simulation schema by following similar steps as defined in default_simulation_schema() (e.g. https://github.com/openforcefield/openff-evaluator/blob/98b50ba870d86c25f1fc9ba88319ca4d495bab9a/openff/evaluator/properties/properties.py#L420-L589 and https://github.com/openforcefield/openff-evaluator/blob/98b50ba870d86c25f1fc9ba88319ca4d495bab9a/openff/evaluator/properties/density.py#L186-L247