[gmx-users] Spurious results with pdb2gmx -chargegrp "yes" option

[sa]

Dear All,

I would like to get little feedback of my findings/experience with recent
simulations I did with the CHARMM force field in GROMACS and the
-chargegroup options with pdb2gmx. I have performed two simulations with a
helical peptide (25 AA) in a cubic box filled with explicit TIP3P water. I
used the GMX 4.5.3. The first simulation was performed with -chargegroup
option "yes" and the second with option "no" (default).

For the first simulation, I used the old version of CHARMM -> gromacs
forcefield conversion given with GMX 4.0.5. I have done this MD because i
would like to compare some "suspect" results obtained with the same peptide
performed 6 months ago with the GMX 4.0.5 version and the old
CHARMM->gromacs files available in this distribution. In case of the second
simulation, I used the most recent charmm27.ff given in GMX 4.5.3 with a
single atom charge groups.

The two simulations were performed during 24 ns with the same protocol (i.e.
with same mdp files, below and the same starting files).

Briefly the protocol i used to equilibrate my system :

EM (1000 steps with steep) -> NVT at 298 K (with berendsen thermostat and
restraints) (100 ps) -> (with nose hoover thermostat and Parinello-Rahman
barostat) 400 ps.

Below md.mdp for the two runs ( i am aware that the Coulomb and
electrostatic parameters in the mdp are not otimals for the CHARMM ff, but
for sake of comparison with others simulations, i didn't change them).

title           = KTM17-TIP3 MD
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 2 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 the use of the two -chargegrp options influence drastically the
MD results especially in my case for the structural properties of the
peptides. Indeed, for the first simulation, the peptide keep in its initial
helical with a RMSD around 2.0 A along all 24 ns, whereas in the second
simulation, the peptide adopt quickly (~10 ns) a unfold conformation with a
RMSD value around 5-6.0 A after only 10 ns. Circular dichroism experiments
have shown that this peptide adopt an unfold structure in water. So I am
sure that the first simulation is not correct (even if the mdp parameters
are not optimal) and that my results confirm that chargegrp "yes" option
should not be used in MD with CHARMM force field in GROMACS, as it was
discussed recently in this mailing list (for example
http://lists.gromacs.org/pipermail/gmx-users/2010-September/054106.html).

SA
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20101114/079ad4c4/attachment.html>