[gmx-users] FEP

[P]

Hi all.
I'm new to gromacs, and I would be glad to get some advice from more experienced usres.
I would like to estimate the difference in free energy between two systems.
First system (A) consists of DPPC membrane, TM-protein and some Zn ions
interacting with the protein. Second (B) has the same elements but Zn ions
are far away from the protein (in the solvent). Is it possible to estimate
free energy difference between those two systems in gromacs?
I thought about using dummy atoms and slow-growth methods for free energy calculations.
During the simulation, Zn ions from system (A) would be perturbated into dummy atoms
(charge=0, mass = 0) and at the same time equal number of dummy atoms in the solvent
would be perturbated to Zn atoms? Is it possible to achieve?
For the beginning I tried to estimate the difference in free energy between two systems.
1) first system consisted of two Cl- and one Zn2+ ions + solvent
2) second only solvent molecules
During the 100 ps MD I perturbated chare and mass of Cl and Zn
ions to zero. The energy values that I've obtained are wrong 960200,46.
I'm new to gromacs so I would appreciate your help a lot.
Here are my input files (I'm using ffgmx )

--------------md.mdp-------------
pp = /lib/cpp
define = -DFLEX_SPC
constraints = none
integrator = md
dt = 0.001
nsteps = 100000
nstxout = 100000
nstvout = 100000
nstfout = 100000
nstlog = 100
nstenergy = 100
nstxtcout = 100
nstlist = 10
ns_type = grid
rlist = 1.0
free_energy = yes
init_lambda = 0
delta_lambda = 0.00001
sc-alpha = 1.5
sc-sigma = 0.3
....
...
....
---------------------------------------
--------ZN.itp-----------------------
[ moleculetype ]
; molname nrexcl
ZnB 1

[ atoms ]
; id at type res nr residu name at name cg nr charge typeB chargeB massB
1 Zn 1 Zn Zn 1 2 65.37000 Zn 0 0

---------------------------------------
-----------Cl.itp-----------------
[ moleculetype ]
; molname nrexcl
ClB 1

[ atoms ]
; id at type res nr residu name at name cg nr charge typeB chargeB massB
1 Cl 1 Cl Cl 1 -1 35.45300 Cl 0 0

----------------------------



THANK YOU FOR HELP
All the best.
:)

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20060323/14eba0fd/attachment.html>

[David Mobley]

Hi,

Hi all.

Post by P
I'm new to gromacs, and I would be glad to get some advice from more experienced usres.
I would like to estimate the difference in free energy between two systems.
First system (A) consists of DPPC membrane, TM-protein and some Zn ions
interacting with the protein. Second (B) has the same elements but Zn ions
are far away from the protein (in the solvent). Is it possible to estimate
free energy difference between those two systems in gromacs?
I thought about using dummy atoms and slow-growth methods for free energy calculations.

During the simulation, Zn ions from system (A) would be perturbated into

Post by P
dummy atoms
(charge=0, mass = 0) and at the same time equal number of dummy atoms in the solvent
would be perturbated to Zn atoms? Is it possible to achieve?

I am rather nervous about doing FEP/TI for disappearing charged molecules,
as it is not at all clear to me that it is possible do this correctly with
current methods. Perhaps someone else may be able to comment more, but at
least with long range electrostatics (PME), systems are required to be
neutral. It IS possible to apply PME to a "charged" system, but I think this
is typically done by spreading a uniform neutralizing charge throughout the
box. It isn't clear that this is really correct for FEP/TI, I don't think,
since in that case you will have a charged molecule overlapping with some of
the uniform neutralizing charge, which could give strange energy artifacts.

Perhaps things are better if you don't use PME electrostatics, but there may
be boundary effects.

I need to personally dig in to the literature on this a bit more, but I
would suggest doing so yourself unless you're very sure that doing FEP/TI
will give you something meaningful on this system.

Also, if you *do* go that route, I would recommend against slow growth, as
it tends to be hysteretic. It's probably better to run a bunch of separate
simulations at different lambda values.

And even further, it is probably *not* a good idea to turn off the charges
and LJ interactions at the same time. It works better to turn off the
electrostatics first (without using soft core, as this generally is fairly
smooth with the charge), and then turn off the LJ interactions using soft
core (recommend sc-alpha=0.5, as sc-alpha=1.5 is usually too large for LJ).

However, if I were doing it, I would probably try to do it by computing the
PMF for pulling the zinc ions (or one zinc ion, which is easier) away from
the protein. You can compute the free energy difference from the PMF. The
unfortunate thing about this is that it involves sampling a lot of stuff you
don't care about, but at least it is clear the electrostatics is correct.
Alternatively, dig into the literature and try and figure out if there is a
correct way to do the electrostatics in this case. (And please let me know
what you find out, if you do!)

Oh, and one more thing: If DO compute the free energy by annhilating the
zinc ions, you probably should use restraints to keep them in a particular
region as you turn off the interactions, otherwise they will have to sample
the entire simulation box in order for you to get converged results.

Good luck.
David

Post by P
For the beginning I tried to estimate the difference in free energy between two systems.
1) first system consisted of two Cl- and one Zn2+ ions + solvent
2) second only solvent molecules
During the 100 ps MD I perturbated chare and mass of Cl and Zn
ions to zero. The energy values that I've obtained are wrong 960200,46.
I'm new to gromacs so I would appreciate your help a lot.
Here are my input files (I'm using ffgmx )
--------------md.mdp-------------
pp = /lib/cpp
define = -DFLEX_SPC
constraints = none
integrator = md
dt = 0.001
nsteps = 100000
nstxout = 100000
nstvout = 100000
nstfout = 100000
nstlog = 100
nstenergy = 100
nstxtcout = 100
nstlist = 10
ns_type = grid
rlist = 1.0
free_energy = yes
init_lambda = 0
delta_lambda = 0.00001
sc-alpha = 1.5
sc-sigma = 0.3
....
...
....
---------------------------------------
--------ZN.itp-----------------------
[ moleculetype ]
; molname nrexcl
ZnB 1
[ atoms ]
; id at type res nr residu name at name cg nr
charge typeB chargeB massB
1 Zn 1 Zn Zn 1 2 65.37000
Zn 0 0
---------------------------------------
-----------Cl.itp-----------------
[ moleculetype ]
; molname nrexcl
ClB 1
[ atoms ]
; id at type res nr residu name at name cg nr charge typeB chargeB massB
1 Cl 1 Cl Cl 1 -1 35.45300
Cl 0 0
----------------------------
THANK YOU FOR HELP
All the best.
:)
_______________________________________________
gmx-users mailing list gmx-users at gromacs.org
http://www.gromacs.org/mailman/listinfo/gmx-users
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-request at gromacs.org.
Can't post? Read http://www.gromacs.org/mailing_lists/users.php

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20060323/e00abe7a/attachment.html>

[Mark Abraham]

Post by David Mobley
I am rather nervous about doing FEP/TI for disappearing charged
molecules, as it is not at all clear to me that it is possible do this
correctly with current methods. Perhaps someone else may be able to
comment more, but at least with long range electrostatics (PME), systems
are required to be neutral. It IS possible to apply PME to a "charged"
system, but I think this is typically done by spreading a uniform
neutralizing charge throughout the box. It isn't clear that this is
really correct for FEP/TI, I don't think, since in that case you will
have a charged molecule overlapping with some of the uniform
neutralizing charge, which could give strange energy artifacts.

In general, PME can be applied to a system with any charge. The PME
section of the gromacs manual is silent on this topic, and I've never
tried to do it, however. You are correct about the uniform neutralizing
charge throughout the box, but this is present for *every* partial
charge in the system, each of which contributes individually in the
Ewald summations. The utility of that as a physical model is a different
question, however.

I can't comment on the effectiveness of PME on such a system with FEP/TI.

Mark

[David Mobley]

Mark,

Post by David Mobley
I am rather nervous about doing FEP/TI for disappearing charged

Post by David Mobley
molecules, as it is not at all clear to me that it is possible do this
correctly with current methods. Perhaps someone else may be able to
comment more, but at least with long range electrostatics (PME), systems
are required to be neutral. It IS possible to apply PME to a "charged"
system, but I think this is typically done by spreading a uniform
neutralizing charge throughout the box. It isn't clear that this is
really correct for FEP/TI, I don't think, since in that case you will
have a charged molecule overlapping with some of the uniform
neutralizing charge, which could give strange energy artifacts.

In general, PME can be applied to a system with any charge. The PME
section of the gromacs manual is silent on this topic, and I've never
tried to do it, however. You are correct about the uniform neutralizing
charge throughout the box, but this is present for *every* partial
charge in the system, each of which contributes individually in the
Ewald summations. The utility of that as a physical model is a different
question, however.

Thanks for the comments. Just to clarify what I meant about energy
artifacts, if you are compute the free energy of turning off a charge which
overlaps with a uniform neutralizing charge, suddenly the density of that
uniform neutralizing charge becomes quite important (at least, it seems like
it would) since the free energy of turning off the charge depends on the
surrounding charges... That's the part that seems concerning. (More so for
when the neutralizing charge overlap with a given atom is constant, as in
PME when FEP/TI is NOT done -- in which case it seems clear it doesn't
contribute to the dynamics or free energies in any way).

Thanks,
David


I can't comment on the effectiveness of PME on such a system with FEP/TI.

Post by David Mobley
Mark
_______________________________________________
gmx-users mailing list gmx-users at gromacs.org
http://www.gromacs.org/mailman/listinfo/gmx-users
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-request at gromacs.org.
Can't post? Read http://www.gromacs.org/mailing_lists/users.php

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20060323/ed15c986/attachment.html>

[paloureiro]

Hi David,

I am not an expert in FEP, but I have used that method to try to answer the
following question: if there was a free energy difference between a
system with
a Na+ ion inserted by the best-potential method or randomly.
I performed a series of MD runs (several lambdas) with annihilation of the ion
(with reaction field electrostatics) and, then, another series with
slow-growth
of the ion, starting with the last frame of the previous series of simulations
(Please, see the attached file - fep.jpg)
The problem was motivated by the fact that a Na+ ion was "inserted" in
a rather
dehydrated environment - and I questioned the physical/biological relevance of
that.
Maybe the errors associated with the rection field method are cancelled
as long
as one is using a delta(deltaG/deltaLambda).
But I agree with you that a PMF will be a better choice for the particular
system under discussion.

