CHARMM c39b2 phmd.doc
Continuous constant pH Molecular Dynamics (PHMD)
Questions and comments regarding PHMD should be directed to
-----------------------------------------------------------
Jana Khandogin (janakhan@scripps.edu)
Charles L. Brooks III (brooks@scripps.edu)
The Scripps Research Institute
References:
1. M.S. Lee, F. R. Salsbury, Jr., and C.L. Brooks III,
Proteins, 56, 738-752 (2004).
2. J. Khandogin and C.L. Brooks III,
Biophys. J., 89, 141-157 (2005).
3. J. Wallace and J. Shen,
J. Chem. Theory and Comput., vXX, 000-000 (2011)
* Menu:
* Description:: Description of the PHMD Commands
* Syntax:: Syntax of the PHMD Commands
* Function:: Purpose of each of the commands
* Format:: Format of parameter file and how to obtain it
* Examples:: Usage examples of the PHMD module
^_
File: PHMD -=- Node: Description
Up: Top -=- Previous: Top -=- Next: Syntax
This module allows one to perform molecular dynamics and simultaneous
titration of specific ionizable residues under specified pH condition.
Titration occurs through the use of lambda variable measuring the
protonation progress of each titrating group. However, only two physical
states exist, namely, lambda = 0 for protonated, and lambda = 1 for
deprotonated states.
The lambda variables, themselves, are functions of theta variables,
through lambda(i) = sin^2[theta(i)]. The thetas can freely propagate without
need for restrictions. When theta = 0, +/- 2n(PI), lambda = 0. When theta = pi
+/- 2n(PI), then lambda = 1. The sin^2 function also provides a natural double
well for quadratic energy functions of lambda.
Analogous to the lambda variables, the x variables measure the tautomer
interconversion progress. The current implementation accounts for two
tautomeric states for either protonated case, such as the carboxylate group, or
deprotonated case, such as the histidine group.
The idea behind titration is that each group has a free energy of
titration when it is an isolated amino acid. in solution In other words, in the
absence of protein, the single group in solvent should spend 50% of the time
protonated and the other 50% of the time deprotonated. To achieve this, a model
energy function has to be derived, which is the potential of mean force of the
model compound titration. In the case of single-site titration (non-tautomer),
the model PMF has a simple quadratic form. In the case of double-site
titration (tautomer), it is a bivariate polynomial (2-d model potential
function), quadratic in both lambda and x variables.
While the same model function can be used for each of the non-termini
groups in the system. We have found that different model functions have to be
used for each possible C- and N-terminus residues. The model compound PMF
parameters are specified in a parameter file, which also serves to select the
desirable titrating groups. Another way to choose or exclude groups from
titration is to use the SELEction keyword in the PHMD command. In the Format
section, the procedure for deriving a model PMF function is explained.
For a double-site titrating group, a new residue type with dummy hydrogens
on both titrating sites has to be defined in the CHARMM topology file. The only
change in the CHARMM parameter file is related to raising the rotation barrier
to the C-O bond to prevent the dummy protons from losing the ability to titrate
once it is rotated to the anti-position (see the example section).
PHMD can be performed with the image facility in CHARMM. In this case,
an image transformation file needs to be read in prior to calling PHMD.
PHMD can be run using explicit solvent conformational sampling, with GBSW or
GBMV to control protonation sampling (see additional documentation in gbsw.doc
& gbmv.doc).
Replica exchange with PHMD has been implemented in REPDSTR replica exchange.
pH or Temperature exchange can be run (see additional documentation in repdstr.doc).
[SYNTAX: PHMD commands]
[starting PHMD]
PHMD { PAR <int> WRIte <int> PH <real> NPRInt <int> MASS <real> PHFRQ <int> BETA <real>
BARR <real> BARTAU <real> TEMP <real> MA1 <real> MA2 <real> MA3 <real>
[THETa] [DERIv] } [SELE atom-selection END]
[test and manipulation commands for PHMD]
used for deriving model PMF parameters
PHTEst { NUM <int> SET <real> }
{ NUM <int> STEP <real> }
{ NUM <int> FORCE <real> POS <real> }
File: PHMD -=- Node: Function
Up: Top -=- Previous: Syntax -=- Next: Format
-----------------------------------------------------------
Parameters of PHMD command
-----------------------------------------------------------
PAR Unit number for PHMD parameter file (input) MUST specify.
WRITE Unit number for PHMD trajectory file (output) MUST specify.
PH Titration pH (default: 1.0)
NPRINT Frequency of writing to PHMD trajectory file (default: 100)
MASS Mass of fictitious theta variable (default: 10)
BARR Quadratic barrier height for each theta variable (default: 2.0)
PHFRQ Frequency of updating theta/lambda variables (default :1)
BETA Friction coefficient (1/ps) for Langevin dynamics of theta variables (default: 0.0)
If not changed from default dynamics is run under Nose' thermostat.
BARTAU Quadratic barrier height for each x variable (default 2.5)
TEMP Nose thermostat temperature for configuation of thetas (default: 298)
MA1,MA2, Masses in Nose-Hoover thermostat multiplied by MASS.
MA3 (defaults: 3,5,7)
LAM Print lambda values in trajectory file. (default: none)
DERI Print theta and dE/dtheta values (not lambdas) in trajectory file.
SELE Use the SELE keyword to manually specify desirable titratable groups
-----------------------------------------------------------
Parameters of PHTE command
-----------------------------------------------------------
NUM Specify titratable group #. Use list generated at beginning
of PHMD output for reference.
SET Set value of theta(NUM)
STEP Increment value of theta(NUM) by STEP
FORCE/POS Place harmonic constraint on theta(NUM) with force constant, FORCE,
at equilibrium position, POS.
File: PHMD -=- Node: Format
Up: Top -=- Previous: Function -=- Next: Examples
Format of the Parameter file and
How to Derive Parameters
The parameter file is a series of entries. Each entry has the format:
1) For single-site titrating groups, such as NTAsp:
(NAME) (EXPERIMENTAL PK_1/2) (A) (B) (BARR)
ATOMTYPE_1 PROT_CHARGE_1 UNPROT_CHARGE_1 [PROT_RAD_1 UNPROT_RAD_1]
ATOMTYPE_2 PROT_CHARGE_2 UNPROT_CHARGE_2 [PROT_RAD_2 UNPROT_RAD_2]
... ... ... ...
ATOMTYPE_N PROT_CHARGE_N UNPROT_CHARGE_N [PROT_RAD_N UNPROT_RAD_N]
NAME : titrating residue name. For C- and N-termini groups, the
name consists of CT (NT) and the terminus residue name,
e.g., CTASP.
A/B : parameters of the PMF function A * ( lambda - B) ^ 2
BARR : barrier for suppressing population of mixed states, or prolonging
residence time for the pure states (default 1.5)
CHARGE : obtained from the CHARMM topology file. Make sure the difference
in the protonated and deprotonated states is 1.
RAD : only needed for the titrating proton: 1.0 for the protonated
form; 0.0 for the deprotonated form
2) for double-site titrating groups, such as ASP, GLU or HIS:
The parameter block for the first tautomer titration has the same look as above:
NAME : One needs to specify parameters for each tautomeric form.
In this case, NAME contains a number, e.g. 1 or 2, which distinguishes
between the first vs. the second tautomer forms.
A/B : coefficients in the quadratic function for the pure tautomeric states.
CHARGE : the dummy atom is assigned with zero charge in both protonated and
deprotonated forms.
Make sure the titrating proton is assigned with VDW radius 1.0 and 0.0.
in the protonated and deprotonated forms, respectively.
The parameter block for the second tautomer titration contains additional
numbers in the first two lines:
(NAME) (EXPERIMENTAL PK_1/2) (A) (B) (BARR) (A10) (B10) (BARTAU)
R1 R2 R3 R4 R5 R6
A10/B10: coefficients in the quadratic function for tautomer interconversion
BARTAU : analog to BARR: barrier for the tautomer interconversion
R1-R6 : parameters for constructing the 2-d model potential function
(bivariate polynomial)
How to derive model parameter values for single-site titration
(using variant of Thermodynamic Integration):
1a) Prepare a coordinate file of the desirable amino acid
capped by ACE and CT3, or un-capped if the terminus
(CT or NT) is to be titrated.
2a) Specify barr = 0, mass=1.E30 and use DERI keyword in the PHMD
input. Use PHTEST command to specify the titrating residue
and its theta value.
Supply a parameter file with pH=exp. pKa and A and B=0.
3a) Run 1ns dynamics at different values of fixed theta between
0 to PI/2. For example, theta = 0.4,0.6,0.8,1.0,1.2, 1.4.
Put corresponding set residue number and
4a) Use trajectory output of derivatives to calculate average
dE/dtheta derivative at each fixed value of theta.
5a) To obtain parameters, A and B, fit the values of dE/dtheta
to the following function, which is dE(model)/dtheta:
2*A*sin(2*theta)*(sin(theta)^2-B)
6a) To verify parameters, run PHMD of the model compound with
parameters plugged into parameter file and check whether the
model system titrates 50% protonated at its experimental pK_1/2.
How to derive model parameter values for double-site titration:
1b) Prepare coordinate file of the model compound with both titrating
sites protonated. A new residue with both sites protonated has to
be defined.
2b) Similar to 2a) except that two "groups" need to be specified
following the command PHTEST. The theta value that follows the
first group corresponds to the titration coordinate lambda while the
theta (or thetax) value that follows the second group corresponds
to the tautomer interconversion coordinate x.
3b) Run 1ns dynamics at different combinations of theta and thetax values
as given in 3a). It is useful to include combinations corresponding to
the pure tautomeric states (thetax=0.0 and PI/2), and the protonated
state (theta=0.0) for carboxyl groups and the deprotonated state
(theta=PI/2) for histidine.
4b) same as in 4a)
5b) Determine A and B as in 5a). For histidine, at theta=PI/2, fit dE/dx to
A10*(x-B10)^2 to obtain A10 and B10.
For carbxyl groups:
Determine R1 R2 and R3 by fitting A(lambda) to R1 lambda^2 + R2 lambda + R3
R4 = 0.5
Determine R5 by fitting A(x) to C1 x^2 + C2 x + R5
Determine R6 by fitting B(x) to C1 x^2 + C2 x + R6
File: PHMD -=- Node: Examples
Up: Top -=- Previous: Format -=- Next: Top
Usage and topology examples
The examples below illustrate how to use PHMD.
See test/phmd.inp for more examples.
--------------------------------------
NOTES TO RUN PHMD
--------------------------------------
1) Parameter file must be specified.
2) Works with GBSW or GBMV
*** note : GBMV does not currently support images,
therefore care should be used when attempting
to use GBMV with hybrid solvent PHMD using
periodic boundary conditions
NOT for use with PME, RDIE.
Example 1
! construct a residue with dummy hydrogens for titration
set name = asp
read sequence @name 1
generate @name first ace last ct3 setup
patch aspp2 @name 1
autogen angles dihed
ic para all
ic seed 1 n 1 ca 1 c
ic build
hbuild
ic gene
ic fill
ic edit
dihe 1 cb 1 cg 1 od1 1 hd1 180.0
dihe 1 cb 1 cg 1 od2 1 hd2 180.0
end
coor init sele type hd2 .or. type hd1 end
ic build
(write out psf and pdb files)
Example 2
! Perform a simple PHMD titration simulation on ASP:
set name = Asp
set barr = 2.25
set bartau = 2.5
set ph = 4.0
set temp = 298.0
set phmdpar = phmd-asp.in
(read in asp_h.psf and asp_h.pdb)
(invoke gbsw)
open unit 23 read form name @phmdpar
open unit 25 write form name @{name}.ph-@{ph}.lambda
PHMD PAR 23 WRI 25 PH @ph NPRI 100 -
BARR @barr BARTAU @bartau TEMP @temp
(dynamics)
Example 3
! Same as above except using Langevin dynamics for theta,
using theta update frequency of 10, and running hybrid solvent phmd.
set name = Asp
set barr = 2.25
set bartau = 2.5
set ph = 4.0
set temp = 298.0
set phmdpar = phmd-asp.in
(read in asp_h_solv.psf and asp_h_solv.pdb)
(setup periodic boundary conditions)
(setup images)
(invoke gbsw or gbmv with keyword hybrid) ***see gbsw.doc and gbmv.doc
open unit 23 read form name @phmdpar
open unit 25 write form name @{name}.ph-@{ph}.lambda
PHMD PAR 23 WRI 25 PH @ph NPRI 100 BETA 5.0 PHFRQ 10 -
BARR @barr BARTAU @bartau TEMP @temp
(dynamics)
Example 4
! Derive model potential function parameters for NtAla
set name = Ntala
set barr = 0.0
set mass = 1.0E30
set ph = 7.5
set temp = 298.0
set phmdpar = phmd-ntala_blank.in
set theta =0.4
(read in ntala_h.psf and ntala_h.pdb)
(invoke gbsw)
open unit 23 read form name @phmdpar
open unit 25 write form name @{name}.ph-@{ph}.lambda
phmd par 23 wri 25 ph @ph npri 100 -
barr @barr temp @temp
phtest num 1 set @theta
(dynamics)
Example 5
! Derive model potential function parameters for Asp
set name = Asp
set barr = 0.0
set bartau = 0.0
set mass = 1.0E30
set ph = 4.0
set temp = 298.0
set phmdpar = phmd-asp_blank.in
set theta =0.4
set thetax = 0.4
(read in asp_h.psf and asp_h.pdb)
(invoke gbsw)
open unit 23 read form name @phmdpar
open unit 25 write form name @{name}.ph-@{ph}.lambda
phmd par 23 wri 25 ph @ph npri 100
barr @barr bartau @bartau temp @temp
phtest num 1 set @theta
phtest num 2 set @thetax
(dynamics)
Example 6
! Do some manipulations of the theta variables:
! Incr theta #1 by 0.1
PHTEST NUM 1 STEP 0.1
! Incr theta #5 to 1.5
PHTEST NUM 5 SET 1.5
! Place harmonic restraint on theta #3 with
! force constant 100.0 kcal/mol and
! equilibrium value 0.5
PHTEST NUM 3 FORCE 100.0 POS 0.5
-----------------------------------------------
Additional patches in the CHARMM topology file
-----------------------------------------------
PRES ASPP2 0.00 ! patch for use in PHMD, proton on od1
GROUP ! and od1 via acetic acid, use in a patch statement
! ANGLes DIHEdrals are given
ATOM CB CT2 -0.21 !
ATOM HB1 HA 0.09 ! HB1 OD1-HD1
ATOM HB2 HA 0.09 ! | /
ATOM CG CC 0.75 ! -CB--CG
ATOM OD1 OC -0.55 ! | \
ATOM OD2 OC -0.61 ! HB2 OD2-HD2
ATOM HD1 H 0.0 HD2!
ATOM HD2 H 0.44 HD1!
BOND OD1 HD1
BOND OD2 HD2
DONOR HD1 OD1
DONOR HD2 OD2
IC HD1 OD1 CG OD2 0.0000 0.0000 0.0000 0.0000 0.0000
IC HD2 OD2 CG OD1 0.0000 0.0000 0.0000 0.0000 0.0000
PRES GLUP2 0.00 ! patch for use in PHMD, proton on od1
GROUP ! and od1 via acetic acid, use in a patch statement
! follow with AUTOGEN
ATOM CG CT2 -0.21 !
ATOM HG1 HA 0.09 ! HG1 OE1-HE1
ATOM HG2 HA 0.09 ! | /
ATOM CD CC 0.75 ! -CG--CD
ATOM OE1 OC -0.55 ! | \
ATOM OE2 OC -0.61 ! HG2 OE2-HE2
ATOM HE1 H 0.0 HE2!
ATOM HE2 H 0.44 HE1!
BOND OE1 HE1
BOND OE2 HE2
DONOR HE1 OE1
DONOR HE2 OE2
IC HE1 OE1 CD OE2 0.0000 0.0000 0.0000 0.0000 0.0000
IC HE2 OE2 CD OE1 0.0000 0.0000 0.0000 0.0000 0.0000
PRES CTRP2 0.00 ! patch for protonated CTER, proton on ot2
GROUP ! use in a patch statement, use AUTOGEN, ignore charges
ATOM C CC 0.72 ! OT1-HC1
ATOM OT1 OC -0.55 ! /
ATOM OT2 OC -0.61 ! -C
ATOM HC1 H 0.00 HC2!\
ATOM HC2 H 0.44 HC1! OT2-HC2
BOND OT1 HC1
BOND OT2 HC2
DONOR HC1 OT1
DONOR HC2 OT2
IC HC1 OT1 C OT2 0.0000 0.0000 0.0000 0.0000 0.0000
IC HC2 OT2 C OT1 0.0000 0.0000 0.0000 0.0000 0.0000
------------------------------------------------------------------
Additional parameters and modification in the CHARMM parameter file
-------------------------------------------------------------------
! additional parameters for CTRP and ASPP2
BONDS
!
!V(bond) = Kb(b - b0)**2
!
!Kb: kcal/mole/A**2
!b0: A
!
!atom type Kb b0
!
OC H 545.000 0.9600 ! ALLOW ALC ARO
! copy of EMB 11/21/89 methanol vib fit
ANGLES
!
!V(angle) = Ktheta(Theta - Theta0)**2
!
!V(Urey-Bradley) = Kub(S - S0)**2
!
!Ktheta: kcal/mole/rad**2
!Theta0: degrees
!Kub: kcal/mole/A**2 (Urey-Bradley)
!S0: A
!
!atom types Ktheta Theta0 Kub S0
!
H OC CC 55.000 115.0000 ! ALLOW ALC ARO PEP POL
! copy ! adm jr. 5/02/91, acetic acid pure solvent
DIHEDRALS
!
!V(dihedral) = Kchi(1 + cos(n(chi) - delta))
!
!Kchi: kcal/mole
!n: multiplicity
!delta: degrees
!
!atom types Kchi n delta
!
X CD OH1 X 3.0000 2 180.00 ! ALLOW PEP POL ARO ALC MSL
! ! adm jr, 10/17/90, acetic acid C-Oh rotation barrier
! ! Kchi can be modified if needed
X CC OC X 3.0000 2 180.00 ! ALLOW PEP POL ARO ALC MSL
! for CTRP ! Kchi can be modified if needed
IMPROPER
!
!V(improper) = Kpsi(psi - psi0)**2
!
!Kpsi: kcal/mole/rad**2
!psi0: degrees
!note that the second column of numbers (0) is ignored
!
!atom types Kpsi psi0
!
!OB X X CD 100.0000 0 0.0000 ! ALLOW ALC ARO POL
! adm jr., 10/17/90, acetic acid vibrations
OH1 OB CT2 CD 100.0000 0 0.0000 ! ALLOW ALC ARO POL
! ASPP1
END
File: PHMD -=- Node: Output
Up: Top -=- Previous: Examples -=- Next: Top
Output format
The only output from PHMD is a file that contains lambda values at specified
trajectory time steps. Following is an example output for the titration of ASP
(from phmd_2.inp in test directory):
# ititr 1 2
# ires 1 1
# itauto 3 4
100 0.86 0.25
line 1: gives the numbering for the titrating groups (runs to the total number)
line 2: gives the titrating residue number as in the PDB file
line 3: gives the type of titrating group:
0 - single-site
1 - titration of histidine
2 - tautomer interconversion in histidine
3 - titration of carboxyl groups
4 - tautomer interconversion in carboxyl groups
This information can be used in collecting statistics of protonation
populations.
line 4: column 1: step number; column 2: lambda value; column 3: x value
When PHTEST and DERI commands are used, dU/dtheta is being output. Following is
an example output for ASP (from phmd_1.inp in the test directory) :
# ititr 1 2
# ires 1 1
# itauto 3 4
100 0.4000 5.0330 0.6000 7.5975
200 0.4000 5.1584 0.6000 6.7531
line 3: Two numbers are printed out for each lambda or x trajectories. The first
is the theta or thetax value and the second is dU/dtheta or dU/dthetax,
respectively.