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. 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 -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
