-
Notifications
You must be signed in to change notification settings - Fork 2
/
prot_lig_openmm.py
172 lines (127 loc) · 6.17 KB
/
prot_lig_openmm.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
# This script is written by quantaosun@gmail.com or https://github.com/quantaosun/openmm.py for OpenMM simulation,2023.
from simtk.openmm import *
from simtk.openmm.app import *
from simtk.unit import *
from sys import stdout
from openmm.app import PDBFile # for final frame pdb generation
# Input Files
prmtop = AmberPrmtopFile('SYS_gaff2.prmtop')
inpcrd = AmberInpcrdFile('SYS_gaff2.crd')
# System Configuration
nonbondedMethod = PME
nonbondedCutoff = 1.0*nanometers
ewaldErrorTolerance = 0.0005
constraints = HBonds
rigidWater = True
constraintTolerance = 0.000001
# Integration Options
dt = 0.002*picoseconds
temperature = 300*kelvin
friction = 1.0/picosecond
pressure = 1.0*atmospheres
barostatInterval = 25
# Simulation Options
# long simulation (default)
#minimizationSteps = 10000
#productionSteps = 50000000 # 100 ns
#equilibrationSteps = 20000000 # 20 ns
##################################################
# short simulation
minimizationSteps = 100000
productionSteps = 10000000 # 10 ns
equilibrationSteps = 5000000 # 5 ns
##################################################
platform = Platform.getPlatformByName('CUDA')
platformProperties = {'Precision': 'single'}
dcdReporter = DCDReporter('prot_lig_prod.dcd', 10000)
minimizationDataReporter = StateDataReporter('stdout', 100, totalSteps=minimizationSteps,
step=True, time=True, speed=True, progress=True, remainingTime=True, potentialEnergy=True, temperature=True, separator='\t')
equilibrationDataReporter = StateDataReporter(stdout, 5000, totalSteps=equilibrationSteps,
step=True, time=True, speed=True, progress=True, remainingTime=True, potentialEnergy=True, temperature=True, separator='\t')
productionReporter = StateDataReporter(stdout, 5000, totalSteps=productionSteps,
step=True, time=True, speed=True, progress=True, remainingTime=True, temperature=True, separator='\t')
dataReporter = StateDataReporter('prot_lig_prod.log', 1000, totalSteps=productionSteps,
step=True, time=True, speed=True, progress=True, remainingTime=True, potentialEnergy=True, temperature=True, separator='\t')
# Prepare the Simulation
print('Starting molecular dynamics simulaiton......')
print('############################################################################')
print('## https://github.com/quantaosun/openmm.py ####')
print('############################################################################')
print('Initialising molecular dynamics....')
print('############################################################################')
print (' This script should only used in GPU platform, CPU not supported.')
print('############################################################################')
import time
# Sleep for 10 seconds
time.sleep(20)
print('############################################################################')
print('Building system...')
topology = prmtop.topology
positions = inpcrd.positions
system = prmtop.createSystem(nonbondedMethod=nonbondedMethod, nonbondedCutoff=nonbondedCutoff,
constraints=constraints, rigidWater=rigidWater, ewaldErrorTolerance=ewaldErrorTolerance)
system.addForce(MonteCarloBarostat(pressure, temperature, barostatInterval))
integrator = LangevinIntegrator(temperature, friction, dt)
integrator.setConstraintTolerance(constraintTolerance)
simulation = Simulation(topology, system, integrator, platform, platformProperties)
simulation.context.setPositions(positions)
if inpcrd.boxVectors is not None:
simulation.context.setPeriodicBoxVectors(*inpcrd.boxVectors)
#########################################################
content = """
###################################################
If you wish to change the simulation time, please modify
these three parameters inside the .py file
#############################################################
Default setting: 10000 steps of minimisation, followed by 10 ns
of equilibration, 50 ns of production.
##############################################################
Results include
###############################################################
########## prot_lig_mini.pdb ##############
########## prot_lig_equil.pdb ##############
########## prot_lig_prod.pdb ##############
########### prot_lig_prod.dcd ##############
##################################################################
# # # # ###### # # # #
# # ## # # # # # #
# # # # # ###### # # # #
# # # ## # # #
##### # # ####### # #
This script was written at The University of New South Wales, Sydney.
"""
print(content)
# Minimize and Equilibrate
print('Performing energy minimization...')
#simulation.minimizeEnergy()
simulation.minimizeEnergy(maxIterations=minimizationSteps)
print('Saving the last minimised frame as PDB file...')
positions = simulation.context.getState(getPositions=True).getPositions()
PDBFile.writeFile(topology, positions, open('prot_lig_mini.pdb', 'w'))
simulation.reporters.clear()
time.sleep(10)
print('##############################################')
print('Starting Equilibrating...')
simulation.context.setVelocitiesToTemperature(temperature)
simulation.reporters.append(equilibrationDataReporter)
simulation.step(equilibrationSteps)
print('Saving the last equilibrated frame as PDB file...')
positions = simulation.context.getState(getPositions=True).getPositions()
PDBFile.writeFile(topology, positions, open('prot_lig_equil.pdb', 'w'))
simulation.reporters.clear()
time.sleep(10)
# Production Simulate
print('##############################################')
print('Production Simulating...')
simulation.reporters.append(productionReporter) # on screen display.
simulation.reporters.append(dcdReporter)
simulation.reporters.append(dataReporter)
simulation.currentStep = 0
simulation.step(productionSteps)
print('##############################################')
# Write the last frame as a PDB file
print('Saving the last frame as PDB file...')
positions = simulation.context.getState(getPositions=True).getPositions()
PDBFile.writeFile(topology, positions, open('prod_lig_prod.pdb', 'w'))
time.sleep(5)
print('Congratulations, your simulation has finished!')