[gmx-users] Re: Re: Discrepancy between -chargegrp and -nochargegrp in simulations with CHARMM ff, Why ?

sa sagmx.mail at gmail.com
Sun Dec 5 18:22:04 CET 2010


Dear all,

I come back to you with my problem between -chargegrp and -nochargegrp with
new results which are contradict my previous
findings. Indeed, in my previous message i have suggested that the
discrepancies between -chargegrp and -nochargegrp in
my simulations of a 25 AA peptide in TIP3P water. I suspect that this
discrepancy came probably from an old version of the CHARMM27 ff and the
used in my simulations. To test this hypothesis i have done six additional
simulations of the same system (a 25 AA peptide in TIP3P water) with the new
released of the CHARMM27 ff and gromacs versions 4.5.1 and 4.5.3

I have done three simulation of with a *.top file generated by pdb2gromacs
v. 4.5.3 and modified by hand to change the charge
groups distribution in the peptide atoms according to the original CHARMM
force field and to mimic the -chargegrp option
available in the previous version of gromacs.
Three others MD were also done with *.top generated with pdb2gromacs of
gromacs 4.5.3 (i.e. with one charge group per atom).

1- gromacs4.5.1 and CHARMM27.ff -> gromacs 4.5.1*
2- gromacs4.5.3 and CHARMM27.ff -> gromacs 4.5.3*
3- gromacs4.5.3 and CHARMM27.ff and MD parameters used by Bjelkmar et al.
(J. Chem. Theory Comput. 2010, 6,
459–466-> gromacs 4.5.3**

non-bonded md parameters used in * simulations

>> 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

To summarize, as shown in the two figures accessible with the two links
below (hosted by hostingpics.net), all the peptides remain in their helix
conformation and no significant differences were observed with the md
parameters, the gromacs version and finally how  the charge group
distribution to compute the force during the simulations.

http://hpics.li/9ed5ee1  --> Secondary structure vs. time

http://hpics.li/7cf1625  --> RMSD vs. time

So i suspect that the errors come from the use a bad version of the CHARMM27
ff

Thank you all for your preceding comments

SA

Francesco Oteri wrote:
> > To see if the problem is force-field related, you could try to run the
> > same simulations using amber-ff.
> > If you will find the same results, probably is a software bug.
> >
>
> Some Amber parameter sets (Amber94, I think) have issues of being overly
> "helix-friendly," but perhaps other force fields, in general, might make
> for
> good comparisons.  But then, too, Gromos96 over-stabilizes extended
> conformations, so buyer beware...
>
> > Maybe the bug has been introduced in the version 4 when the Domain
> > Decomposition has been introduced.
> > You can check if it is a software problem using the 3.3.3 version.
> >
>
> I doubt that would work.  There have been many code changes in order to
> implement CHARMM.  If anything, try mdrun -pd to use the old particle
> decomposition mode, but I would be very hesitant to blindly say that DD
> might be
> causing this behavior.
>
> -Justin
>
> >
> >
> >
> > Il 25/11/2010 23:42, sa ha scritto:
> >> 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/20101205/5858c509/attachment.html>


More information about the gromacs.org_gmx-users mailing list