[gmx-users] "B" state for free energy corrupts minimization?

[mernst at tricity.wsu.edu]

>>From: mernst at tricity.wsu.edu
>>Reply-To: Discussion list for GROMACS users <gmx-users at gromacs.org>
>>To: gmx-users at gromacs.org
>>Subject: [gmx-users] "B" state for free energy corrupts minimization?
>>Date: Tue, 31 Jan 2006 00:16:11 -0800 (PST)
>>
>>I am attempting to perform a free energy simulation that involves the
>>transformation of the lesioned DNA base 8-oxoguanine to guanine. Before I
>>set up the B state in my .itp file, I was able to simulate 8-oxoguanine
>>with MD and its behavior looked to be in reasonable agreement with
>>simulations I have performed in another package.
>>
>>Now that I have defined B-state parameters, MD and even energy
>>minimization give strange and undesirable results, even before I begin
>> the
>>free energy calculations or set lambda to a non-zero value.
>>
>>I've narrowed it down to a simple test case where I supply a B state for
>>just one atom.
>>
>>If I use this line in my .itp file, which changes nothing between A and
>> B,
>>the energy minimization gives a reasonable-looking structure:
>>
>>46 amber99_35 2 8oG  N7   46   -0.5129    14.01   amber99_35 -0.5129
>> 14.01
>>
>>This is what the output looks like:
>>
>>http://www.sciencemadness.org/scipics/goodtrimer.png
>>
>>If I use this line, using a different atom type in the B state,
>>
>>46 amber99_35 2 8oG  N7   46   -0.5129    14.01   amber99_36 -0.5129
>> 14.01
>>
>>the minimization gives me a hydrogen that sticks out at a strange angle:
>>
>>http://www.sciencemadness.org/scipics/badtrimer.png
>>
>>The amber99_35 atom is AMBER type NA, corresponding to sp2 N in 5
>>memb.ring w/H atom (HIS). The amber99_36 atom is AMBER type NB,
>>corresponding to sp2 N in 5 memb.ring w/LP (HIS,ADE,GUA). In the A state
>>there is a hydrogen that becomes a dummy in the B state, so I would like
>>to go from NA to NB. The AMBER force field port to Gromacs is from the
>>Pande group.
>>
>>I cannot figure out why this effect shows up when lambda=0, and how I
>>might set up the free energy simulation so that it is not horribly
>>distorted and prone to crashing in the equilibrium period.
>
> This sounds very strange.
> B-state parameters should have no influence on the forces at lambda=0.
> I have done many free-energy simulations with the Amber force-field
> without problems.
>
> Three questions:
>
> Have you tried the free-energy topology with free-energy in the mdp
> file turned on and off, and is there a difference?
>
> Could you compare the normal and free-energy potential energies
> for the first step of the minimization?
> This should tell you in which term the difference is located.
>
> You can compare tpr files with gmxcheck -s1 ... -s2 ...
> You should do this for your normal and free-energy tpr.
> You might get a lot of output if the 36 atom type is not at all
> present in the A-state, but we are probably looking for
> differences in the bonded interactions.
>
> Berk.
Berk,

Thank you for your hints. I have tried setting free_energy=yes and free_energy=no in the
mdp file and saw no difference.

When the nitrogen N7 changes to type amber99_36, this is the energy at the first
minimization step:

          Step           Time         Lambda
              0        0.00000        0.00000

   Energies (kJ/mol)
           Bond          Angle    Proper Dih. Ryckaert-Bell.          LJ-14
    6.95106e+02    1.95648e+02    3.17244e-01    5.42812e+02    2.67505e+02
     Coulomb-14        LJ (SR)        LJ (LR)   Coulomb (SR)   Coul. recip.
   -1.91244e+03    2.46390e+04   -4.41360e+02   -1.59749e+05   -1.50672e+04
      Potential    Kinetic En.   Total Energy    Temperature Pressure (bar)
   -1.50830e+05    0.00000e+00    0.00000e+00    0.00000e+00    0.00000e+00

When no atoms change type, this is the energy at the first minimization step:

           Step           Time         Lambda
              0        0.00000        0.00000

   Energies (kJ/mol)
           Bond          Angle    Proper Dih. Ryckaert-Bell.          LJ-14
    6.96918e+02    1.99483e+02    3.17244e-01    5.42812e+02    2.67505e+02
     Coulomb-14        LJ (SR)        LJ (LR)   Coulomb (SR)   Coul. recip.
   -1.91244e+03    2.46390e+04   -4.41360e+02   -1.59749e+05   -1.50672e+04
      Potential    Kinetic En.   Total Energy    Temperature Pressure (bar)
   -1.50824e+05    0.00000e+00    0.00000e+00    0.00000e+00    0.00000e+00

