[gmx-users] Doubts about g_lie and g_bar

[Varvdekar Bhagyesh Rajendra]

Dear all,

I have tried to find Delta G of Binding using g_lie and g_bar from the following systems (Protein+Ligand with ligand = chain B):
PDB Ids: 1CDL, 1CKK, 1NIW, 2L7L, 2o60.

I have encountered following 2 doubts while doing so:

1. This doubt may have been answered in the past threads but I couldn't find a precise answer for it. So, I would appreciate if this is answered here again in more details. I have used PME in the production runs of my systems for finding Delta G_Binding using g_lie. But many threads on Gromacs forum warn against it and advise to rerun the whole trajectory using Reaction field Zero, but the cutoff for this rerun is not clearly given. Some threads say a cutoff with large cutoff using RF-0 is necessary. Can anyone please answer as to what cutoff(Relative to the system) should be given while re-running the trajectory and what all changes have to be done in the original mdp file(which uses PME). Also, a brief explanation for the whole process for the need of the re-run would be beneficial. Following are the changes I made after studying the past threads:
------------------------------------------------------------------------------------------------------------------------
;; With PME  ;;;                                     | ;; Re-Run without PME ;;
-------------------------------------------------------------------------------------------------------------------------
 ; Neighborsearching                                 |  ; Neighborsearching
  ns_type     = grid      ; search neighboring grid  |  ns_type     = grid      ; search neighboring grid cells
  nstlist     = 5         ; 10 fs                    |  nstlist     = 5         ; 10 fs
  rlist       = 0.9       ; short-range neighborlist |  rlist       = 1.2       ; short-range neighborlist cutoff (nm)                                                                                            
  rcoulomb    = 0.9       ; short-range electrostatic|  rcoulomb    = 0.9       ; short-range electrostatic cutoff (in nm)
  rvdw        = 1.4       ; short-range van der Waals|  rvdw        = 1.4       ; short-range van der Waals cutoff (in nm)
                                                     |
; Electrostatics                                     | ; Electrostatics
coulombtype     = PME                                | coulombtype  =  Reaction-Field-zero 
                                                     | epsilon_rf      =  0
---------------------------------------------------------------------------------------------------------------------------

2. Free Energy tutorials of Gromacs suggest using oriental restraints and/or distance restraints on protein-ligand systems for calculating Delta G_Binding using g_bar. But the exact details are not given on the website. One can see that in my systems the protein is fully surrounding the ligand except two entrances. So, I ask if distance and/or oriental constraint is necessary to keep from floating it away. If yes, can anyone suggest some simple ways for doing so in Gromacs. I found some research papers who employ gruesome techiques of finding the correct atoms to constraint in the ligand (for eg.Theoretical Study of the Binding Profile of the Allosteric Modulator NS-1738 with a Chimera Structure of the α 7 Nicotinic Acetylcholine Receptor).

I would be grateful if you could answer the above doubts.

Thanking you in anticipation,

Best Regards,

Bhagyesh Varvdekar, 
Undergrad student,
IIIT Hyderabad, India.