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

[Steven Neumann]

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 ??

Profiles do not vary slightly - the maximum pulling force (breaking
point) varies from 290 to 500 kJ/mol nm2 which is really a lot.

Steven

>
> 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 list    gmx-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 to gmx-users-request at gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists