Combining Espaloma and Amber

[amin-sagar]

Dear Espaloma developers,
Thanks for this awesome work.
I am running some simulations of protein-peptide complexes where I want to treat protein with Amber and peptide with espaloma as it has some non canonical amino acids. With my current script, I am able to run the simulations but since this is the first time I am doing such simulations, I would like to be sure that I am doing this correctly. The following is my script. Does this seem like the correct way?

#import libraries
import os
import torch
import espaloma as esp
from openff.toolkit.topology import Molecule
from openmm.app import *
from openmm import *
from openmm.unit import *
from sys import stdout
from openmm.app import Modeller
from simtk import unit
from openmmforcefields.generators import EspalomaTemplateGenerator
from copy import copy
import shutil
import sys

#Load the molecules
##Load the peptide
molecule = Molecule.from_file("Pep.mol")
## create an Espaloma Graph object to represent the molecule of interest
molecule_graph = esp.Graph(molecule)
# load pretrained model
espaloma_model = esp.get_model("latest")
espaloma_model(molecule_graph.heterograph)
openmm_system = esp.graphs.deploy.openmm_system_from_graph(molecule_graph)
amber_forcefields = ['amber/protein.ff14SB.xml', 'amber14/tip3pfb.xml']
forcefield = ForceField(*amber_forcefields)
espaloma_generator = EspalomaTemplateGenerator(molecules=molecule, forcefield='espaloma-0.3.1')
openmm_topology = molecule.to_topology().to_openmm()
openmm_positions =  molecule.conformers[0].to_openmm()
forcefield.registerTemplateGenerator(espaloma_generator.generator)
#Load protein
protein = PDBFile('protein_fixed.pdb')
protein_modeller = Modeller(protein.topology, protein.positions)
protein_modeller.addHydrogens(forcefield)
protein_modeller.add(openmm_topology, openmm_positions)
#Add solvent
protein_modeller.addSolvent(forcefield, model='tip3p', padding=1*nanometer,boxShape='dodecahedron',ionicStrength=0.15*molar)
print('System has %d atoms' % protein_modeller.topology.getNumAtoms())

output_complex = 'Protein_Pep_solv_ions.pdb'
with open(output_complex, 'w') as outfile:
        PDBFile.writeFile(protein_modeller.topology, protein_modeller.positions, outfile)
system = forcefield.createSystem(protein_modeller.topology, nonbondedMethod=PME,nonbondedCutoff=0.9*nanometer, constraints=HBonds,hydrogenMass=1.5*amu)
#Set simulation parameters
temperature = 300 * unit.kelvin
collision_rate = 1.0 / unit.picoseconds
timestep = 4 * unit.femtoseconds
integrator = openmm.LangevinMiddleIntegrator(temperature, collision_rate, timestep)
platform = Platform.getPlatformByName('CUDA')

NVTsteps=50000
NPTsteps=1000000
num_steps=NVTsteps+NPTsteps

integrator = copy(integrator)
simulation = Simulation(protein_modeller.topology, system, integrator, platform)
simulation.context.setPositions(protein_modeller.positions)
simulation.minimizeEnergy(maxIterations=5000)
print('Saving...')
positions = simulation.context.getState(getPositions=True).getPositions()
pdbfilename = 'Protein_Bic_min_run.pdb'
PDBFile.writeFile(simulation.topology, positions, open(pdbfilename, 'w'))
print('Minimization Done')
dcdname = 'Protein_Bic_min_NVT_NPT.dcd'
simulation.reporters.append(DCDReporter(dcdname,5000))
logfilename = 'Protein_Bic_min_NVT_NPT.log'
simulation.reporters.append(StateDataReporter(logfilename, 5000, step=True, potentialEnergy=True, temperature=True, progress=True, remainingTime=True, speed=True,totalSteps=num_steps))
print('Running NVT\n')
simulation.step(NVTsteps)
system.addForce(MonteCarloBarostat(1*bar, 300*kelvin))
simulation.context.reinitialize(preserveState=True)
print("Running NPT")
simulation.step(NPTsteps)

[mikemhenry]

@amin-sagar good question, I will tag in @ijpulidos and @yuanqing-wang here since they may have the answer. It is a good question, but we have limited bandwidth to answer questions like this.

[yuanqing-wang]

Hi @amin-sagar , would it be possible if you could show us the minimal code to reproduce the error, which is not specific to your system?

[amin-sagar]

Thanks @yuanqing-wang and @mikemhenry .
Actually, there is no error. The simulations run perfectly fine. But because of my lack of experience in setting up simulations like this, I wanted to check with more experienced users/developers that this is the correct way.
I think a minimal code for setting up the system combining espaloma and amber is as follows.

#Load the molecules
##Load the peptide
molecule = Molecule.from_file("Pep.mol")
## create an Espaloma Graph object to represent the molecule of interest
molecule_graph = esp.Graph(molecule)
# load pretrained model
espaloma_model = esp.get_model("latest")
espaloma_model(molecule_graph.heterograph)
openmm_system = esp.graphs.deploy.openmm_system_from_graph(molecule_graph)
amber_forcefields = ['amber/protein.ff14SB.xml', 'amber14/tip3pfb.xml']
forcefield = ForceField(*amber_forcefields)
espaloma_generator = EspalomaTemplateGenerator(molecules=molecule, forcefield='espaloma-0.3.1')
openmm_topology = molecule.to_topology().to_openmm()
openmm_positions =  molecule.conformers[0].to_openmm()
forcefield.registerTemplateGenerator(espaloma_generator.generator)
#Load protein
protein = PDBFile('protein_fixed.pdb')
protein_modeller = Modeller(protein.topology, protein.positions)
protein_modeller.addHydrogens(forcefield)
protein_modeller.add(openmm_topology, openmm_positions)
#Add solvent
protein_modeller.addSolvent(forcefield, model='tip3p', padding=1*nanometer,boxShape='dodecahedron',ionicStrength=0.15*molar)

system = forcefield.createSystem(protein_modeller.topology, nonbondedMethod=PME,nonbondedCutoff=0.9*nanometer, constraints=HBonds,hydrogenMass=1.5*amu)

simulation = Simulation(protein_modeller.topology, system, integrator, platform)
simulation.context.setPositions(protein_modeller.positions)
simulation.minimizeEnergy(maxIterations=5000)

Best,
Amin

[jarvist]

Thanks for your minimum working example @amin-sagar - it helped us work around #194 :)

[amin-sagar]

Great! I am glad the example was useful. Since there are no new comments and the simulations seem to give results in agreement with the experiment, I am closing this issue.