[gmx-users] Pulling ligand - Different Profiles (Force vs time)

[Thomas Schlesier]

As a side note:

The rupture process is a stochastic process, so a single rupture force 
is meaningless, since it is a distributed property. So you need to do 
many simulations to get the distribution / average rupture force.
It that same like equilibrium properties, one doesn't determine them 
from only a single frame from a trajectory.


Am 27.06.2012 15:52, schrieb gmx-users-request at gromacs.org:
> On 6/27/12 9:36 AM, Steven Neumann wrote:
>> >  On Wed, Jun 27, 2012 at 1:51 PM, Justin A. Lemkul<jalemkul at vt.edu>  wrote:
>>> >>
>>> >>
>>> >>  On 6/27/12 7:48 AM, Steven Neumann wrote:
>>>> >>>
>>>> >>>  Dear Gmx Users,
>>>> >>>
>>>> >>>  I obtained a protein-ligand complex from 100ns simulation. Now I am
>>>> >>>  pulling my ligand away from the protein after the energy minimzation
>>>> >>>  in water and equilibration of 100ps (two coupling baths: Protein,
>>>> >>>  LIG_Water_and_ions).
>>>> >>>  Then I proceed my pulling :
>>>> >>>
>>>> >>>  grompp -f pull.mdp -c npt.gro -p topol.top -n index.ndx -t npt.cpt -o
>>>> >>>  pull.tpr
>>>> >>>
>>>> >>>  mdrun -s pull.tpr -deffnm pull
>>>> >>>
>>>> >>>
>>>> >>>  title       = Umbrella pulling simulation
>>>> >>>  define      = -DPOSRES
>>>> >>>  ; Run parameters
>>>> >>>  integrator  = md
>>>> >>>  dt          = 0.002
>>>> >>>  tinit       = 0
>>>> >>>  nsteps      = 500000    ; 1 ns
>>>> >>>  nstcomm     = 10
>>>> >>>  ; Output parameters
>>>> >>>  nstxout     = 0
>>>> >>>  nstvout     = 0
>>>> >>>  nstfout     = 500
>>>> >>>  nstxtcout   = 1000       ; every 1 ps
>>>> >>>  nstenergy   = 500
>>>> >>>  ; Bond parameters
>>>> >>>  constraint_algorithm    = lincs
>>>> >>>  constraints             = all-bonds
>>>> >>>  continuation            = yes       ; continuing from NPT
>>>> >>>  ; Single-range cutoff scheme
>>>> >>>  nstlist     = 5
>>>> >>>  ns_type     = grid
>>>> >>>  rlist       = 1.4
>>>> >>>  rcoulomb    = 1.4
>>>> >>>  rvdw        = 1.2
>>>> >>>  vdwtype     = Switch
>>>> >>>  rvdw-switch = 1.0
>>>> >>>  ; PME electrostatics parameters
>>>> >>>  coulombtype     = PME
>>>> >>>  fourierspacing  = 0.12
>>>> >>>  fourier_nx      = 0
>>>> >>>  fourier_ny      = 0
>>>> >>>  fourier_nz      = 0
>>>> >>>  pme_order       = 4
>>>> >>>  ewald_rtol      = 1e-5
>>>> >>>  optimize_fft    = yes
>>>> >>>  ; Temperature coupling is on
>>>> >>>  tcoupl      = V-rescale                                  ; modified
>>>> >>>  Berendsen thermostat
>>>> >>>  tc_grps     = Protein LIG_Water_and_ions   ; two coupling groups - more
>>>> >>>  accurate
>>>> >>>  tau_t       = 0.1   0.1                                 ; time constant,
>>>> >>>  in ps
>>>> >>>  ref_t       = 298   298                                  ; reference
>>>> >>>  temperature, one for each group, in K
>>>> >>>  ; Pressure coupling is on
>>>> >>>  Pcoupl          = Parrinello-Rahman
>>>> >>>  pcoupltype      = isotropic
>>>> >>>  tau_p           = 2.0
>>>> >>>  compressibility = 4.5e-5
>>>> >>>  ref_p           = 1.0
>>>> >>>  ; Generate velocities is off
>>>> >>>  gen_vel     = no
>>>> >>>  ; Periodic boundary conditions are on in all directions
>>>> >>>  pbc     = xyz
>>>> >>>  ; Long-range dispersion correction
>>>> >>>  DispCorr    = EnerPres
>>>> >>>  ; Pull code
>>>> >>>  pull            = umbrella
>>>> >>>  pull_geometry   = distance  ; simple distance increase
>>>> >>>  pull_dim        = N N Y
>>>> >>>  pull_start      = yes       ; define initial COM distance>  0
>>>> >>>  pull_ngroups    = 1
>>>> >>>  pull_group0     = Protein
>>>> >>>  pull_group1     = LIG
>>>> >>>  pull_rate1      = 0.004      ; 0.004 nm per ps = 4 nm per ns
>>>> >>>  pull_k1         = 500      ; kJ mol^-1  nm^-2
>>>> >>>
>>>> >>>  I run 3 pulling simulations with the same mdp  and I obtain 3
>>>> >>>  different profiles (Force vs time). Then I used 2xlonger pulling with
>>>> >>>  the same pulling distance and I run 3 simulations again. Each time I
>>>> >>>  obtain different profile. Can anyone explain me this? I am using
>>>> >>>  velocities from npt simulation as above (gen_vel = no and continuation
>>>> >>>  = yes) so I presume the output should be similar. Please, advice.
>>>> >>>
>>> >>
>>> >>  I assume you're passing a checkpoint file to grompp?  If you're relying on
>>> >>  velocities from the .gro file, they are of insufficient precision to
>>> >>  guarantee proper continuation.
>> >
>> >  Thank you Justin. I am using according to your tutorial:
>> >
>> >  grompp -f pull.mdp -c npt.gro -p topol.top -n index.ndx -t npt.cpt -o pull.tpr
>> >  mdrun -s pull.tpr -deffnm pull
>> >
>> >  Would you suggest:
>> >
>> >  grompp -f pull.mdp -c npt.gro -p topol.top -n index.ndx -t npt.cpt -o pull.tpr
>> >  mdrun -s pull.tpr -cpi npt.cpt -deffnm pull ??
>> >
> No, I would not, especially if the NPT run uses position restraints, in which
> case the two phases are different.  I missed the command line in the earlier
> message.  What you are doing makes sense.
>
>> >  Profiles do not vary slightly - the maximum pulling force (breaking
>> >  point) varies from 290 to 500 kJ/mol nm2 which is really a lot.
>> >
> Consult the points below and watch your trajectories.  If you're getting
> different forces, your ligands are experiencing different interactions.  SMD is
> a path-dependent, non-equilibrium process.  Good sampling and a justifiable path
> are key.
>
> -Justin
>
>>> >>
>>> >>  Small variations are inherent in any simulation set, and in the case of
>>> >>  pulling, small changes (though intentional) are the basis for Jarzynski's
>>> >>  method.  In any case, all MD simulations are chaotic and so it depends on
>>> >>  what your definition of "different" is in the context of whether or not
>>> >>  there are meaningful changes imparted through the course of each simulation.
>>> >>     Also note that in the absence of the -reprod flag, the same .tpr file may
>>> >>  result in a slightly different outcome.  The implications of these outcomes
>>> >>  are limited by sampling; the ensemble should converge with sufficient time
>>> >>  and/or replicates.  For non-equilibrium processes like pulling, convergence
>>> >>  is probably harder, but again you have to ask whether the differences are
>>> >>  meaningful.
>>> >>
>>> >>  http://www.gromacs.org/Documentation/Terminology/Reproducibility
>>> >>
>>> >>  -Justin
>>> >>
>>> >>  --
>>> >>  ========================================
>>> >>
>>> >>  Justin A. Lemkul, Ph.D.
>>> >>  Research Scientist
>>> >>  Department of Biochemistry
>>> >>  Virginia Tech
>>> >>  Blacksburg, VA
>>> >>  jalemkul[at]vt.edu | (540) 231-9080
>>> >>  http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
>>> >>
>>> >>  ========================================
>>> >>
>>> >>
>>> >>  --
>>> >>  gmx-users mailing listgmx-users at gromacs.org
>>> >>  http://lists.gromacs.org/mailman/listinfo/gmx-users
>>> >>  * Only plain text messages are allowed!
>>> >>  * Please search the archive at
>>> >>  http://www.gromacs.org/Support/Mailing_Lists/Search  before posting!
>>> >>  * Please don't post (un)subscribe requests to the list. Use the www
>>> >>  interface or send it togmx-users-request at gromacs.org.
>>> >>  * Can't post? Readhttp://www.gromacs.org/Support/Mailing_Lists
>> >
> --