Regards.

Pedro.

Post by David Mobley
I am rather nervous about doing FEP/TI for disappearing charged molecules,
as it is not at all clear to me that it is possible do this correctly with
current methods. Perhaps someone else may be able to comment more, but at
least with long range electrostatics (PME), systems are required to be
neutral. It IS possible to apply PME to a "charged" system, but I think this
is typically done by spreading a uniform neutralizing charge throughout the
box. It isn't clear that this is really correct for FEP/TI, I don't think,
since in that case you will have a charged molecule overlapping with some of
the uniform neutralizing charge, which could give strange energy artifacts.
Perhaps things are better if you don't use PME electrostatics, but there may
be boundary effects.
I need to personally dig in to the literature on this a bit more, but I
would suggest doing so yourself unless you're very sure that doing FEP/TI
will give you something meaningful on this system.
Also, if you *do* go that route, I would recommend against slow growth, as
it tends to be hysteretic. It's probably better to run a bunch of separate
simulations at different lambda values.
And even further, it is probably *not* a good idea to turn off the charges
and LJ interactions at the same time. It works better to turn off the
electrostatics first (without using soft core, as this generally is fairly
smooth with the charge), and then turn off the LJ interactions using soft
core (recommend sc-alpha=0.5, as sc-alpha=1.5 is usually too large for LJ).
However, if I were doing it, I would probably try to do it by computing the
PMF for pulling the zinc ions (or one zinc ion, which is easier) away from
the protein. You can compute the free energy difference from the PMF. The
unfortunate thing about this is that it involves sampling a lot of stuff you
don't care about, but at least it is clear the electrostatics is correct.
Alternatively, dig into the literature and try and figure out if there is a
correct way to do the electrostatics in this case. (And please let me know
what you find out, if you do!)
Oh, and one more thing: If DO compute the free energy by annhilating the
zinc ions, you probably should use restraints to keep them in a particular
region as you turn off the interactions, otherwise they will have to sample
the entire simulation box in order for you to get converged results.
Good luck.
David

----------------------------------------------------------------
This message was sent using IMP, the Internet Messaging Program.

-------------- next part --------------
A non-text attachment was scrubbed...
Name: fep.JPG
Type: image/pjpeg
Size: 30425 bytes
Desc: not available
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20060323/7ce546e7/attachment.bin>

[Chris Oostenbrink]

Hello,

In principle, I think one can do this with FEP and TI, but there will be
several things to take into account. What you are basically doing is
calculating the free energy of 'solvating' the ions in water and in a
protein environment. In water, the solvation of a singly charged sodium
ion will give you roughly 480 kJ/mol. In the protein, I would expect you
get something of the same order, and then you subtract these two big
numbers and it will be difficult to get an accurate estimate of the
difference.

There will indeed be many electrostatic effects. I have seen a nice
manuscript by Kastenholz and Huenenberger which is in press with J.
Chem. Phys. (you'll find it on www.igc.ethz.ch -> publications). They
have established exactly which corrections one would have to apply to
calculate the solvation free energy of an ion using different methods
(Straight cutoff, reaction field or lattice sum methods). Some of these
correction terms may cancel if you look at a free energy difference,
others may not...

As a general comment: don't use slow growth, as that is a
non-equilibrium method and can by definition not give correct results
(if you do a single run). Depending on how you would do the
PMF-calculation, you will run into the same problem. If you use soft
atoms, then you do not really have to split up the Van der Waals and the
electrostatics.

Good luck,

Chris

[Michal Kolinski]

Hi all.
1)

Thank you all for comments and help on FEP issue.
I still have problems with setting up the system.
How can I obtain topology for DPPC lipids to use with GROMOS96
force field for protein membrane simulations. Should I
stay with old ffgmx force field, is it the best choice.
Is it possible to adopt old lipid parameters (lipid.itp
for ffgmx) for new GROMS96 force field.

2)

During FEP calculations I need to
slowly turn off charges and LJ interactions for Zn2+ ion.
In my dual topology for Zn2+ I have:


[ moleculetype ]
; molname nrexcl
ZnB 1

[ atoms ]
; id at type res nr residu name at name cg nr charge typeB chargeB massB
1 Zn 1 Zn Zn 1 2 65.37000 Zn (?) 0 0

and that will only perturbate charges on the Zn atom.
In order to turn off Zn ion LJ interactions I need to
change it into dummy atom. I'm using ffgmx forcefield and
which atom type should I use for above Zn(?). I checked
ffgmx.atp and there is no atom type DUM (like in GROMOS96).
Does it mean I have to add new atom type DUM in all ffgmx.atm,
ffgmxnbon.itp.. files. Please give my some suggestions.
Thank you for your time and help.
All best.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20060324/fd5182c1/attachment.html>