[gmx-users] charmm and oplsaa simulation results don't match

[Chris Neale]

If a user could be assured that they would obtained the same microstates, even in an infinite simulation, then we could 
finally get rid of the multitude of ff's and pick 'one' since they would all agree. However, this is not the case. Therefore
your assumption that you *should* get the same answer with different ff's is hopefully possible in the future, but currently false.
Second, are you even sure that you are converged? Try cutting the 100ns in 3 pieces and treating them as repeats. You will then be
able to estimate how much the difference depends on the different ff and how much depends on chance.

Chris.

-- original message --

Hi everybody,
                    I had run a long 100ns simulation on a tripeptide using opls/aa force field and 
spc water model.
I got a few microstates.
Then I tried to run simulation on the same system using charmm force field and tip3p water 
model.
However I couldn't find the microstates. Within the first 10 ns the system adopts the final
 configuration that was observed in the simulation using oplsaa force field.
 
 
 The process I used to run gromacs with charmm and tip3p water model is
 
 pdb2gmx -f opls.gro -ff charmm -ter -water tip3p
 
 charmm27-run-addtop-gromacs-dihe.sh topol.top new.top
  
  
 mv new.top topol.top
 
 
 #### I did modify the top file to include "tip3p-charmm.itp" 
 for water topology.
 
 
 editconf -f conf.gro -o box.gro -bt cubic -d 0.75
 
 genbox -cs spc216.gro -cp box.gro -p topol.top -o solvated.gro
 
 grompp -zero -f em.mdp -c solvated.gro -o em.tpr
 
 
 mdrun -v -s em.tpr
 
 
 Then I generated the "tpr" file for 100ns.
 
 The em.mdp file includes :
 
 
 
;       User spoel (236)
;       Wed Nov  3 17:12:44 1993
;       Input file
;
cpp                 = /usr/bin/cpp
constraints         =  none
integrator          =  steep
nsteps              =  100
;
;       Energy minimizing stuff
;
emtol               =  2000
emstep              =  0.01



nstcomm             =  1
ns_type             =  grid
rlist               =  1
rcoulomb            =  1.0
rvdw                =  1.0
Tcoupl              =  no
Pcoupl              =  no
gen_vel             =  no


The 100ns.mdp file includes:


;
;       User spoel (236)
;       Wed Nov  3 17:12:44 1993
;       Input file
;
title               =  Yo
cpp                 =  /usr/bin/cpp
constraints         =  all-bonds
integrator          =  md
dt                  =  0.002    ; ps !
nsteps              =  50000000 ; total 10000 ps.
nstcomm             =  1
nstxout             =  250
nstvout             =  1000
nstfout             =  0
nstlog              =  250
nstenergy           =  250
nstlist             =  10
ns_type             =  grid
rlist               =  1.0
rcoulomb            =  1.0
rvdw                =  1.0
; Berendsen temperature coupling is on in two groups
Tcoupl              =  berendsen
tc-grps             =  Protein  SOL
tau_t               =  0.1      0.1
ref_t               =  300      300
; Energy monitoring
energygrps          =  Protein  SOL
; Isotropic pressure coupling is now on
Pcoupl              =  berendsen
Pcoupltype          = isotropic
tau_p               =  0.5
compressibility     =  4.5e-5
ref_p               =  1.0
; Generate velocites is off at 300 K.
gen_vel             =  no
gen_temp            =  300.0
gen_seed            =  173529


Can anyone suggest why am I not getting the microstates with "charmm " force field.
Is there anything more that I needed to do and I missed. I am sorry for such a lengthy query


Thanks in advance

Sarbani