[Fwd: Re: [gmx-users] Energies rising]

[Sampo Karkola]

Hi,

thanks for the answer.

David van der Spoel wrote:
> Sampo Karkola wrote:
>> Dear list,
>>
>> earlier I posted a message on the list about rising energies in md 
>> simulation and I still have not found a solution (see a thread with 
>> the header "incomparable results with files    produced with 
>> gromacs3.2.1 and 3.3.1"). I observed that the energies rise in any 
>> simulation I'm trying, i.e. any ligand or protein whether simulated in 
>> complex or alone. I got one suggestion that the temperature (entropy) 
>> could be the reason for weird behaviour of the potential energy. I'm 
>> simulating at 310K, and I studied the effect of the temperature and 
>> found the following:
>>
>> -at 50 K the potential energy decreases as it should  and as it has 
>> done earlier at 310K
>>
>> -at 200K the potential energy rises immediately at the beginning, but 
>> after 0.2 ps, the energy decreases below the starting point
>>
>> -at 250 the potential energy rises and stays high
>>
>> I would like to simulate the system at body temperature. Any 
>> suggestions how to do that?
>>
> Your results can not be judged from this information. Rising energies 
> from the minimum are fine, as long as they equilibrate to some value 
> (see tutor example).
> 
> But you also say that you get different results with 3.2.1 and 3.3.1, is 
> that correct?

With 3.2.1 I did an md test for a ligand I was going to use in a
simulation for a theoretical model of a protein with a ligand (self made 
topology). I observed decreasing energy in md simulation for both the 
ligand alone and for the complex. While doing a similar simulation for 
the complex of the same model and the same ligand in 3.3.1, I observed 
decreasing energy.
Then, I started to work on a crystal structure of crystal structure 
protein. I tested another self made ligand and the protein alone and 
observed rising energies in both cases. I guess this should not be the 
case, even if it is fine to have rising and equilibrating potential 
energy. With a theoretical protein model and with Gromacs 3.3.1, the 
energy decreases nicely. That's why I expected a similar behaviour for a 
crystal structure as well.

> 
> In that case you can check things by doing two tests:
> A: making a tpr file in 3.2.1 and then running it with 3.3.1 and 3.2.1
> B: making a tpr in both versions and compare it using gmxcheck -s2 -s1

I do not have 3.2.1 installed anymore and do not feel like re-installing
it if it is not absolutely necessary, but I do have the files from 3.2.1
simulation left.

A) the tpr file from 3.2.1 simulated earlier in 3.2.1 and now in 3.3.1 
give identical, decreasing energies whether run in 3.2.1 or 3.3.1. So it 
is probably not a bug.

B) the comparison between the old, working tpr and the new non-working 
one gives the following differences:
comparing inputrec
inputrec->nkx (32 - 22)
inputrec->nky (32 - 24)
inputrec->nkz (32 - 24)
inputrec->sc_power (0 - 2)
inputrec->dihre_fc (1.000000e+03 - 1.401298e-42)
inputrec->grpopts.nrdf[1] (1.025102e+04 - 3.783054e+03)
.
.
comparing ilist SETTLE
SETTLE->nr[0] (3418 - 1262)
SETTLE->multinr[0] (3418 - 1262)
Comparing radically different topologies - SETTLE->iatoms is wrong
.
.
comparing atoms
atoms->nr (5151 - 1917)
atom.type[1917] (7 - 33708)
atom.ptype[1917] (0 - 33712)
atom.resnr[1917] (632 - 135103412)
atom.m[1917] (1.599940e+01 - 0.000000e+00)
atom.q[1917] (-8.200000e-01 - 1.075216e-41)
atom.typeB[1917] (7 - 2061)
atom.mB[1917] (1.599940e+01 - 4.258539e-34)
atom.qB[1917] (-8.200000e-01 - 4.258541e-34)
atom.grpnr(0)[1917] (1 - 184)
atom.grpnr(1)[1917] (1 - 131)
atom.grpnr(2)[1917] (0 - 13)
atom.grpnr(3)[1917] (0 - 8)
atom.grpnr(4)[1917] (0 - 188)
atom.grpnr(5)[1917] (0 - 131)
atom.grpnr(6)[1917] (0 - 13)
atom.grpnr(7)[1917] (0 - 8)
atom.grpnr(8)[1917] (0 - 192)
atom.grpnr(9)[1917] (0 - 131)

plus a lot more these differences with "atom."-prefix. The number of 
atoms is different, because I increased the box size (and the number of 
water molecules) in the test run with 3.3.1. I've already tried 
different box sizes and used same mdp-files for both simulations but it 
does not seem to matter.

I'm just confused because the procedure I followed [pdb2gmx, editconf, 
genbox, grompp (for em), mdrun (for em), grompp (for md), mdrun (for 
md)] was similar with both versions and with 3.3.1 I get weird behaviour.

Any ideas? I'd be mostly grateful if you could point me to the right 
direction.

Kind regards,

Sampo


> 
> If A gives you different results then it is a problem, submit a bugzilla.
> If B gives you different tpr then check your input.
>