-
Notifications
You must be signed in to change notification settings - Fork 3
Expand file tree
/
Copy path5_gphmd.tex
More file actions
379 lines (353 loc) · 22.3 KB
/
Copy path5_gphmd.tex
File metadata and controls
379 lines (353 loc) · 22.3 KB
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
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
\section{Workflow for implicit-solvent CpHMD in Amber}
Here we discuss how to run GBNeck2-CpHMD \cite{Huang_Shen_2018_J.Chem.Inf.Model.,Harris_Shen_2019_J.Chem.Inf.Model.} simulations in Amber.
The related files can be found in the GitLab directory
\href{https://gitlab.com/shenlab-amber-cphmd/cphmd-tutorial/-/tree/main/gphmd_Amber}{cphmd-tutorial/gphmd\_Amber}.
\subsection{Preparation of structure and input files}
\begin{checklist}{Files generated by cphmd\_prep.sh}
\begin{itemize}
\item Cleaned up PDB file (.pdb)
\item Coordinates and topology file (.rst7)
\item CpHMD job control file (.phmdin)
\item CpHMD parameter file (.parm)
\item Minimization input file (\_mini.mdin)
\item Equilibration input files (\_equil*.mdin)
\item Template production file (template*prod.mdin)
\item Instructions and commands to run CpHMD (README)
\item Python script to run asynchronous pH replica exchange CpHMD if 'async' is specified (apHrex.py)
\item Amber group file to run pH replica-exchange CpHMD if 'reex' is specified (cphmd.groupfile)
\end{itemize}
\end{checklist}
To simplify system preparation for GBNeck2-CpHMD simulations,
we recommend the user to download and install a tool set \href{https://gitlab.com/shenlab-amber-cphmd/cphmd-prep/}{CpHMD-prep}.
Once installed, the shell script \href{https://gitlab.com/shenlab-amber-cphmd/cphmd-prep/-/tree/master/cphmd_tools/cphmd_prep.sh}{cphmd\_prep.sh} can be called to prepare the structure and input files for running CpHMD.
Before explaining \href{https://gitlab.com/shenlab-amber-cphmd/cphmd-prep/-/tree/master/cphmd_tools/cphmd_prep.sh}{cphmd\_prep.sh}, we briefly explain the important CpHMD specific files and options in Amber.
Similar to CpHMD simulations in CHARMM (section 3.1), we first convert the PDB file to the Amber format, cap the terminal groups, and change
the names of titratable Asp, Glu, and His residues to AS2, GL2, and HIP, respectively.
We then add dummy hydrogens to AS2 and GL2.
Besides the standard Amber force field files, tleap requires two CpHMD specific files for building topology and parameters:
\href{https://gitlab.com/shenlab-amber-cphmd/cphmd-prep/-/blob/master/Files/frcmod.phmd}{frcmod.phmd},
which specifies the modifications of bonded parameters for AS2 and GL2,
and \href{https://gitlab.com/shenlab-amber-cphmd/cphmd-prep/-/blob/master/Files/phmd.lib}{phmd.lib}, which
contains the definitions of AS2 and GL2.
For GBNeck2-CpHMD, set the PBradii to mbondi3 in tleap (see later discussion of modifications).
For running CpHMD, an additional file
\href{https://gitlab.com/shenlab-amber-cphmd/cphmd-prep/-/blob/master/Files/gbneck2_input.parm}{gbneck2$\_$input.parm} is needed, which contains the
partial charges and van der Waals energy flags of the protonated/deprotonated forms of the titratable residues as well as the parameters of the model titration PMFs (see section 7).
To run CpHMD, the \textit{pmemd.cuda} command should contain several CpHMD options.
The option
\textbf{-phmdparm} reads in
the aforementioned CpHMD parameter file \href{https://gitlab.com/shenlab-amber-cphmd/cphmd-prep/-/blob/master/Files/gbneck2_input.parm}{gbneck2\_input.parm}.
The option \textbf{-phmdin} reads in
a CpHMD job control file, which specifies various options for CpHMD simulations.
We note, except for
\textbf{MaskTitrRes(:)} (titratable residue types, e.g., 'AS2','GL2','HIP','CYS','LYS')
and \textbf{MaskTitrResTypes} (number of titratable residues, e.g., 5), the default settings in the phmdin file should be kept.
Finally, if specified the option \textbf{-phmdstrt} reads in a CpHMD restart file.
We should also point out that
for CpHMD runs, the Amber mdin files are
similar to those for the fixed-charge MD using GBNeck2 implicit solvent and Langevin dynamics.
The two CpHMD specific flags are \textbf{iphmd} and \textbf{solvph}.
For GBNeck2-CpHMD, \textbf{iphmd} should be set to 1 and \textbf{solvph} should be set to the desired pH condition.
Now we explain the shell script
\href{https://gitlab.com/shenlab-amber-cphmd/cphmd-prep/-/tree/master/cphmd_tools/cphmd_prep.sh}{cphmd\_prep.sh}.
\begin{lstlisting}
cphmd_prep.sh
[-pdb|--pdbid] | [-fil|--filename]
[-cha|--chainid]
[-mod|--modelid]
[-con|--conc]
[-tim|--time]
[-tem|--temp]
[-dph]
[-phl|--phlow]
[-phu|--phup]
[-res|--restype]
[-run|--runmode]
[-h|--help]
\end{lstlisting}
%[-tes|--test]
%[-ot|--outfiletype]
If specified, the option
\textbf{-pdb|--pdbid} reads in a PDB ID and downloads the corresponding structure from RCSB.
If no PDB ID is given, the option
\textbf{-fil|--filename} reads in a PDB file.
The option
\textbf{-cha|--chainid} specifies the chain ID to include in simulations (default 'A').
If specified, the option \textbf{-mod|--modelid} reads in a rotamer model ID (e.g., A or B) in the PDB file.
The following options are related to the simulation:
\textbf{-con|--conc} specifies the ionic strength in M (default 0.15);
\textbf{-tim|--time} specifies the simulation length in ns for each pH condition (default 10);
\textbf{-tem|--temperature} specifies the simulation temperature in K (default 300);
\textbf{-phl|--phlow} specifies the lowest pH (default 6.5);
\textbf{-phu|--phup} specifies the highest pH (default 9.0);
\textbf{-phi|--pHintvl} specifies the pH interval for calculating the individual
simulation pH conditions between phlow and phup;
\textbf{-res|--restype} specifies the amino acid residue types allowed to be titratable (default 'Asp Glu His Cys Lys'); and
\textbf{-rm|--runmode} specifies the type of simulation ('async' for asynchronous pH replica exchange,
'rex' for synchronous pH replica exchange,
and 'ind' for independent pH simulations).
Finally,
\textbf{-t|--test} offers a test option ('t' means temporary files will be kept for checking) and
\textbf{-h|--help} displays the usage and options.
\href{https://gitlab.com/shenlab-amber-cphmd/cphmd-prep/-/tree/master/cphmd_tools/cphmd_prep.sh}{cphmd\_prep.sh}
accepts either a PDB ID
or a user-prepared PDB file and calls PDBFixer in
OpenMM\cite{Eastman_Pande_2017_PLoSComput.Biol.}
to extract the desired chain, add missing heavy atoms
or residues and the terminal capping groups
(CH$_3$CO and NH$_2$).
Note, ions, ligands, waters,
and unspecified chains from the PDB file are removed.
% The option -mod|--modelid (default 'A') can be used to specify a particular sidechain rotamer model in the PDB file.
Once the PDB file is cleaned up, tleap in Amber (or Ambertools) is used to build the positions of hydrogen as well as dummy hydrogen atoms according to the specified titratable residue types.
As required by GBNeck2-CpHMD, \cite{Huang_Shen_2018_J.Chem.Inf.Model.,Harris_Shen_2019_J.Chem.Inf.Model.}
the GB input radii of HD1/HE2 of His, SG of Cys, OD1/OD2 of Asp, and OE1/OE2 of Glu are changed to 1.17, 2.00, 1.40, and 1.40 {\AA} respectively, and the interactions between the dummy carboxylate hydrogens in Asp/Glu are excluded.
\href{https://gitlab.com/shenlab-amber-cphmd/cphmd-prep/-/tree/master/cphmd_tools}{cphmd\_prep.sh}
automatically generates input files as well.
Below we give two examples.
\paragraph{Example 1. Predict the {\pka's} of the EGFR kinase given a PDB ID.}
Here we wish to predict
the {\pka's} and
the protonation states
of His, Cys, and Lys at the physiological pH 7.4.
%
\begin{lstlisting}
$ cphmd_prep.sh -pdb 5U8L -mod A -phl 5.5
-phu 9.5 -res 'His Cys Lys'
\end{lstlisting}
%
Here 5U8L is the PDB ID.
We omit chain ID, because 5U8L contains only one chain (default is chain A), but we specify the model A (there are two rotamer models in the PDB file).
For PDB file containing NMR models, a specific model can be selected using the -m flag.
The pH range is 5.5 to 9.5, which extends 2 pH
units above and below the physiological pH 7.4.
Since the (default) pH interval is 0.5 units, the total number of simulation pH conditions is 9.
In this example, simulation length (-l|--simlength), temperature (-tp|--temperature), and ionic strength (-i|--ionic) are set to the respective default values of 10 ns, 300 K, and 0.15 M.
Upon execution, cphmd\_prep.sh generates a file
`cphmd\_prep.log' and a folder
\href{https://gitlab.com/shenlab-amber-cphmd/cphmd-tutorial/-/tree/main/gphmd_Amber/gphmd_kinase_aphrex}{5U8L}, which
contains the files listed in the aforementioned checklist.
`cphmd\_prep.log' contains the information of the job control parameters and also the names of the files generated for the current job.
Inside the folder \href{https://gitlab.com/shenlab-amber-cphmd/cphmd-tutorial/-/tree/main/gphmd_Amber/gphmd_kinase_aphrex}{5U8L},
there are multiple files with names like 5U8L\_chainA\_A\_TAG.EXT with different tags TAG and extensions EXT.
'5U8L' is the PDB ID provided, 'chainA' refers to chain A, and 'A' is the rotamer model ID provided. Tags include mini for minimization,
equil* for different equilibration steps,
and prod for production.
Extensions include pdb for the PDB file,
rst7 for the Amber rst7 coordinate file,
parm for the Amber parm7 parameter file,
mdin for the Amber input file,
and phmdin for CpHMD control parameter file.
We can also specify multiple chains.
For example, if we want to include both chain A and chain B,
the command option -cha 'A B' can be used,
and the file names will start with 5U8L\_chainAB.
All files in the folder `5U8L' are also zipped as
\href{https://gitlab.com/shenlab-amber-cphmd/cphmd-tutorial/-/tree/main/gphmd_Amber/gphmd_kinase_aphrex}{cphmd\_inputs.zip}.
\paragraph{Example 2. Predict the {\pka's} of the protein Snase given a PDB file.}
%
\begin{lstlisting}
$ cphmd_prep.sh -fil 3BDC.pdb -mod A -phl 1.0
-phu 12.0 -tim 4
\end{lstlisting}
%
Here the script takes a PDB filename and the simulations are conducted in a wider pH range from 1.0 to 12.0.
With a default interval of 0.5 units, this
amounts to 23 simulation pH conditions.
Since the model {\pka's} of Asp/Glu are 3.7/4.3 and
the model {\pka} of Lys is 10.4, this pH range
typically allows the titration of
all default residue types 'Asp Glu His Cys Lys'.
For this reason, the residue type (-type) is omitted.
Another difference from Example 1 is that a much shorter simulation length of 4 ns (-tim 4) is specified.
This is because only a rough estimate of the {\pka} values is desired.
Similar files are generated as in Example 1; however,
there are 14 extra Amber input files due to the additional pH conditions.
All related files can be found here \href{https://gitlab.com/shenlab-amber-cphmd/cphmd-tutorial/-/tree/main/gphmd_Amber/gphmd_snase_aphrex/3BDC}{3BDC}.
Several additional files are generated by cphmd\_prep.sh.
\href{https://gitlab.com/shenlab-amber-cphmd/cphmd-tutorial/-/tree/main/gphmd_Amber/gphmd_kinase_aphrex/5U8L/apHREX.py}{apHREX.py}
is a Python script for running asynchronous pH replica exchange simulations if -run is not specified or specified as async.
The related file
\href{https://gitlab.com/shenlab-amber-cphmd/cphmd-tutorial/-/tree/main/gphmd_Amber/gphmd_snase_phrex/3BDC/cphmd.groupfile}{cphmd.groupfile}
is an Amber style for specifying individual replica files in pH replica-exchange simulations.
Finally, the file \href{https://gitlab.com/shenlab-amber-cphmd/cphmd-tutorial/-/tree/main/gphmd_Amber/gphmd_snase_phrex/3BDC/README}{README} contains all the commands for successfully running GBNeck2-CpHMD simulations for the current job.
\subsection{Running GBNeck2-CpHMD simulations}
Here we discuss the protocol of running GBNeck2-CpHMD simulations using the two examples discussed above.
All commands for minimization, equilibration, and production can also be found in the file
\href{
https://gitlab.com/shenlab-amber-cphmd/cphmd-tutorial/-/tree/main/gphmd_Amber/gphmd_snase_phrex/3BDC/README}{README}.
\paragraph{Energy minimization}
The energy minimization is performed
with the input file ending with `\_mini.mdin',
e.g.,
\href{https://gitlab.com/shenlab-amber-cphmd/cphmd-tutorial/-/blob/main/gphmd_Amber/gphmd_snase_phrex/3BDC/3BDC_chainA_A_mini.mdin}{3BDC\_chainA\_A\_mini.mdin} for SNase.
During minimization, a harmonic restraint with the
force constant of 50 kcal mol$^{-1}$ \AA$^{-2}$ is placed
on all protein heavy atoms, except for the
terminal and capping residues where all atoms are allowed to move.
The minimization lasts 1000 steps, using the steepest descent and
conjugate gradient algorithms for the first and second 500 steps, respectively.
To run minimization, we call \textbf{pmemd.MPI} or \textbf{pmemd.cuda} if Amber 20 \cite{Case_Kollman_2020} or a later version is used.
\begin{lstlisting}
$ export name='3BDC_chainA_A'
$ mpirun -n 4
$AMBERHOME/bin/pmemd.MPI -O
-i ${name}_mini.mdin -c ${name}.rst7
-p ${name}.parm7 -ref ${name}.rst7
-r ${name}_mini.rst7 -o ${name}_mini.out
\end{lstlisting}
% The latter requires at least an Nvidia GPU installed and configured properly as discusses above.
Here, '-n 4' specifies 4 CPU cores. Note, a large number of processors may be specified if hardware allows.
The `-ref' option
specifies the reference coordinate file
(3BDC\_ModelA\_A\_fixed.rst7).
The '-l' option specifies the initial coordinate file. It is omitted here, since it the same as the reference file.
If minimization runs error free, a timing info is given out
and two files are generated:
the file 3BDC\_ModelA\_A\_fixed\_mini.out
contains the printed energies, and
*\_mini.rst7 is the restart coordinate file for the next step equilibration.
\paragraph{Equilibration}
Following minimization, we perform a four-stage equilibration
to relax the protein structure at physiological pH 7.5 (or the crystallization pH).
Note, for implicit-solvent simulations heating is not required.
During the four-stage equilibration (2000 steps each),
the harmonic force constant on the heavy atoms is gradually reduced from 5.0, 2.0, 1.0, to 0.0 kcal/mol/\AA$^2$.
CpHMD is turned on by setting the iphmd keyword to 1
(GBNeck2-CpHMD)
in the Amber input file.
Using SNase as an example, the first equilibration stage is run with the following command:
\begin{lstlisting}
$ export name='3BDC_chainA_A'
$ $AMBERHOME/bin/pmemd.cuda -O
-i ${name}_equil1.mdin -c ${name}_mini.rst7
-p ${name}.parm7 -ref ${name}_mini.rst7
-r ${name}_equil1.rst7 -o ${name}_equil1.out
-x ${name}_equil1.nc -phmdin ${name}.phmdin
-phmdparm gbneck2_input.parm
-phmdout ${name}_equil1.lambda
-phmdrestrt ${name}_equil1.phmdrst
$ sed -i 's/QPHMDStart = .true.,/QPHMDStart
= .false.,/g' ${name}.phmdin
$ sed -i 's/PHMDRST/PHMDSTRT/g'
${name}.phmdrst
\end{lstlisting}
Here, 3BDC\_chainA\_A\_mini.rst7 generated by the minimization step is used as the initial (-c) and reference (-ref) coordinates.
The Amber job input file 3BDC\_chainA\_A\_equil1.mdin file specifies simulation length,
restraints, and other MD control parameters.
3BDC\_chainA\_A.phmdin contains job control parameters
for GBNeck2-CpHMD.
gbneck2\_input.parm, which is a common file
for GBNeck2-CpHMD simulations,
contains the model {\pka} values
and titration parameters.
3BDC\_chainA\_A\_equil1.rst7 is the restart file that contains coordinates and velocities.
3BDC\_chainA\_A\_equil1.phmdrst is the
corresponding CpHMD restart file that contains
the $\lambda$ values and velocities.
3BDC\_chainA\_A\_equil1.out contains the energy information.
3BDC\_chainA\_A\_equil1.lambda contains the $\lambda$ coordinates, although it is of little use since they are from the equilibration run.
The two 'sed' commands modify the names
for the next step.
For the rest of the equilibration stages,
the same commands are used
but file names should be changed. For example, we should replace 'equil1' and 'mini1' with 'equil2' and 'mini2' for the second equilibration step, and use 3BDC\_chainA\_A\_equil1.rst7 as the initial and reference coordinate file.
In addition, we need to use a new -phmdstrt flag and use 3BDC\_chainA\_A\_equil1.phmdrst as a CpHMD restart file.
After all four equilibration stages are completed,
files with the names *equil4* are generated.
The above example uses pmemd.cuda, but we can also use pmemd.MPI if no GPU is available. The only modification is to replace
`\$AMBERHOME/bin/pmemd.cuda' with `mpirun -n 16 \$AMBERHOME/bin/pmemd.MPI' in the command, assuming 16 cores can be used for the simulation.
\paragraph{Production}
For production runs, we can use either independent pH,
synchronous pH replica exchange, or asynchronous pH replica exchange.
%(\href{https://gitlab.com/shenlab-amber-cphmd/async_ph_replica_exchange}{Async-pHREX}).
For independent pH runs, the simulation at different pH
conditions are run independently either on the same or different machines.
The running order is arbitrary, and we collect all the output files together for data analysis.
The input files can be generated using
\href{https://gitlab.com/shenlab-amber-cphmd/cphmd-prep/-/tree/master/cphmd_tools/cphmd_prep.sh}{cphmd\_prep.sh}
with \textbf{-run ind}.
For SNase, the production run
is invoked with the following command:
\begin{lstlisting}
$ export name='3BDC_chainA_A'
$ $AMBERHOME/bin/pmemd.cuda -O
-i ${name}_pH7.5_prod.mdin
-c ${name}_equil4.rst7 -p ${name}.parm7
-r ${name}_pH7.5_prod1.rst7
-o ${name}_pH7.5_prod1.out
-x ${name}_pH7.5_prod1.nc
-phmdin ${name}.phmdin
-phmdparm gbneck2_input.parm
-phmdstrt ${name}_equil4.phmdrst
-phmdout ${name}_pH7.5_prod1.lambda
-phmdrestrt ${name}_pH7.5_prod1.phmdrst
-inf pH7.5.mdinfo
\end{lstlisting}
The restart files generated by equilibration step 4 are used as starting files for the production run. No reference file is needed.
Along with *.out, *.rst7, *.nc (trajectory),
and *.phmdrst files,
$\lambda$ coordinates are saved in 3BDC\_chainA\_A\_fixed\_pH*\_prod1.lambda.
For the command above, we again use pmemd.cuda
but we can also run on CPUs with pmemd.MPI.
All input files for this example can be found in \href{https://gitlab.com/shenlab-amber-cphmd/cphmd-tutorial/-/tree/main/gphmd_Amber/gphmd_snase_phind}{gphmd\_snase\_phind}.
pH replica-exchange simulations usually converge
much faster and thus require less simulation time than independent pH simulations \cite{Wallace_Shen_2011_J.Chem.TheoryComput.,Harris_Shen_2019_J.Chem.Inf.Model.,Harris_Shen_2020_J.Chem.TheoryComput.}.
The traditional replica-exchange algorithm
allows each replica to run in parallel and
data exchange is synchronized.
The pH replica-exchange protocol is executed
similarly as the temperature replica-exchange protocol in Amber,
i.e., we need a group file that contains the commands and inputs for every pH condition and then mpirun is called.
The commands in the group file are quite similar to those for
the independent pH runs.
The input files can be generated using
\textbf{cphmd\_prep.sh} with the option \textbf{-run rex}.
For example,
\begin{lstlisting}
mpirun -np 64
$AMBERHOME/bin/pmemd.MPI -ng 16
-groupfile cphmd.groupfile
\end{lstlisting}
Note that we use 4 processes for a single pH replica and a total of 64 processes are needed (-np 64) given the number of pH conditions (-ng 16) in the Snase example command.
For using mpirun for CpHMD, we should make sure that each pH replica runs on a single node so that there is only minimal data exchange between nodes if the whole job spreads across multiple node. Supposing we have a CPU cluster consisting of 64-core nodes, we should make sure the number of pH replicas is a factor of 64.
This example uses pmemd.MPI to run on CPU clusters.
If we have a GPU cluster, we can replace the '-np 64
\$AMBERHOME/bin/pmemd.MPI -ng 16' with '-np 1
\$AMBERHOME/bin/pmemd.cuda.MPI -ng 16'. However, if the number of GPUs is less than 16, it might not work properly.
The input files can be found in the \href{https://gitlab.com/shenlab-amber-cphmd/cphmd-tutorial/-/tree/main/gphmd_Amber/gphmd_snase_phrex}{gphmd\_snase\_phrex} folder.
The traditional replica-exchange algorithm was designed for CPU based HPC clusters, where a large number of processors are available such that all replicas can be run in parallel.
However, the latter condition cannot be always met for today's GPU workstation or even HPC cluster, i.e., the number of GPUs is smaller than the number of replicas. Thus, we implemented an asynchronous pH replica-exchange algorithm, \cite{Henderson_Shen_2020_J.Chem.Phys.}
in which each replica is run consecutively and data exchange is performed after all replicas are run.
The users can download the related Python program \href{https://gitlab.com/shenlab-amber-cphmd/async_ph_replica_exchange}{async\_ph\_replica\_exchange} that allows an arbitrary number of pH replicas to be run on an arbitrary number of GPUs.
A Python script apHREX.py generated by
\href{https://gitlab.com/shenlab-amber-cphmd/cphmd-prep/-/tree/master/cphmd_tools/cphmd_prep.sh}{cphmd\_prep.sh}
with `-run async'
is called to run the asynchronous pH replica-exchange simulations for 3BDC.
\begin{lstlisting}
$ export name='3BDC_chainA_A'
$ python apHREX.py 2000
template_${name}_prod.mdin ${name}.parm7
gbneck2_input.parm ${name}.phmdin
${name}_cphmd ${name}_equil4.rst7
1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5
5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5
9.0 9.5 10.0 10.5 11.0 11.5 12.0
\end{lstlisting}
Output files, including the standard Amber outputs, lambda files, and trajectory files, are concatenated on the fly according to the pH conditions instead of replicas in the pHREX method.
These files are stored in the 3BDC\_chainA\_A\_cphmd subdirectory.
The input files for this example can be found in \href{https://gitlab.com/shenlab-amber-cphmd/cphmd-tutorial/-/tree/main/gphmd_Amber/gphmd_snase_aphrex}{gphmd\_snase\_aphrex}.
\paragraph{Other relevant settings}
For optimum performance in {\pka} prediction,
the default settings in GBNeck2-CpHMD should be followed.
Currently, the Amber ff14SB force field
is recommended for proteins.
A Langevin thermostat with the collision frequency of 1 ps$^{-1}$ is used to maintain the temperature for $\lambda$ dynamics.
A non-bond cutoff of 999 {\AA} (i.e., no cutoff) is used.
Files should be saved infrequently to avoid IO that substantially slows down simulations on the GPUs.
The default frequency of saving $\lambda$ files
every 1000 MD steps
is designed for simulations of 10 ns per replica. For longer simulations, less frequent saving should be used.
If analysis of conformational dynamics is desired, a saving frequency of every 10,000 MD steps (20 ps with a 2-fs time step) is recommended
for a simulation of 10 ns per replica.
However, if only {\pka's} are desired, trajectory saving should be set to the frequency of restart file writing, e.g., 500,000 MD steps (1 ns with a 2-fs time step) or even larger depending on the wall clock time of the simulation (aiming at writing restart file once or twice a day). Writing to the Amber log file should be avoided or as infrequent as possible.