generated from intbio/Armeev_et_al_2021
-
Notifications
You must be signed in to change notification settings - Fork 0
/
MD_production_protocol.mdp
176 lines (155 loc) · 7.6 KB
/
MD_production_protocol.mdp
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
173
174
175
176
;MDP-file production run: 10 000 ns.
; this protocol does not generate velocities (gen-vel = no), so velocities should be provided to grompp, e.g. from .cpt file
; VARIOUS PREPROCESSING OPTIONS
include =
define = ; NO RESTRAINTS DEFINED HERE
; will trigger the inclusion of posre.itp into your topology, used for implementing position restraints.
; here we assume that in posres.itp file all force constants have been changed to POSRES_FC
; see this link for details http://mdsquad.wikia.com/wiki/Change_position_restraint_force_constant_in_MDP
;1000 kJ mol-1 nm-2 is the default constant
;to match AMBER we need 100 kcal/mol/A2, which is 418 kJ/mol/A2 - we set 500
; RUN CONTROL PARAMETERS
integrator = md ; Steepest descent integrator
nsteps = 5000000000 ; 5 000 000 000 = 10 microsecond
init-step = 0 ; For exact run continuation or redoing part of a run
simulation-part = 1 ; Part index is updated automatically on checkpointing (keeps files separate)
comm-mode = Linear ; mode for center of mass motion removal
nstcomm = 100 ; number of steps for center of mass motion removal
comm-grps = System ; group(s) for center of mass motion removal default is the whole system;
dt = 0.002 ; this is only for MD
; ENERGY MINIMIZATION OPTIONS
; Force tolerance and initial step-size
emtol = 100.0
emstep = 0.01
; OUTPUT CONTROL OPTIONS
; Output frequency for coords (x), velocities (v) and forces (f)
nstxout = 0
nstvout = 0
nstfout = 0
; Output frequency for energies to log file and energy file
nstlog = 1000
nstcalcenergy = 100 ; why less -> looks like this is needed for thermostats/barostats
nstenergy = 1000
; Output frequency and precision for .xtc file
nstxout-compressed = 500000 ; every 1 ns
compressed-x-precision = 1000
; This selects the subset of atoms for the compressed
; trajectory file. You can select multiple groups. By
; default, all atoms will be written.
compressed-x-grps =
; Selection of energy groups
energygrps =
; NEIGHBORSEARCHING PARAMETERS
; cut-off scheme (Verlet: particle based cut-offs, group: using charge groups)
cutoff-scheme = Verlet
; nblist update frequency
nstlist = 10
; ns algorithm (simple or grid)
ns-type = Grid
; Periodic boundary conditions: xyz, no, xy
pbc = xyz
periodic-molecules = no
; Allowed energy error due to the Verlet buffer in kJ/mol/ps per atom,
; a value of -1 means: use rlist
verlet-buffer-tolerance = 0.005
; nblist cut-off
rlist = 1.2 ; actually will be ignore and calculated from verlet-buffer-tolerance
; OPTIONS FOR ELECTROSTATICS AND VDW
; Method for doing electrostatics
coulombtype = PME
coulomb-modifier = None ; change from default Potential-shift-Verlet because Amber does not use shift.
rcoulomb-switch = 0
rcoulomb = 0.8
; Relative dielectric constant for the medium and the reaction field
epsilon-r = 1
epsilon-rf = 0
; Method for doing Van der Waals
vdwtype = Cut-off
vdw-modifier = None ; this matches Amber
; cut-off lengths
rvdw_switch = 0.8
rvdw = 0.8
; Apply long range dispersion corrections for Energy and Pressure
DispCorr = EnerPres ; This is default in AMBER, ??? AllEnerPres what is it?
; Extension of the potential lookup tables beyond the cut-off
table-extension = 1
; Separate tables between energy group pairs
energygrp-table =
; Spacing for the PME/PPPM FFT grid
fourierspacing = 0.12
; FFT grid size, when a value is 0 fourierspacing will be used
fourier-nx = 0
fourier-ny = 0
fourier-nz = 0
; EWALD/PME/PPPM parameters
pme-order = 4
ewald-rtol = 1e-05
ewald-geometry = 3d
epsilon-surface = 0
; OPTIONS FOR BONDS
constraints = all-bonds ; the same as used in AMBER
; Type of constraint algorithm
constraint_algorithm = LINCS ; this is better than SHAKE so we retain it
; Do not constrain the start configuration
continuation = no
; Highest order in the expansion of the constraint coupling matrix
lincs-order = 4
; Number of iterations in the final step of LINCS. 1 is fine for
; normal simulations, but use 2 to conserve energy in NVE runs.
; For energy minimization with constraints it should be 4 to 8.
lincs-iter = 1 ; in MD we will set this to 1
; Lincs will write a warning to the stderr if in one step a bond
; rotates over more degrees than
lincs-warnangle = 30
; Below is only relevant for MD
; GENERATE VELOCITIES FOR STARTUP RUN
gen-vel = no
gen-temp = 300
gen-seed = -1
; OPTIONS FOR WEAK COUPLING ALGORITHMS
; Temperature coupling
tcoupl = v-rescale
nsttcouple = -1 ; the frequency of coupling the temperarure, -1 is automatic = nstlist for md integrator
; Groups to couple separately
tc-grps = System ; NOTE: Group name "System" is defined if you are NOT using an index file. But better to generate a default index file anyway.
; Time constant (ps) and reference temperature (K)
tau-t = 1 ; 1 ps as used in Amber paper, but we might need to research and use different for large scale dynamics!
ref-t = 300
; pressure coupling
pcoupl = Parrinello-Rahman
pcoupltype = Isotropic
nstpcouple = -1
; Time constant (ps), compressibility (1/bar) and reference P (bar)
tau-p = 1
compressibility = 4.5e-5 ; for water
ref-p = 1 ; 1 bar
; Scaling of reference coordinates, No, All or COM
refcoord-scaling = all ; The reference coordinates are scaled with the scaling matrix of the pressure coupling.
; ??? check if com is better?
; ENERGY GROUP EXCLUSIONS
; Pairs of energy groups for which all non-bonded interactions are excluded
energygrp-excl =
; Pull code
pull = yes
pull_ncoords = 2 ; only one reaction coordinate
pull_ngroups = 4 ; two groups defining one reaction coordinate
pull_group1_name = a_11
pull_group2_name = a_9326
pull_group3_name = a_4653
pull_group4_name = a_4684
pull_coord1_type = umbrella ; harmonic potential
pull_coord1_geometry = distance ; simple distance increase
pull_coord1_dim = Y Y Y ; pull along z
pull_coord1_groups = 1 2 ; groups 1 (Chain A) and 2 (Chain B) define the reaction coordinate
pull_coord1_start = no ; define initial COM distance > 0
pull_coord1_rate = 0.00 ; 0.01 nm per ps = 10 nm per ns
pull_coord1_k = 1000 ; kJ mol^-1 nm^-2, 4 kT per A which is ok
pull-coord1-init = 0.8431168557081486
pull_coord2_type = umbrella ; harmonic potential
pull_coord2_geometry = distance ; simple distance increase
pull_coord2_dim = Y Y Y ; pull along z
pull_coord2_groups = 3 4 ; groups 1 (Chain A) and 2 (Chain B) define the reaction coordinate
pull_coord2_start = no ; define initial COM distance > 0
pull_coord2_rate = 0.00 ; 0.01 nm per ps = 10 nm per ns
pull_coord2_k = 1000 ; kJ mol^-1 nm^-2, 4 kT per A which is ok
pull-coord2-init = 0.8669820894546378