How to do pulling simulations

[aminsagar]

Dear OpenSMOG developers,

Thanks for this awesome work. I would like to do some pulling simulations with the SMOG model to simulate protein unfolding. Are there any examples/documentation about how to do this? I would be really grateful for any suggestions.

Best, Amin.

[Whitford]

Amin, Thank you for a nice comment. Yes, we have some examples floating around, but we haven't organized it for distribution (yet). I've tagged Esteban Dodero-Rojas on this thread, because I know he has done a lot of simulations with pulling forces. He may have a simple example to get you started.

[ed29rice]

Amin, Here is a simple example of how to introduce external forces into SBMs:

from OpenSMOG import SBM from openmm.app import from openmm import from openmm.unit import *

Initialize simulation as normal

sbm = SBM(name='sim', time_step=0.002, collision_rate=1.0, r_cutoff=0.65, temperature=1) sbm.setup_openmm(platform='cuda',GPUindex='default') sbm.saveFolder('output')

Input files

sbm_grofile = './sim.gro' sbm_topfile = './sim.top' sbm_xmlfile = './sim.xml'

Load input files into the sbm class

sbm.loadSystem(Grofile=sbm_grofile, Topfile=sbm_topfile, Xmlfile=sbm_xmlfile)

Define the external force using OpenMM's CustomExternalForce function

force = CustomExternalForce('-K*z') #Constant force force.addGlobalParameter('K', 1)

Add particles to the force

particle_indexes = [0,1,2,3] for index in particle_indexes: force.addParticle(int(index))

Add force to the SBM system

sbm.system.addForce(force)

Run simulation

sbm.createSimulation() sbm.minimize() sbm.run(nsteps=106, report=True, interval=103)

Notice that the force needs to be included before createSimulation(). Let us know if this answers your question.

[aminsagar]

Thanks @Whitford and @ed29rice

I have used your script to write the following script to restrain the other end of the protein.

from OpenSMOG import SBM
from openmm.app import *
from openmm import *
from openmm.unit import *

def load_positions(grofile):
    with open(grofile, 'r') as f:
        lines = f.readlines()
    positions = []
    for line in lines[2:-1]:  # Skip the first two lines and the last line
        parts = line.split()
        x = float(parts[3]) * nanometers
        y = float(parts[4]) * nanometers
        z = float(parts[5]) * nanometers
        positions.append((x, y, z))
    return positions

sbm_AA = SBM(name='EC1-pull', time_step=0.002, collision_rate=1.0, r_cutoff=0.65, temperature=0.5, cmm=False)
sbm_AA.setup_openmm(platform='cuda',GPUindex='default')
sbm_AA.saveFolder('output_EC1_pull_AA_10p8_3/')
sbm_AA_grofile = 'Rank68_orient_correct_noH_AA.gro'
sbm_AA_topfile = 'Rank68_orient_correct_noH_AA.top'
sbm_AA_xmlfile = 'Rank68_orient_correct_noH_AA.xml'
sbm_AA.loadSystem(Grofile=sbm_AA_grofile, Topfile=sbm_AA_topfile, Xmlfile=sbm_AA_xmlfile)

# Load positions
positions = load_positions(sbm_AA_grofile)

force = CustomExternalForce('-K*x') #Constant force
force.addGlobalParameter('K', -0.72*kilocalories_per_mole/angstrom)
#Add particles to the force
pulled_indexes = [4490]
for index in pulled_indexes:
    force.addParticle(int(index))
#Restrain one end
restraint = CustomExternalForce('k*((x-x0)^2+(y-y0)^2+(z-z0)^2)')
restraint.addGlobalParameter('k', 100*kilocalories_per_mole/nanometers**2)
restraint.addPerParticleParameter('x0')
restraint.addPerParticleParameter('y0')
restraint.addPerParticleParameter('z0')
fixed_indexes = [6]
for index in fixed_indexes:
    x0, y0, z0 = positions[index]
    restraint.addParticle(int(index), [x0, y0, z0])

#Add force to the SBM system
sbm_AA.system.addForce(force)
sbm_AA.system.addForce(restraint)

sbm_AA.createSimulation()
sbm_AA.createReporters(trajectory=True, energies=True, energy_components=True, interval=10**3)
sbm_AA.minimize()
sbm_AA.run(nsteps=(10**8)*5, report=True, interval=10**3)

Does this seem correct?

I have another question. What are the units of force here? Do I need to use reduced units? Here, I am trying to pull at 50 pN (50*0.0144) = 0.72 kcal/mol/Å. Is this the correct way of specifying the pulling force.

Best, Amin.

[Whitford]

Also, make sure to turn off center of mass motion removal, whenever you apply an external force.

[aminsagar]

Thanks. Sorry, I pasted the wrong script earlier. I have updated it now.

If it looks okay, can you please help me with the units of force.

[Whitford]

This is a good question. There is no perfect way of estimating the energy scale, but here is my preferred approach.

First, you need to determine how a reduced temperature and relates to temperature in K. As a general guide, with our models you can typically estimate that 0.5-0.6 will generate native-basin fluctuations that are consistent with those observed at ~300K. This comparison has been done for several systems in Jackson, Nguyen and Whitford, Int. J. Mol. Sci. 2015, 16(4), 6868-6889; https://doi.org/10.3390/ijms16046868

Next, you can convert the reduced energy scale using the ratio of kBT in the desired units and kBT in reduced units. So, if you using a constant force of the form F=Kx, where K is given in units of (reduced energy unit)/nm, then you can convert K to any desired units by multiplying by the ratio. For example, the force constant in pN would be K(4.114 pNnm)/(0.6). In this example, I'm using 0.6 as the reference reduced temperature.

Hope that helps.

[aminsagar]

Thanks @Whitford . This is very helpful.

Amin.