[gmx-users] Discrepancy between -chargegrp and -nochargegrp in simulations with CHARMM ff, Why ?
Justin A. Lemkul
jalemkul at vt.edu
Sat Nov 27 04:44:15 CET 2010
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
>
--
========================================
Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
========================================
More information about the gromacs.org_gmx-users
mailing list