[gmx-users] Measuring the Influence of Harmonic Restraints

[Natalie Nguyen]

Dear All,


I am interested in measuring the change in free energy associated with applying and removing a harmonic potential to a ligand (glucose) as part of a protein-ligand system (GGBP + Glucose). I am using thermodynamic integration to measure this change; however I am having trouble setting up the production run to allow to output of dhdl files to use g_bar analysis on.



The following is a selection of the mdp used for the production run:


; Free energy parameters
free-energy              = yes
; charges as lambda goes to 0
sc-power                 = 1
sc-sigma                 = 0.3
sc-alpha                 = 0.5
; we still want the molecule to interact with itself at lambda=0 and lambda=1
couple-intramol          = no
couple-lambda0        = vdw-q
couple-lambda1        = vdw-q
init-lambda-state        = $LAMBDA$  (here I used a script to replace the lambda state value for each window)
; These are the lambda states at which we simulate
restraint-lambdas              = 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
foreign_lambda        = 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
nstdhdl            = 100
separate-dhdl-file    = yes
dhdl-derivatives    = yes
calc-lambda-neighbors    = -1

; Pull code
pull    = umbrella
pull_geometry    = distance
pull_dim    = N N Y
pull_start    = yes
pull_init1    = 0
pull_ngroups    = 1
pull_group0    = Protein
pull_group1    = BGC
pull_k1        = 0
pull_kB1        = 1000



Things I have tried to change:

-changing nstdhdl value

-editing pull code k values as well as using pull_start = yes and no

-removing dhdl-derivatives option

-setting a delta-lambda value


I am using the 4.6.5 Version of Gromacs.


If there is any other kind of information I need to provide, feel free to let me know.


Thanks to all ahead of time!


-Natalie Nguyen