GROMACS version: 4.6.7
GROMACS modification: No
Here post your question
Dear all,
Sorry to bother you. I have been working on studying the interaction energy of ligand to its dna target. The idea is using the potential of mean force to evaluate the dissociation of the ligand along the reaction coordinate (the distance between dna and its binding ligand). To enhance the sampling, angle restraints were introduced in the system to define the orientation and position of the binding ligand (via 5 restraint angles formed by two sets of atoms, each set on the ligand and dna molecule, respectively). Before the PMF, I need to estimate the free energy change when applying these angle restraints in the bound state of the complex. That’s where the problem takes place… I use 1000 kJ/mol/rad force constant to restraint these angles in the free energy calculation. Part of the .mdp file is like this:
; Neighborsearching
cutoff-scheme = Verlet
ns_type = grid ; search neighboring grid cells
nstlist = 10 ; 20 fs, largely irrelevant with Verlet scheme
rcoulomb = 1.0 ; short-range electrostatic cutoff (in nm)
rvdw = 1.0 ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
; Temperature coupling
tc-grps = system
tau_t = 1 ; time constant, in ps
ref_t = 300 ; reference temperature, one for each group, in K
; Pressure coupling is off
pcoupl = no
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Dispersion correction
DispCorr = EnerPres ; account for cut-off vdW scheme
; Velocity generation
gen_vel = yes ; Velocity generation is off
gen_temp = 300
gen_seed = -1
;Free energy control stuff
free_energy = yes
init_lambda_state = 0
;Vectors of lambda specified here
;Each combination is an index that is retrieved from init_lamba_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 21 22 23 24 25
restraint_lambdas =0.00 0.02 0.05 0.08 0.1 0.13 0.15 0.2 0.25 0.28 0.30 0.32 0.36 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1.0
nstdhdl = 100
sc-alpha = 0.5
sc-sigma = 0.3
sc-power = 1
sc-coul = no
and the angle restraints were specified in the topology file,
[angle_restraints]
; ai aj ak al type theta0A fcA mult theta0B fcB
; mixed state – need to have MULT listed twice
730 342 730 761 1 78.14 0.0 1 78.14 1000 1
761 730 761 760 1 89.53 0.0 1 89.53 1000 1
[dihedral_restraints]
; ai aj ak al type phi dphiA kfac phi dphi kfac
342 730 761 760 1 23.53 0 0 23.53 0 1000
730 761 760 759 1 -100.5 0 0 -100.5 0 1000
502 342 730 761 1 -142 0 0 99.3 0 1000
The problem is, the DNA system was severely distorted when applying these angle restraints in the free energy simulations, which is not spotted in the PMF simulations where the force constant is at its full strength (1000 kJ/mol/rad).
If I reduce the value of the force constant, say 85 kJ/mol/rad, the distortion will also reduced, is that make sense? I think is kind of weird because this structural problem only happens in the calculation of the free energy change (when the strength of the force constant gradually increase from 0 to 1)…
It would be highly appreciated if anyone can help me with that.