[gmx-users] "B" state for free energy corrupts minimization?
Berk Hess
gmx3 at hotmail.com
Wed Feb 1 11:14:53 CET 2006
>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: Wed, 1 Feb 2006 00:55:44 -0800 (PST)
>
> >>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).
Indeed the zeros in the bonds (and probably also in the angles?) are the
problem.
This is caused by a bug that I have fixed in November,
but we did not release 3.3.1 yet:
http://www.gromacs.org/pipermail/gmx-revision/2005-November/000117.html
Apparently some bonded interaction parameters are not defined
for bonded interactions that involve your 36 atom type in the B-state.
You should add parameters for these bonded interactions,
either directly after the definition of the bonded interaction,
or in the force-field bondtypes and angletypes sections.
Berk.
More information about the gromacs.org_gmx-users
mailing list