1.数据数据请参考前面文章:基于colab 分子仿真学习系列-allosteric binding 分子.
2.Colab ipynb数据链路: 包含分子接合运算+分子动力学
3.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
4.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...')
5.存为amber 可擦写轨迹档格式方便ambertools后处理:
simulation.reporters.append(NetCDFReporter('amber.nc', 1000)) # <-- AMBER compatible
6.下载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')