基于colab 分子仿真学习系列-allosteric binding 分子接合＋分子动力学运算

１．数据数据请参考前面文章：基于colab 分子仿真学习系列-allosteric binding 分子．

２．Colab ipynb数据链路: 包含分子接合运算+分子动力学

３．ambertools: top and crd files自动形成python代码．

!antechamber -fi mol2 -fo prepi -i dock.mol2 -o dock.prepi -c bcc -j 4 -at gaff
!parmchk2 -i dock.prepi -o dock.frcmod -f prepi
##leaprc.water.tip3p
!rm -rf tleapin

with open("tleapin", "w") as output:
    output.write('source /usr/local/dat/leap/cmd/leaprc.protein.ff14SB'+'\n')
    output.write('source /usr/local/dat/leap/cmd/leaprc.gaff'+'\n')
    output.write('source /usr/local/dat/leap/cmd/leaprc.water.tip3p'+'\n')
    output.write('loadamberprep dock.prepi'+'\n')
    output.write('mod = loadAmberParams dock.frcmod'+'\n')
    output.write('set UNL head null'+'\n')
    output.write('set UNL tail null'+'\n')
    output.write('aa = loadpdb dock-out.pdb '+'\n')
    output.write('bb = loadpdb rec.pdb '+'\n')
    output.write('cc = combine {aa bb}'+'\n')
    output.write('saveAmberParm cc  z.top z.crd'+'\n')
    output.write('addions cc Cl- 0'+'\n')
    output.write('addions cc Na+ 0'+'\n')
    output.write('solvatebox cc TIP3PBOX 20.0'+'\n')
    output.write('saveAmberParm cc  zw.top zw.crd'+'\n')
    output.write('quit'+'\n')

!tleap -s -f tleapin
!ls -l

４．OPENMM cMD simulations (Explicit solvent models+CUDA)自动形成python代码：

from simtk.openmm import *
from simtk.openmm.app import *
from simtk.unit import *
from mdtraj.reporters import NetCDFReporter

# Input Files

prmtop = AmberPrmtopFile('zw.top')
inpcrd = AmberInpcrdFile('zw.crd')

# System Configuration

nonbondedMethod = PME
nonbondedCutoff = 1.0*nanometers
ewaldErrorTolerance = 0.0005
constraints = HBonds
rigidWater = True
constraintTolerance = 0.000001
hydrogenMass = 1.5*amu

# Integration Options

dt = 0.004*picoseconds
temperature = 300*kelvin
friction = 1.0/picosecond
pressure = 1.0*atmospheres
barostatInterval = 25

# Simulation Options

steps = 500000
equilibrationSteps = 10000
platform = Platform.getPlatformByName('CUDA')
platformProperties = {'Precision': 'mixed'}
dcdReporter = DCDReporter('trajectory.dcd', 1000)
dataReporter = StateDataReporter('log.txt', 1000, totalSteps=steps,
    step=True, time=True, speed=True, progress=True, potentialEnergy=True, temperature=True, density=True, separator='\t')
checkpointReporter = CheckpointReporter('checkpoint.chk', 1000)

# Prepare the Simulation

print('Building system...')
topology = prmtop.topology
positions = inpcrd.positions
system = prmtop.createSystem(nonbondedMethod=nonbondedMethod, nonbondedCutoff=nonbondedCutoff,
    constraints=constraints, rigidWater=rigidWater, ewaldErrorTolerance=ewaldErrorTolerance, hydrogenMass=hydrogenMass)
system.addForce(MonteCarloBarostat(pressure, temperature, barostatInterval))
integrator = LangevinMiddleIntegrator(temperature, friction, dt)
integrator.setConstraintTolerance(constraintTolerance)
simulation = Simulation(topology, system, integrator, platform)
simulation.context.setPositions(positions)
if inpcrd.boxVectors is not None:
    simulation.context.setPeriodicBoxVectors(*inpcrd.boxVectors)

# Write XML serialized objects

with open("system.xml", mode="w") as file:
    file.write(XmlSerializer.serialize(system))
with open("integrator.xml", mode="w") as file:
    file.write(XmlSerializer.serialize(integrator))

# Minimize and Equilibrate

print('Performing energy minimization...')
simulation.minimizeEnergy()
print('Equilibrating...')
simulation.context.setVelocitiesToTemperature(temperature)
simulation.step(equilibrationSteps)

# Simulate

print('Simulating...')
#simulation.reporters.append(
#        NetCDFReporter('amber.nc', 100)
#)
simulation.reporters.append(dcdReporter)
simulation.reporters.append(dataReporter)
simulation.reporters.append(checkpointReporter)
simulation.currentStep = 0
simulation.reporters.append(NetCDFReporter('amber.nc', 1000))   # <-- AMBER compatible
simulation.step(steps)

# Write file with final simulation state

simulation.saveState("final_state.xml")
print('Simulating end...')

５．存为amber 可擦写轨迹档格式方便ambertools后处理：

simulation.reporters.append(NetCDFReporter('amber.nc', 1000))   # <-- AMBER compatible

６．下载cMD所有计算档案并完成计算:

import os
with open("cpptrajdcd", "w") as output:
    output.write('parm zw.top'+'\n')
    output.write('trajin amber.nc'+'\n')
    output.write('strip :WAT'+'\n')
    output.write('strip :Cl-'+'\n')
    output.write('strip :Na+'+'\n')
    output.write('trajout amber.x nobox'+'\n')

cmd ='cpptraj -i cpptrajdcd'
os.system(cmd)

from google.colab import files
files.download('z.top')
files.download('z.crd')
files.download('amber.x')