1,4 scaling factors (lj+coulomb) overwritten by forcefield.createSystem()

[candidechamp]

Hello,

I found a bug which might be a bit tricky to fix, and will likely not affect most users. Essentially I am trying to run protein/ligand simulations with the implicit solvent model GB99dms (https://github.com/greener-group/GB99dms) and the GAFF force field (will later try OpenFF).

Unlike all amber force fields, the 1,4 scaling factors (scee and scnb variables in the code) are set to slightly different values in GB99dms. Hence to use a small molecule force field, the values must be set to these new values (by modifying the <NonbondedForce> line in the GAFF xml file).

When preparing a system (see code snippet below), the values are reset back at one point to the standard amber scaling factors due to a conflict with the dihedral types scee/scnb. i.e. the force field created does not match the xml given as input.

Specifically, the forcefield.createSystem() ends up calling the GAFFTemplateGenerator's generate_residue_template() function. In this function, the line params = parmed.amber.AmberParameterSet.from_leaprc(leaprc) leads to the bug, as the dihedral types are given a value for scee and scnb which is inconsistent with the params.default_scee and params.default_scnb.

I have found an ugly workaround to get this to work (but sets this unusual scaling factor for all openmm simulations in this environment), but it might be worth fixing for use with force fields from different families (e.g., CHARMM which have different scaling factors).

I was using openmm 8.2.0 (and 8.3.0) as well as openmmforcefields 0.14.2

example code:

import os
import sys
from math import ceil

from openff.toolkit import Molecule
from openmm import Platform, Vec3, XmlSerializer
from openmm.app import CutoffPeriodic, ForceField, HBonds, Modeller, PDBFile
from openmm.unit import amu, kelvin, nanometer, picosecond
from openmmforcefields.generators import GAFFTemplateGenerator

platform = Platform.getPlatformByName('CUDA')

pdb_filepath = "minimization.pdb"
sdf_filepath = "RL2.sdf"
out_prefix   = "debug_"

# We want to have periodic boundary conditions !!!
nonbondedMethod = CutoffPeriodic


dt = 0.004 # ps
n_steps       = int(ceil(1e6 / dt)) # 1 μs
n_steps_equil = int(ceil(500 / dt)) # 500 ps
n_steps_save  = int(ceil(100 / dt)) # 100 ps
temperature = 300.0 * kelvin
kappa = 0.7 / nanometer
friction = 1 / picosecond
nonbondedCutoff = 2 * nanometer
restraint_fc = 500.0 #* unit.kilojoule_per_mole # kJ/mol

pdb = PDBFile(pdb_filepath)
forcefield = ForceField("/home/cchampion/work/simulations/a_synuclein/binding/GB99dms.xml")
traj_fp       = out_prefix + ".xtc"
checkpoint_fp = out_prefix + ".chk"

# Open the ligand:
molecule = Molecule.from_file(sdf_filepath)
gaff = GAFFTemplateGenerator(molecules=molecule, forcefield="gaff-2.11")

forcefield.registerTemplateGenerator(gaff.generator)

openmm_top = molecule.to_topology().to_openmm()
lig_positions = [Vec3(x[0], x[1], x[2]) for x in molecule.conformers[0].m_as("nanometer")] * nanometer

# Combine
modeller = Modeller(pdb.topology, pdb.positions) # pdb.positions are a list of Vec3 in nm.
modeller.add(openmm_top, lig_positions)

# Print this out
with open(f"{out_prefix}_init.pdb", "w") as f:
    PDBFile.writeFile(modeller.topology, modeller.positions, f, keepIds=True)

# Back to the original

system = forcefield.createSystem(
    modeller.topology,
    nonbondedMethod=nonbondedMethod,
    nonbondedCutoff=nonbondedCutoff,
    constraints=HBonds,
    hydrogenMass=2*amu,
    implicitSolventKappa=kappa,
)