# [gmx-users] Question about gmx bar

Justin Lemkul jalemkul at vt.edu
Wed Feb 28 00:07:00 CET 2018

```
On 2/27/18 4:13 PM, Thanh Le wrote:
> Hi everyone,
> I just finished running 2 sets of simulations (ligand in water and RNA+ligand in water) for my system to learn about its binding energy. Using the parameter for BAR method, I ran 20 simulations for ligand in water and 30 simulations for RNA+ligand in water with different lambda stages. What is interesting/confusing to me is the 0 energy from stages 0-10. Based on what I have been reading and the set up of my prod.mdp, these stages should give Coulomb Energy. If you can tell me why or how to fix it, I would greatly appreciate your help.

You've defined only a bonded transformation in those stages:

bonded-lambdas           = 0.0 0.01 0.025 0.05 0.075 0.1 0.2 0.35 0.5 0.75 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.0 1.0 1.0 1.0 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0
coul-lambdas             = 0.0 0.00 0.000 0.00 0.000 0.0 0.0 0.00 0.0 0.00 0.0 0.25 0.5 0.75 1.0 1.00 1.0 1.0 1.0 1.0 1.0 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0
vdw-lambdas              = 0.0 0.00 0.000 0.00 0.000 0.0 0.0 0.00 0.0 0.00 0.0 0.00 0.0 0.00 0.0 0.05 0.1 0.2 0.3 0.4 0.5 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1.0

Your electrostatic terms remain in state A (lambda = 0) and there is no difference in bonded energies, so the first windows are essentially doing nothing.

-Justin

> point      0 -      1,   DG  0.00 +/-  0.00
> point      1 -      2,   DG  0.00 +/-  0.00
> point      2 -      3,   DG  0.00 +/-  0.00
> point      3 -      4,   DG  0.00 +/-  0.00
> point      4 -      5,   DG  0.00 +/-  0.00
> point      5 -      6,   DG  0.00 +/-  0.00
> point      6 -      7,   DG  0.00 +/-  0.00
> point      7 -      8,   DG  0.00 +/-  0.00
> point      8 -      9,   DG  0.00 +/-  0.00
> point      9 -     10,   DG  0.00 +/-  0.00
> point     10 -     11,   DG 1198.03 +/- 11.41
> point     11 -     12,   DG 711.32 +/-  1.58
> point     12 -     13,   DG 258.19 +/-  3.80
> point     13 -     14,   DG -116.21 +/-  3.81
> point     14 -     15,   DG 24.24 +/-  0.17
> point     15 -     16,   DG 25.40 +/-  0.18
> point     16 -     17,   DG 51.57 +/-  0.42
> point     17 -     18,   DG 51.20 +/-  0.57
> point     18 -     19,   DG 48.42 +/-  0.97
> point     19 -     20,   DG 46.55 +/-  0.56
> point     20 -     21,   DG 44.53 +/-  0.23
> point     21 -     22,   DG 20.56 +/-  0.12
> point     22 -     23,   DG 18.23 +/-  0.11
> point     23 -     24,   DG 13.64 +/-  0.16
> point     24 -     25,   DG -3.58 +/-  0.84
> point     25 -     26,   DG -44.60 +/-  0.29
> point     26 -     27,   DG -33.45 +/-  0.07
> point     27 -     28,   DG -15.22 +/-  0.02
> point     28 -     29,   DG  0.81 +/-  0.12
>
> total      0 -     29,   DG 2299.63 +/-  5.17
>
> Here is the prod.mdp for complex:
>
> ;====================================================
> ; Production simulation
> ;====================================================
>
> ;----------------------------------------------------
> ; RUN CONTROL
> ;----------------------------------------------------
> integrator   = sd            ; stochastic leap-frog integrator
> nsteps       = 50000000        ; 2 * 500,000 fs = 1000 ps = 1 ns
> dt           = 0.002         ; 2 fs
> comm-mode    = Linear        ; remove center of mass translation
> nstcomm      = 100           ; frequency for center of mass motion removal
>
> ;----------------------------------------------------
> ; OUTPUT CONTROL
> ;----------------------------------------------------
> nstxout                = 0          ; don't save coordinates to .trr
> nstvout                = 0          ; don't save velocities to .trr
> nstfout                = 0          ; don't save forces to .trr
> nstxout-compressed     = 1000       ; xtc compressed trajectory output every 1000 steps (2 ps)
> compressed-x-precision = 1000       ; precision with which to write to the compressed trajectory file
> nstlog                 = 1000       ; update log file every 2 ps
> nstenergy              = 1000       ; save energies every 2 ps
> nstcalcenergy          = 100        ; calculate energies every 100 steps
>
> ;----------------------------------------------------
> ; BONDS
> ;----------------------------------------------------
> constraint_algorithm   = lincs      ; holonomic constraints
> constraints            = all-bonds  ; hydrogens only are constrained
> lincs_iter             = 1          ; accuracy of LINCS (1 is default)
> lincs_order            = 4          ; also related to accuracy (4 is default)
> lincs-warnangle        = 30         ; maximum angle that a bond can rotate before LINCS will complain (30 is default)
> continuation           = yes        ; formerly known as 'unconstrained-start' - useful for exact continuations and reruns
>
> ;----------------------------------------------------
> ; NEIGHBOR SEARCHING
> ;----------------------------------------------------
> cutoff-scheme   = Verlet
> ns-type         = grid   ; search neighboring grid cells
> nstlist         = 10     ; 20 fs (default is 10)
> rlist           = 1.0    ; short-range neighborlist cutoff (in nm)
> pbc             = xyz    ; 3D PBC
>
> ;----------------------------------------------------
> ; ELECTROSTATICS
> ;----------------------------------------------------
> coulombtype      = PME      ; Particle Mesh Ewald for long-range electrostatics
> rcoulomb         = 1.0      ; short-range electrostatic cutoff (in nm)
> ewald_geometry   = 3d       ; Ewald sum is performed in all three dimensions
> pme-order        = 6        ; interpolation order for PME (default is 4)
> fourierspacing   = 0.10     ; grid spacing for FFT
> ewald-rtol       = 1e-6     ; relative strength of the Ewald-shifted direct potential at rcoulomb
>
> ;----------------------------------------------------
> ; VDW
> ;----------------------------------------------------
> vdw-type                = PME
> rvdw                    = 1.0
> vdw-modifier            = Potential-Shift
> ewald-rtol-lj           = 1e-3
> lj-pme-comb-rule        = Geometric
> DispCorr                = EnerPres
>
> ;----------------------------------------------------
> ; TEMPERATURE & PRESSURE COUPL
> ;----------------------------------------------------
> tc_grps          = System
> tau_t            = 1.0
> ref_t            = 300
> pcoupl           = Parrinello-Rahman
> pcoupltype       = isotropic            ; uniform scaling of box vectors
> tau_p            = 2                    ; time constant (ps)
> ref_p            = 1.0                  ; reference pressure (bar)
> compressibility  = 4.5e-05              ; isothermal compressibility of water (bar^-1)
>
> ;----------------------------------------------------
> ; VELOCITY GENERATION
> ;----------------------------------------------------
> gen_vel      = no       ; Velocity generation is off
>
> ;----------------------------------------------------
> ; FREE ENERGY CALCULATIONS
> ;----------------------------------------------------
> free-energy              = yes
> couple-moltype           = Protein_chain_B
> couple-lambda0           = vdw-q
> couple-lambda1           = none
> couple-intramol          = yes
> separate-dhdl-file       = yes
> sc-alpha                 = 0.5
> sc-power                 = 1
> sc-sigma		             = 0.3
> init-lambda-state        = 0
> bonded-lambdas           = 0.0 0.01 0.025 0.05 0.075 0.1 0.2 0.35 0.5 0.75 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.0 1.0 1.0 1.0 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0
> coul-lambdas             = 0.0 0.00 0.000 0.00 0.000 0.0 0.0 0.00 0.0 0.00 0.0 0.25 0.5 0.75 1.0 1.00 1.0 1.0 1.0 1.0 1.0 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0
> vdw-lambdas              = 0.0 0.00 0.000 0.00 0.000 0.0 0.0 0.00 0.0 0.00 0.0 0.00 0.0 0.00 0.0 0.05 0.1 0.2 0.3 0.4 0.5 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1.0
> nstdhdl                  = 100
> calc-lambda-neighbors    = -1
>
>
> Here is the prod.mdp for ligand
>
> ;====================================================
> ; Production simulation
> ;====================================================
>
> ;----------------------------------------------------
> ; RUN CONTROL
> ;----------------------------------------------------
> integrator   = sd            ; stochastic leap-frog integrator
> nsteps       = 50000000        ; 2 * 500,000 fs = 1000 ps = 1 ns
> dt           = 0.002         ; 2 fs
> comm-mode    = Linear        ; remove center of mass translation
> nstcomm      = 100           ; frequency for center of mass motion removal
>
> ;----------------------------------------------------
> ; OUTPUT CONTROL
> ;----------------------------------------------------
> nstxout                = 0          ; don't save coordinates to .trr
> nstvout                = 0          ; don't save velocities to .trr
> nstfout                = 0          ; don't save forces to .trr
> nstxout-compressed     = 1000       ; xtc compressed trajectory output every 1000 steps (2 ps)
> compressed-x-precision = 1000       ; precision with which to write to the compressed trajectory file
> nstlog                 = 1000       ; update log file every 2 ps
> nstenergy              = 1000       ; save energies every 2 ps
> nstcalcenergy          = 100        ; calculate energies every 100 steps
>
> ;----------------------------------------------------
> ; BONDS
> ;----------------------------------------------------
> constraint_algorithm   = lincs      ; holonomic constraints
> constraints            = all-bonds  ; hydrogens only are constrained
> lincs_iter             = 1          ; accuracy of LINCS (1 is default)
> lincs_order            = 4          ; also related to accuracy (4 is default)
> lincs-warnangle        = 30         ; maximum angle that a bond can rotate before LINCS will complain (30 is default)
> continuation           = yes        ; formerly known as 'unconstrained-start' - useful for exact continuations and reruns
>
> ;----------------------------------------------------
> ; NEIGHBOR SEARCHING
> ;----------------------------------------------------
> cutoff-scheme   = Verlet
> ns-type         = grid   ; search neighboring grid cells
> nstlist         = 10     ; 20 fs (default is 10)
> rlist           = 1.0    ; short-range neighborlist cutoff (in nm)
> pbc             = xyz    ; 3D PBC
>
> ;----------------------------------------------------
> ; ELECTROSTATICS
> ;----------------------------------------------------
> coulombtype      = PME      ; Particle Mesh Ewald for long-range electrostatics
> rcoulomb         = 1.0      ; short-range electrostatic cutoff (in nm)
> ewald_geometry   = 3d       ; Ewald sum is performed in all three dimensions
> pme-order        = 6        ; interpolation order for PME (default is 4)
> fourierspacing   = 0.10     ; grid spacing for FFT
> ewald-rtol       = 1e-6     ; relative strength of the Ewald-shifted direct potential at rcoulomb
>
> ;----------------------------------------------------
> ; VDW
> ;----------------------------------------------------
> vdw-type                = PME
> rvdw                    = 1.0
> vdw-modifier            = Potential-Shift
> ewald-rtol-lj           = 1e-3
> lj-pme-comb-rule        = Geometric
> DispCorr                = EnerPres
>
> ;----------------------------------------------------
> ; TEMPERATURE & PRESSURE COUPL
> ;----------------------------------------------------
> tc_grps          = System
> tau_t            = 1.0
> ref_t            = 300
> pcoupl           = Parrinello-Rahman
> pcoupltype       = isotropic            ; uniform scaling of box vectors
> tau_p            = 2                    ; time constant (ps)
> ref_p            = 1.0                  ; reference pressure (bar)
> compressibility  = 4.5e-05              ; isothermal compressibility of water (bar^-1)
>
> ;----------------------------------------------------
> ; VELOCITY GENERATION
> ;----------------------------------------------------
> gen_vel      = no       ; Velocity generation is off (if gen_vel is 'yes', continuation should be 'no')
> gen_seed     = -1       ; Use random seed
> gen_temp     = 300
>
> ;----------------------------------------------------
> ; FREE ENERGY CALCULATIONS
> ;----------------------------------------------------
> free-energy              = yes
> couple-moltype           = Protein_chain_B
> couple-lambda0           = vdw-q
> couple-lambda1           = none
> couple-intramol          = yes
> separate-dhdl-file       = yes
> sc-alpha                 = 0.5
> sc-power                 = 1
> sc-sigma		             = 0.3
> init-lambda-state        = 0
> coul-lambdas             = 0.0 0.25 0.5 0.75 1.0 1.00 1.0 1.0 1.0 1.0 1.0 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0
> vdw-lambdas              = 0.0 0.00 0.0 0.00 0.0 0.05 0.1 0.2 0.3 0.4 0.5 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1.0
> nstdhdl                  = 100
> calc-lambda-neighbors    = -1
>

--
==================================================

Justin A. Lemkul, Ph.D.
Assistant Professor
Virginia Tech Department of Biochemistry

303 Engel Hall
340 West Campus Dr.
Blacksburg, VA 24061

jalemkul at vt.edu | (540) 231-3129
http://www.biochem.vt.edu/people/faculty/JustinLemkul.html

==================================================

```