[gmx-users] free energy perturbation

[Yuhua Song]

Hi, 

I am doing some FEP work with Gromacs. According to one reply in the
gromacs mailing list as attached below, I do the corresponding work, but
the results seem not right. Here is the procedure that I did my work,
could you please give a look to see what might be wrong?

1) First I run 2ns production for a NA binding to a protein. I get the
trajectory file: 1.trj

2) to study the free energy perturbation, we replace the binding NA with
LI. the modified part of the topology file is as following ( I use oplsa
force field, opls_406 is LI atom):

  4724   opls_272    295    GLU     O2   1496       -0.8    15.9994   ;
qtot 5
  4725   opls_407    296     NA     NA   1497          1    22.9898  
opls_406     1      6.9410

3) from the new topology file perturbation by LI, and the *.gro from
production run, I get the new *.tpr file

3) Then I do the mdrun -rerun, first rerun the *.tpr from original tpr
with NA binding, then do the rerun with the *.tpr with LI perturbation,
so I get two *.edr file 

4) then I do g_energy -f 1.edr -f2 2.edr -fee ....
The results show that difference of energy is only 0.53

5) Because I am not sure the results is right, then I repeat the above
procedure, but this time I do the free energy perturbation of K to
replace of NA instead of LI, but the result that I get is totally same
as the case of LI, So, I feel that something might be wrong here.

Could you please to take a look to see what might be wrong about the
procedure that I had done?

Thank you.

Yuhua 



 

________________________________

[gmx-users] g_energy -fee -f2 
Anton Feenstra  gmx-users at gromacs.org
Fri Sep 19 09:24:00 2003

      * Previous message: [gmx-users] g_energy -fee -f2
      * Next message: [gmx-users] g_energy -fee -f2
      * Messages sorted by: [ date ] [ thread ] [ subject ] [ author ]

________________________________________________________________________
Xavier Periole wrote:
> One changes the topology according to the perturbation and generate a>
new file.tpr and rerun. That what I would do. But I guess that if you
change> the number of atoms (I don't know about the type though) the
trajectory> will no longer correspond to the topology defined in the tpr
and the program> will complain.> At that point I think somebody with the
experience of FEP calculation> within gromacs would be of better use
than my comments.
You're right, FEP doesn't work if you change the number of atoms. There
are possibilities involving changing atoms to a 'non-interacting' type
(aka 'dummy' atoms), but that involves corrections to the entropy that
I'm not familiar with.

If you don't need to add or remove atoms, it's fairly simple, as I've
already briefly explained in my other mail. Make two topologies, that
e.g. differ in charges, atom types, bond lengths, angles, dihedrals.
Run your simulation with the one (non-perturbed) .tpr, then do a re-run
with the non-perturbed and a rerun with the perturbed .tpr. Feed those
energy files to g_energy. You do need to re-run the non-perturbed also
due to slight differences in 'original' and 're-run' energies that are
large enough to show up in the free energy estimate...

For additional accuracy (and an estimate of the error), you could also
run the reverse: simulate the perturbed system and then re-run both and
calculate the backward FEP.