[gmx-users] Calculation of free energy of distance restraints alone

[ZACHARY GVILDYS]

Hello Everybody,

I am trying to determine the free energy of restraining and
"un-restraining" a protein-RNA complex while leaving on vdw and coulomb
interactions. My current free energy parameters are the following for one
of the mdp files in the lambda_state;
; Free energy control stuff
free_energy              = yes
init_lambda_state        = 18
delta_lambda             = 0
calc_lambda_neighbors    = 1        ; only immediate neighboring windows
; Vectors of lambda specified here
; Each combination is an index that is retrieved from init_lambda_state for
each simulation
; init_lambda_state        0    1    2    3    4    5    6    7    8    9
 10   11   12   13   14   15   16   17   18   19   20
vdw_lambdas              = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
coul_lambdas             = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
bonded_lambdas           = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
restraint_lambdas        = 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40
0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00
mass_lambdas             = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
temperature_lambdas      = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
couple-moltype           = RNA  ; name of moleculetype to decouple
couple-lambda0           = vdw-q     ; only van der Waals interactions
couple-lambda1           = vdw-q     ; turn off everything, in this case
only vdW
couple-intramol          = yes
nstdhdl                  = 10


I have couple-inramol turned on because the RNA has distance restraints
along its backbone. What I am curious about is what the couple-lambda0 and
couple-lambda1 values should be. There are no restraint options, as far
that I know of, so I have both vdw and coulomb interactions on in both
states. However, when I ran this simulation and performed the gmx bar
function, I get a change in free energy of 0 and a warning that the initial
and final states are the same. Is there something that I am missing here?
Is it possible to determine the free energy from restraints and nothing
else using the built in GROMACS software? I am comfortable using my own
methods to calculate the free energy of restraining the system, however it
is less accurate than using the BAR method gromacs has implemented.

Thanks in advance for all the help!

Best,
Zach Gvildys