[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
==================================================
More information about the gromacs.org_gmx-users
mailing list