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

[Mark Abraham]

Sampo Karkola wrote:
> Dear list,

Thanks for the good problem description!

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

So far everything looks good to me.

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

gmxcheck -s1 old.tpr -s2 new.tpr is the recommended way to check for 
this sort of thing.

> 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
> plus a lot of these atom definitions.

These seem strange to me - the "type" values in the above arrays should 
be indices into an atomtype array earlier in the .tpr file, but if the 
above values are right, and the two atoms 4 should correspond, then the 
atomtype array should also be different... but diff is not reporting 
this, and I do not understand why. Does gmxcheck differ? (Recall that 
this is not the fourth atom, since numbering starts at zero).

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

Evidently, since the starting configuration is embedded in the .tpr and 
the endpoint of the previous minimization was different.

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

Not much joy from me, I'm afraid!

Mark