[gmx-users] Discrepancy between -chargegrp and -nochargegrp in simulations with CHARMM ff, Why ?
sa
sagmx.mail at gmail.com
Thu Nov 25 23:42:06 CET 2010
Dear All,
In a previous message (
http://lists.gromacs.org/pipermail/gmx-users/2010-November/055839.html), I
described the results obtained with MD performed with the CHARMM27 ff and
the chargegrp "yes" and "no" options of a peptide in TIP3P water simulated
gromacs. Since these results puzzled me a lot, i would like to share with
you others results obtained from the gromacs community advices to explain
these results.
In few words, the context of these simulations. One of my labmate did, 8
months ago (march/april), several simulations of a peptide (25 AA) with the
CHARMM27 ff (and CMAP). The peptide is a transmembrane segment (TM) and
belongs to a large membrane protein. This TM segment has an initial helical
conformation. The simulations were performed in a cubic box filled with app.
14000 TIP3P water (Jorgensen's model) with 2 Cl ions. To construct the
topology file of the system, -chargegrp "yes" with pdb2gmx and the MD were
done with the gromacs 4.0.5. For some reasons, he had to left the lab, and
my boss asked me to continue his work. When I checked their results, i was
very intrigued by these MD results because he found that the peptide keep
along all the simulation time (100 ns) its initial helical conformation.
This results are not in agreement with circular dichroism experiments which
are shown that the same peptide in water has no helix segment and is
completely unfold. I am aware that the simulation time is short compared to
experiment time scale, however since i haven't seen any unfolding events in
this simulation, so I was not very confident about these results.
To explain this inconsistency, I have suspected that the error came probably
of the use of the default -chargegrp with CHARMM ff in these simulations
since i have read several recent threads about the charge groups problems in
the CHARMM ff implementation in gromacs. To examine this hypothesis I have
done two simulations with last gromacs version (4.5.3) and two top files
containing charge groups and no charge groups for the peptide residus. I
used the *same* initial pdb file, box size and simulations parameters. The
two simulations were carried out during 24 ns in the NPT ensemble with the
md.mdp parameters described below after energy minimisation, NVT and NPT
equilibration steps.
constraints = all-bonds
integrator = md
nsteps = 12000000 ; 24000ps ou 24ns
dt = 0.002
nstlist = 10
nstcalcenergy = 10
nstcomm = 10
continuation = no ; Restarting after NPT
vdw-type = cut-off
rvdw = 1.0
rlist = 0.9
coulombtype = PME
rcoulomb = 0.9
fourierspacing = 0.12
fourier_nx = 0
fourier_ny = 0
fourier_nz = 0
pme_order = 4
ewald_rtol = 1e-05
optimize_fft = yes
nstvout = 50000
nstxout = 50000
nstenergy = 20000
nstlog = 5000 ; update log file every 10 ps
nstxtcout = 1000 ; frequency to write coordinates to xtc
trajectory every 2 ps
Tcoupl = nose-hoover
tc-grps = Protein Non-Protein
tau-t = 0.4 0.4
ref-t = 298 298
; Pressure coupling is on
Pcoupl = Parrinello-Rahman
pcoupltype = isotropic
tau_p = 3.0
compressibility = 4.5e-5
ref_p = 1.0135
gen_vel = no
I found that with charge groups, the peptide remains in its initial helical
conformation, whereas with no charge group, the peptide unfolds quickly and
has a random coil conformation. I have shown these results to my boss but I
was not able to explain why we observe these differences between the two
simulations. Indeed since i use PME in the MD, chargegroup should not affect
the dynamic results (correct ?) . He asked to do others simulations with
different versions of gromacs to see if is not a bug with charge group
implementation in gromacs. For testing i have done four others MD wit the
*same* initial pdb file, box size and simulations parameters and with
previous GMX version 4.5.0 and pre 4.5.2 and with the two types of top
files.
In addition since i have done simulations with non optimal parameter for
the vdW et electrostatic, i have also tested the influence of the vdW et
electrostatic MD parameters on the MD by performing two additional
simulations with the parameters used in the Bjelkmar et al . paper (J.
Chem. Theory Comput. 2010, 6, 459–466) (labeled GMX4.5.3 CHARMM in the
figures)
; Run parameters
integrator = md ; leap-frog integrator
nsteps = 12000000 ; 2 * 1000000 = 20000 ps
dt = 0.002 ; 2 fs
nstcomm = 10
nstcalcenergy = 10
; Output control
nstxtcout = 1000 ; frequency to write coordinates to xtc
trajectory every 10 ps
nstxout = 50000
nstvout = 50000 ; save velocities every 0.2 ps
nstenergy = 25000 ; save energies every 0.2 ps
nstlog = 5000 ; update log file every 0.2 ps
energygrps = protein Non-Protein
; Bond parameters
continuation = yes ; continuation after NVT
constraint_algorithm = lincs ; holonomic constraints
constraints = all-bonds ; all bonds (even heavy atom-H bonds)
constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
; Neighborsearching
ns_type = grid ; search neighboring grid cels
nstlist = 10 ; 20 fs
rlist = 1.2
;
coulombtype = PME
rcoulomb-switch = 0
rcoulomb = 1.2
vdw-type = Switch
; cut-off lengths
rvdw-switch = 1.0
rvdw = 1.2
DispCorr = EnerPres
table-extension = 5
fourierspacing = 0.14
fourier_nx = 0
fourier_ny = 0
fourier_nz = 0
pme_order = 4
ewald_rtol = 1e-05
ewald_geometry = 3d
epsilon_surface = 0
optimize_fft = no
; Temperature coupling is on
tcoupl = nose-hoover ;
tc-grps = protein Non-Protein ; one coupling groups - more
accurate
tau_t = 0.4 0.4 ; time constant, in ps
ref_t = 297 297 ; reference temperature, one for each group,
in K
; Pressure coupling is off
pcoupl = Parrinello-Rahman ; pressure coupling in NPT
pcoupltype = isotropic
tau_p = 3.0
compressibility = 4.5e-5
ref_p = 1.0135
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
gen_vel = no
The results are presented in the three figures accessibles with the three
links below
http://img4.hostingpics.net/pics/586637ResultatsVSGMX1jpg.jpg ->
Conformation
http://img4.hostingpics.net/pics/392578ResultatsVSGMX2jpg.jpg -> RMSD vs.
time
http://img4.hostingpics.net/pics/271597ResultatsVSGMX3jpg.jpg -> Secondary
structure vs. time
As you can see in all runs with no charge groups option, the peptide quickly
unfolds and have a random coil conformation during the end of the
simulation times, whereas with charge group option, the peptide remains in
its initial (helix) conformation whatever the gromacs version used.
To finish I saw that in all the simulation seems to be well conserved
CHARGE GROUP 4.5.3 (along the 24 ns of the run)
Energy Average Err.Est. RMSD Tot-Drift
-------------------------------------------------------------------------------
Potential -561452 5.6 774.891 -1.37923
(kJ/mol)
Total Energy -455901 5.6 959.234 -1.42286
(kJ/mol)
NO_CHARGE GROUP 4.5.3 (along the 24 ns of the run)
Energy Average Err.Est. RMSD Tot-Drift
-------------------------------------------------------------------------------
Potential -563934 18 774.632 -113.02
(kJ/mol)
Total Energy -458347 18 960.371 -112.951
(kJ/mol)
CHARGE GROUP pre4.5.2 (along the 24 ns of the run)
Energy Average Err.Est. RMSD Tot-Drift
-------------------------------------------------------------------------------
Potential -562415 5.1 777.12 -33.8615
(kJ/mol)
Total Energy -456864 5.1 964.365 -33.7504
(kJ/mol)
NO_CHARGE_GROUP_PRE_4.5.2
Energy Average Err.Est. RMSD Tot-Drift
-------------------------------------------------------------------------------
Potential -564899 4.6 778.002 -3.00182
(kJ/mol)
Total Energy -459312 4.6 963.702 -2.58102
(kJ/mol)
NO_CHARGE_GROUP_PRE_4.5.0
Energy Average Err.Est. RMSD Tot-Drift
-------------------------------------------------------------------------------
Potential -563933 30 777.279 -202.641
(kJ/mol)
Total Energy -458346 29 960.748 -202.578
(kJ/mol)
CHARGE_GROUP_PRE_4.5.0
Energy Average Err.Est. RMSD Tot-Drift
-------------------------------------------------------------------------------
Potential -561442 8 773.826 -49.7596
(kJ/mol)
Total Energy -455891 8 958.546 -49.7573
(kJ/mol)
USE_CHARGE_GROUP_AND_GMX_4.5.3 AND CHARMM parameters settings
Energy Average Err.Est. RMSD Tot-Drift
-------------------------------------------------------------------------------
Potential -566048 5.4 765.184 -31.6659
(kJ/mol)
Total Energy -460851 5.4 950.934 -31.7008
(kJ/mol)
USE_NO_CHARGE_GROUP_AND_GMX_4.5.3 AND CHARMM parameters settings
Energy Average Err.Est. RMSD Tot-Drift
-------------------------------------------------------------------------------
Potential -568520 24 770.264 -165.498
(kJ/mol)
Total Energy -463287 24 955.755 -165.768
(kJ/mol)
I am particularly eager to obtain from you some comments and advices about
these findings. Thanks you so much for your help.
SA
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20101125/80dfe4e4/attachment.html>
More information about the gromacs.org_gmx-users
mailing list