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?
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
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 HELPAll the best.:)
_______________________________________________
gmx-users mailing list [email protected]
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 [EMAIL PROTECTED].
Can't post? Read http://www.gromacs.org/mailing_lists/users.php
_______________________________________________ gmx-users mailing list [email protected] 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 [EMAIL PROTECTED] Can't post? Read http://www.gromacs.org/mailing_lists/users.php
