[gmx-users] incomparable results with files produced with gromacs3.2.1 and 3.3.1

[Sampo Karkola]

Dear list,

I am trying to simulate a crystal structure of a protein with a cofactor 
(NADP, I've used NDPP topology) and a (self-made) ligand in a box of 
water. The preparation of the simulation went normally: pdb2gmx, 
editconf, genbox, grompp (for ions), genion, grompp (for em) mdrun (em), 
grompp (for md), mdrun (md). Now, in the md simulation, the potential 
and total energies of the system rise immediately from the beginning, 
stabilize to a certain level and later the run crashed with LINCS 
warnings. Earlier I have succesfully simulated a model of an another 
enzyme with similar procedure and in that run the potential energy 
decreased immediately and stabilized with no crashes. I tried to 
simulate the protein alone, the cofactor alone and the ligand alone and 
the energies rise in all cases.

Luckily, I had an older simulation (with gromacs 3.2.1, now I'm using 
3.3.1) with the same ligand (except one atom type is changed from NR in 
the old topology to N in new topology) alone and there the energies 
decreased as they should. To study the differences in the ligand-alone 
simulations, I performed an mdrun with the old tpr file (mdrun -s 
old.tpr -o new_run -c new_run...) and the energies decreased again, as 
they should. Next, I grompped a new tpr-file from the old gro, top and 
mdp files, got a new tpr file and simulated it resulting in decreasing 
energies again. But when I try to do the whole process from the 
beginning using the same pdb file as in the older run something goes 
differently and the mdrun produces rising energies. The energy 
minimisation with the new run takes over 3000 steps while the old run 
was minimised in 34 steps. I've checked that the gro and top files 
between the old (working) and the new (non-working) runs are identical. 
The only differences I find are in the tpr files of the energy 
minimisation and md simulation and they look like this (compared with diff):

Energy minimisation (hm2 is the new bad behaving simulation and hm5 is 
the old, good one)

1c1
< hm2_em.tpr:
---
 > hm5_em.tpr:
78c78
<       sc_power             = 0
---
 >       sc_power             = 2
89c89
<       dihre-fc             = 1000
---
 >       dihre-fc             = 1.4013e-42
142c142
<       scalefactor          = 1
---
 >       scalefactor          = 0
22983c22983
<             atom[     4]={type=  2, typeB=  2, ptype=    Atom, m= 
1.40067e+01, q= 4.80000e-02, mB= 1.40067e+01, qB= 4.80000e-02, resnr= 
  0} grpnrs=[ 0 0 0 0 0 0 0
  0 0 0 ]}
---
 >             atom[     4]={type=  1, typeB=  1, ptype=    Atom, m= 
1.40067e+01, q= 4.80000e-02, mB= 1.40067e+01, qB= 4.80000e-02, resnr= 
  0} grpnrs=[ 0 0 0 0 0 0 0
  0 0 0 ]}
22985c22985
<             atom[     6]={type=  3, typeB=  3, ptype=    Atom, m= 
1.59994e+01, q=-5.56000e-01, mB= 1.59994e+01, qB=-5.56000e-01, resnr= 
  0} grpnrs=[ 0 0 0 0 0 0 0
  0 0 0 ]}
---
 >             atom[     6]={type=  2, typeB=  2, ptype=    Atom, m= 
1.59994e+01, q=-5.56000e-01, mB= 1.59994e+01, qB=-5.56000e-01, resnr= 
  0} grpnrs=[ 0 0 0 0 0 0 0
  0 0 0 ]}

plus a lot of these atom definitions.


MD simulation

1c1
< hm2_md.tpr:
---
 > hm5_md.tpr:
78c78
<       sc_power             = 0
---
 >       sc_power             = 2
89c89
<       dihre-fc             = 1000
---
 >       dihre-fc             = 1.4013e-42
143c143
<       scalefactor          = 1
---
 >       scalefactor          = 0
160,11566c160,11566
<       x[    0]={ 2.47200e+00,  2.49800e+00,  2.51500e+00}
<       x[    1]={ 2.37100e+00,  2.47900e+00,  2.61100e+00}
<       x[    2]={ 2.39700e+00,  2.50800e+00,  2.74600e+00}
<       x[    3]={ 2.52300e+00,  2.55600e+00,  2.78600e+00}
<       x[    4]={ 2.62400e+00,  2.58300e+00,  2.69000e+00}
<       x[    5]={ 2.59800e+00,  2.54700e+00,  2.55500e+00}
<       x[    6]={ 2.69600e+00,  2.51500e+00,  2.48800e+00}
<       x[    7]={ 2.42800e+00,  2.46600e+00,  2.38500e+00}
<       x[    8]={ 2.29600e+00,  2.42000e+00,  2.38200e+00}
<       x[    9]={ 2.21300e+00,  2.41600e+00,  2.54400e+00}
<       x[   10]={ 2.49800e+00,  2.49100e+00,  2.26500e+00}
<       x[   11]={ 2.45600e+00,  2.39700e+00,  2.14900e+00}
<       x[   12]={ 2.31700e+00,  2.37400e+00,  2.14600e+00}
<       x[   13]={ 2.24100e+00,  2.36700e+00,  2.26400e+00}
<       x[   14]={ 2.73500e+00,  2.67900e+00,  2.71500e+00}
<       x[   15]={ 2.82900e+00,  2.63400e+00,  2.82800e+00}

plus a lot more coordination, velocity, atom type and interaction 
definitions.

I have no dihedral restraints applied and those sc_power, dihre-fc and 
scalepower definitions are different even in the old, working simulation 
and in the one which was grompped again and producing correct results 
with decreasing energies.

So the problem is that the old gro, top and mdp files give correct 
results even with grompped again, but a new procedure starting from the 
pdb file give bad results.

I was hoping that someone could direct me to the correct direction and 
find out what's wrong. All suggestions are highly appreciated and tested!


Kind regards,

Sampo Karkola


__________________________________

Sampo Karkola
Doctoral Student

Laboratory of Organic Chemistry
Department of Chemistry
Faculty of Science
PO Box 55 (A.I. Virtasen aukio 1)
FIN-00014 University of Helsinki
Finland

email: sampo.karkola at helsinki.fi

tel +358 9 19150369
fax +358 9 19150366