[gmx-users] Discrepancy between -chargegrp and -nochargegrp in simulations with CHARMM ff, Why ?
Francesco Oteri
francesco.oteri at gmail.com
Fri Nov 26 11:13:24 CET 2010
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.
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.
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
More information about the gromacs.org_gmx-users
mailing list