Example 1. protein sh3 based on one bead model

To use this example on your local computer, you need to clone or download the folder example from the GitHub repository of openmicron.

# import required packages
import openmm as mm
from openmm import app
from openmm import unit
import numpy as np
import pandas as pd
import mdtraj as md
from openmicron.forcefield.aicgmodel import AICG2Model
from openmicron import utils

Firstly, initialize the topology from the PDB and PSF files (if the PDB file includes the CONNECT section, you do not need to load the PSF to complete the topology in the PDB). Use the ParserNinfo from utils to parse the native information file to obtain the force field parameters, and store them in the ParserNinfo class variable.

# Initialize
T = 300
tot_simu_steps = 5
report_period = 1000
friction = 1
timestep = 20

# load pdb and psf files
pdb = md.load('../input/sh3_clementigo.pdb')
psf = md.load_psf('../input/sh3_clementigo.psf')

# complement the bond in topology
bonds = psf._bonds
top = pdb.topology.to_openmm()
rdtop = utils.RedefineTopology()
rdtop.redefine_bond(top,bonds)

# load native information
ParserNinfo=utils.ParserNinfo()
ParserNinfo.get_ninfo('../input/sh3_clementigo.ninfo')

After initializing the topology and native information file, create an AICG model using the AICG2Model class, and then create an OpenMM system based on the model. The specific parameters include the initialized top, use_pbc to control whether nonbonded interactions are used, and nonbondedMethod to set the properties of nonbonded interactions, including whether they are truncated and whether periodic boundary conditions are used.

model = AICG2Model()
model.create_system(top,use_pbc=True,
                 box_a=100, box_b=100, box_c=100,
                 nonbondedMethod=app.CutoffPeriodic,
                 remove_cmmotion=False)

After creating the system, we use the append_ff_params function to input the parsed force field parameters from ParserNinfo. For non-bonded interaction calculations, besides excluding atom pairs with residue indices , native contacts are also excluded, controlled by the exclude_nat_con parameter of the get_exclusion function. Then, add the various forces in sequence. It is important to assign a different force group to each force for easier subsequent force field analysis.

# append native information to model
model.append_ff_params(ParserNinfo)
np.save('../output/native_contact.npy',model.protein_intra_contact[['a1','a2','sigma']].to_numpy())
# get exclusion for nonbonded interaction
model.get_exclusion(exclude_nat_con=True)
# add force to system
model.add_protein_bond(force_group=0)
model.add_protein_harmonic_angle(force_group=1)
model.add_protein_native_dihedral(force_group=2)
model.add_protein_native_pair(force_group=3)
model.add_excluded(force_group=4)

Use OpenMM to create a Langevin integrator, then select the computing platform, such as CPU or CUDA, and build the OpenMM simulation.

# create a integrator
integrator = mm.LangevinIntegrator(T*unit.kelvin,friction/unit.picosecond,timestep*unit.femtosecond)
init_coord = pdb.xyz[0,:,:] * unit.nanometer
# create a simulation
model.set_simulation(integrator, platform_name='CPU',init_coord=init_coord)
model.move_COM_to_box_center(use_pbc=False)
model.simulation.context.setVelocitiesToTemperature(T*unit.kelvin)
model.simulation.minimizeEnergy()

Call the add_reporters function inherited to simulationsystem to set the input and output. tot_simu_steps specifies the total number of simulation steps, report_period controls the interval at which information is recorded, output_traj_name sets the output file path and name, report_traj_format sets the output trajectory format, with options of “dcd” and “xtc”. report_traj and report_state_log respectively set whether to output the trajectory file and the log file. Finally, run the simulation.

model.add_reporters(tot_simu_steps, report_period,
                 output_traj_name='../output/sh3_clementigo',report_traj_format='dcd'
                 ,report_traj=True,report_state_log=True)
print('Running simulation')
model.simulation.step(tot_simu_steps)