[gmx-users] Umbrella sampling - chain A does not get pulled

[Davide Mercadante]

Dear Justin, 

Thank you for your reply. To be sure I carefully ran the simulations again
using the following .mdp file that includes the pull lines and
-DPOSRES_B to restrain chain B from being pulled when chain A is being
pulled (it is the .mdp file given in the tutorial):

title       = Umbrella pulling simulation
define      = -DPOSRES_B
; Run parameters
integrator  = md
dt          = 0.002
tinit       = 0
nsteps      = 250000    ; 500 ps
nstcomm     = 10
; Output parameters
nstxout     = 5000      ; every 10 ps
nstvout     = 5000 
nstfout     = 500
nstxtcout   = 500       ; 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.4
; 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
; Berendsen temperature coupling is on in two groups
Tcoupl      = Nose-Hoover
tc_grps     = Protein   Non-Protein
tau_t       = 0.5       0.5
ref_t       = 310       310
; Pressure coupling is on
Pcoupl          = Parrinello-Rahman
pcoupltype      = isotropic
tau_p           = 1.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     = Chain_B
pull_group1     = Chain_A
pull_rate1      = 0.01      ; 0.01 nm per ps = 10 nm per ns
pull_k1         = 1000      ; kJ mol^-1 nm^-2


I get slightly different distances but it is still now working. From the
load of the trajectories chain B is not getting pulled. However, after
running the distance.pl script I get a COM distance within the range 2.3nm
to 6.0nm). Can you please help me to understand what I am doing wrong?

My command line is:

Grompp -f pull.mdp -c after_npt.gro -p topol.top -o pull.tpr -t state.cpt
(checkpoint file from the previous npt step) -n index.ndx (which includes
groups 19 and 20, chains A and B respectively).

Mdrun then uses pull.tpr...

Please don't hesitate to ask me any other information you may need to
figure this out. I am not sure what I am doing wrong..

Thank you for your time.

Davide

On 12/12/12 2:38 AM, "Justin Lemkul" <jalemkul at vt.edu> wrote:

>
>
>On 12/11/12 2:40 AM, Davide Mercadante wrote:
>> Dear Justin,
>>
>> I have been practicing umbrella sampling simulations following your
>>tutorial
>> step by step. I have just finished to perform the pull simulations to
>> identify the configurations to use in the umbrella runs. I have used the
>> distances.pl script to run iteratively g_dist and the resulting file
>> (summary_distances.dat) is showing a distance between 2.62nm (for
>>conf0.gro)
>> and 5.3nm (for conf500.gro). I don't seem to have the values reported
>>in the
>> tutorial at all (in the tutorial is reported that the COM distance
>>would be
>> something of an order of magnitude lowerŠ) and I am not sure if
>>something is
>> gone wrong.
>>
>
>The distances shown in the tutorial are just an illustrative example and
>don't 
>necessarily correspond to anything you might obtain.  A distance of 2.62
>nm for 
>the initial configuration does not seem right to me though; it should be
>on the 
>order of 0.5 nm, coincident with interpeptide spacing in the beta-sheet.
>
>> I have also tried to visualize the trajectories but I am not able to see
>> chain_A moving away from the rest of the fibril. I have set, at the end
>>of
>> the topol_Protein_chain_B.itp file, the lines suggested to restrain
>>chain B.
>>
>> #ifdef POSRES_B
>> #include "posre_Protein_chain_B.itp"
>> #endif
>>
>> I have used all the input files given in the tutorial and the
>>commands/files
>> extensions suggested.
>>
>> Can you please help me to understand if I am doing something wrong?
>>
>
>What you're reporting is not consistent - your distances increase over
>time, but 
>then you say peptide A does not move away from the others?  That's not
>possible. 
>  If peptide A doesn't move, then the COM distance should remain
>relatively 
>constant.  Are you loading the right trajectory (i.e., the actual SMD,
>not a 
>previous equilibration or something else)?
>
>-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
>* 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