Dear users, I come across with an issue when I try to do free energy calculations. The issue is about the roto-translational motions of the ligand in the decoupled state. I mean the ligand doesn't stay stable in the binding pocket as in the coupled state. It seems that the restraints that are applied on the atoms of the protein and ligand under [ intermolecular_interactions ] part (an example of it is below) in the compex top file aren't sufficient to keep the ligand from repositioning/rotation with respect to the protein in the decoupled state. Even one, two and three sets of restraints couldn't solve the issue.
[ intermolecular_interactions ] [ bonds ] ; ai aj type bA kA bB kB 629 3 6 0.597 0.0 0.597 4184.0 [ angles ] ; ai aj ak type thA fcA thB fcB 281 629 3 1 37.5 0.0 37.5 41.84 629 3 21 1 121.5 0.0 121.5 41.84 [ dihedrals ] ; ai aj ak al type thA fcA thB fcB 249 281 629 3 2 -147.4 0.0 -147.4 41.84 281 629 3 21 2 -60.5 0.0 -60.5 41.84 629 3 21 16 2 -153.9 0.0 -153.9 41.84 The mdp file used for production run with GROMACS2016.2 is here. ; RUN CONTROL ;---------------------------------------------------- integrator = sd ; leap-frog integrator nsteps = 20000000 dt = 0.001 ; 2 fs comm-mode = Linear ; remove center of mass translation nstcomm = 100 ; frequency for center of mass motion removal comm-grps = System ; OUTPUT CONTROL ;---------------------------------------------------- nstxout = 20000 ; save coordinates to .trr every 100 ps nstvout = 20000 ; don't save velocities to .trr nstfout = 0 ; don't save forces to .trr nstxout-compressed = 20000 ; xtc compressed trajectory output every 2 ps compressed-x-precision = 1000 ; precision with which to write to the compressed trajectory file nstlog = 20000 ; update log file every 2 ps nstenergy = 20000 ; save energies every 2 ps nstcalcenergy = 100 ; calculate energies every 100 steps energygrps = Protein ligand ; BONDS ;---------------------------------------------------- constraint_algorithm = lincs ; holonomic constraints constraints = h-bonds ; hydrogens only are constrained lincs_iter = 2 ; accuracy of LINCS (1 is default) lincs_order = 12 ; 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 ; 'group' is default in 4.6 - from GMX5 on 'Verlet' is default ns-type = grid ; search neighboring grid cells nstlist = 10 ; 20 fs (default is 10) pbc = xyz ; 3D PBC ; ELECTROSTATICS & EWALD ;---------------------------------------------------- coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics rcoulomb = 1.2 ; 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.12 ; grid spacing for FFT ewald-rtol = 1e-6 ; relative strength of the Ewald-shifted direct potential at rcoulomb ; VAN DER WAALS ;---------------------------------------------------- vdwtype = Cut-off vdw_modifier = Potential-switch rvdw = 1.1 ; short-range van der Waals cutoff (in nm) rvdw-switch = 1.0 ; where to start switching the LJ force DispCorr = EnerPres ; apply long range dispersion corrections for Energy and Pressure ; TEMPERATURE COUPLING (SD ==> Langevin dynamics) ;---------------------------------------------------- tc_grps = Protein_ligand Water_and_ions tau_t = 1.0 1.0 ref_t = 300 300 ; PRESSURE COUPLING ;---------------------------------------------------- pcoupl = no pcoupltype = isotropic ; uniform scaling of box vectors tau_p = 1.0 ; 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 on (if gen_vel is 'yes', continuation should be 'no') gen_seed = -1 ; Use random seed gen_temp = 300 ;---------------------------------------------------- ; FREE ENERGY free-energy = yes init-lambda = 1 delta-lambda = 0 sc-alpha = 0.3 sc-power = 1 sc-sigma = 0.25 sc-coul = yes couple-moltype = ligand couple-intramol = no couple-lambda0 = vdw-q couple-lambda1 = none nstdhdl = 100 I would try to freeze the ligand in the decoupled state in the canonical ensemle with the above restrains under [ intermolecular_interactions ] but I am not sure whether that makes sense or not? Justin says ( https://mailman-1.sys.kth.se/pipermail/gromacs.org_gmx-users/2013-August/083647.html): "...Anything that is frozen, by definition, never has its position updated. Under the influence of pressure coupling, other particles around the frozen group can have their positions scaled and thus collide with the frozen group, which has remained in its original location". I think I can use the frozen ligand in both coupled and decoupled state? And I should take into consideration the effect of the frozen ligand on the free energy calculation, right? Any comment is appreciated. -- Ahmet Yildirim -- 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.
