[gmx-users] Problem with Position Restraints (before running SMD)!

[Justin A. Lemkul]

Talal E. AlOtaibi wrote:
> Hi guys,
>  
> When I did my SMD simulation by applying load to a certain group of 
> atoms and fixing certain atoms by modifying the topology file, the 
> points that i have fixed moved from there original position. I don't 
> know what the problem is!
>  
> Any suggestion!
>  

Position restraints disfavor motion, they don't prevent it.  If the pulling 
force is sufficient to overcome the restraint bias, then the atoms move.  Two 
options:

1. Use a stronger position restraint force constant
2. Use freezegrps instead - these don't move

-Justin

>  
> To make what I was doing clear, I started the equilibrium by creating a 
> box and filling it with water. Then I did the energy minimization. After 
> that I heated up and pressurized the structure. Every thing is good so 
> far. After that, I modified the topology file by adding four points that 
> i want to fix as following:
> ; Position restraint for each water oxygen
> [ position_restraints ]
> ; i funct       fcx        fcy        fcz
>      1 1 1000 1000 1000
> 1330 1 1000 1000 1000              <---------   my 1st point
> 1426 1 1000 1000 1000              <---------   my 2nd point
> 1437 1 1000 1000 1000              <---------   my 3rd point
> 1548 1 1000 1000 1000              <---------   my 4th point
> #endif
> This only what I did to fix the points.
>  
> The results have shown that the fixid points moved from their original 
> position.
>  
>  
>  
> Here is my topology file:
>  
> ;
> ;       File '3gm1a.top' was generated
> ;       By user: onbekend (0)
> ;       On host: onbekend
> ;       At date: Sat Jul  2 14:03:14 2011
> ;
> ;       This is a standalone topology file
> ;
> ;       It was generated using program:
> ;       pdb2gmx - VERSION 4.5.3
> ;
> ;       Command line was:
> ;       pdb2gmx -f 3gm1a.pdb -o 3gm1a.gro -p 3gm1a.top
> ;
> ;       Force field data was read from:
> ;       /opt/apps/pgi7_2/mvapich1_1_0_1/gromacs/4.5.3/share/gromacs/top
> ;
> ;       Note:
> ;       This might be a non-standard force field location. When you use 
> this topology, the
> ;       force field must either be present in the current directory, or 
> the location
> ;       specified in the GMXLIB path variable or with the 'include' mdp 
> file option.
> ;
>  
> ; Include forcefield parameters
> #include "gromos53a6.ff/forcefield.itp"
>  
> ; Include chain topologies
> #include "3gm1a_Protein_chain_A.itp"
> #include "3gm1a_Protein_chain_E.itp"
> #include "3gm1a_Protein_chain_F.itp"
>  
> ; Include water topology
> #include "gromos53a6.ff/spc.itp"
> #ifdef POSRES_WATER
> 
> ; Position restraint for each water oxygen
> [ position_restraints ]
> ;  i funct       fcx        fcy        fcz
>      1    1       1000       1000       1000
> 1330    1       1000       1000       1000
> 1426    1       1000       1000       1000
> 1437    1       1000       1000       1000
> 1548    1       1000       1000       1000
> 
> #endif
>  
> ; Include topology for ions
> #include "gromos53a6.ff/ions.itp"
>  
> [ system ]
> ; Name
> PROTEIN TYROSINE KINASE 2 BETA; 5 RESIDUES 861-1009; PAXILLIN in water
>  
> [ molecules ]
> ; Compound        #mols
> Protein_chain_A     1
> Protein_chain_E     1
> Protein_chain_F     1
> SOL         8788
> NA          25
> CL          18
> 
>  
> And Here is my mdp file:
>  
> title           = smd_120 ; 120 pN
> ; this is loosely based off of the VT pulling tutorial; heavily modified
> ; Run parameters
> integrator      = md
> dt              = 0.002
> tinit           = 0
> nsteps          = 2500000       ; 5 ns
> ; Output parameters
> nstxout         = 1000          ; every  2 ps
> nstvout         = 1000
> nstfout         = 5000
> nstxtcout       = 5000          ; every ps
> nstenergy       = 1000
> ; Bond parameters
> constraint_algorithm    = lincs
> constraints             = hbonds
> lincs_iter      = 1              ; accuracy of LINCS
> lincs_order     = 4              ; also related to accuracy
> ; 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.16
> 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.2 0.2
> 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
> ; COM motion removal
> ; These options remove comm motion of the constraint / freeze group
> nstcomm         = 1
> comm_mode       = Linear
> comm_grps       = System
> ; pull parameters
> pull            = constant_force
> pull_geometry   = direction
> pull_nstxout    = 500 ; will print the c.o.m. coordinates
> pull_nstfout    = 500 ; forces on group
> pull_ngroups    = 1
> pull_group0     = Protein ;
> pull_group1     = group_B ;
> pull_pbcatom1   = 100 ; the CA of the 10th residue
> pull_vec1       = -2.63 -16.48 -14.95 ; direction of pull, will be 
> normalized
> pull_k1         = 72.29 ; pull_k1*1.66 = pN;  units: [kJ / (mol * nm^2)]
>  
>  
>  
>  
>  
> What should I do?
>  
>  
> Thanks,
> Talal
>  
> 

-- 
========================================

Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin

========================================