It seems most of the difference is in potential, followed by angle and finally bond.

When I run gmxcheck I get a large output which can be seen here:

http://www.sciencemadness.org/scipics/trimer-topology-differences.txt

Here's part of the output:

idef->ntypes (458 - 462)
idef->iparam1: b0A= 1.35850e-01, cbA= 4.42667e+05, b0B= 1.35850e-01, cbB= 4.42667e+05
idef->iparam2: b0A= 1.35850e-01, cbA= 4.42667e+05, b0B= 1.30400e-01, cbB= 4.42667e+05
idef->iparam1: b0A= 1.41900e-01, cbA= 3.74050e+05, b0B= 1.41900e-01, cbB= 3.74050e+05
idef->iparam2: b0A= 0.00000e+00, cbA= 0.00000e+00, b0B= 0.00000e+00, cbB= 0.00000e+00
idef->iparam1: b0A= 1.38800e-01, cbA= 3.49782e+05, b0B= 1.38800e-01, cbB= 3.49782e+05
idef->iparam2: b0A= 1.41900e-01, cbA= 3.74050e+05, b0B= 1.41900e-01, cbB= 3.74050e+05
idef->iparam1: b0A= 1.38100e-01, cbA= 3.57314e+05, b0B= 1.38100e-01, cbB= 3.57314e+05
idef->iparam2: b0A= 1.38800e-01, cbA= 3.49782e+05, b0B= 1.38800e-01, cbB= 3.49782e+05
idef->iparam1: b0A= 1.36500e-01, cbA= 3.74886e+05, b0B= 1.36500e-01, cbB= 3.74886e+05
idef->iparam2: b0A= 1.38100e-01, cbA= 3.57314e+05, b0B= 1.38100e-01, cbB= 3.57314e+05
idef->iparam1: b0A= 1.38300e-01, cbA= 3.54803e+05, b0B= 1.38300e-01, cbB= 3.54803e+05
idef->iparam2: b0A= 1.36500e-01, cbA= 3.74886e+05, b0B= 1.36500e-01, cbB= 3.74886e+05
idef->iparam1: b0A= 1.35000e-01, cbA= 4.59403e+05, b0B= 1.35000e-01, cbB= 4.59403e+05
idef->iparam2: b0A= 1.38300e-01, cbA= 3.54803e+05, b0B= 1.38300e-01, cbB= 3.54803e+05
idef->iparam1: b0A= 1.51000e-01, cbA= 2.65266e+05, b0B= 1.51000e-01, cbB= 2.65266e+05
idef->iparam2: b0A= 1.35000e-01, cbA= 4.59403e+05, b0B= 1.35000e-01, cbB= 4.59403e+05
idef->iparam1: b0A= 1.44400e-01, cbA= 3.43088e+05, b0B= 1.44400e-01, cbB= 3.43088e+05
idef->iparam2: b0A= 1.51000e-01, cbA= 2.65266e+05, b0B= 1.51000e-01, cbB= 2.65266e+05
idef->iparam1: b0A= 1.43300e-01, cbA= 3.57314e+05, b0B= 1.43300e-01, cbB= 3.57314e+05
idef->iparam2: b0A= 1.44400e-01, cbA= 3.43088e+05, b0B= 1.44400e-01, cbB= 3.43088e+05
idef->iparam1: b0A= 1.35800e-01, cbA= 3.82418e+05, b0B= 1.35800e-01, cbB= 3.82418e+05
idef->iparam2: b0A= 1.43300e-01, cbA= 3.57314e+05, b0B= 1.43300e-01, cbB= 3.57314e+05

All the iparam2 values correspond to the system where the nitrogen atom type
transformation takes place. I am not very familiar with gmxcheck output, but the line of
parameters with all zeros (4th down) looks odd to me. At first glance it looks like many
other parameters changed too, it seems that gmxcheck was just thrown off in pairing up
parameters due to the atom type change. A similar line of all-zero values shows up later
beneath
idef->functype[320] (7 - 0). These lines could well indicate the source of my problem,
since the "A" state values were all zero as well. I don't know how these lines are being
generated and how I might instead force generation of correct values (whatever those may
be).

The only warning I have received was when using genbox:

WARNING: masses will be determined based on residue and atom names,
         this can deviate from the real mass of the atom type

Matt Ernst
Washington State University