I tried to do COM-pulling simulation (using Gromacs 5 .0.4) to pull a molecule away from sodium ion that is bound to acidic site of protein. My reference group is protein backbone. I do that in presence of solvent. The main problem is that my pulled group is not pulled in linear fashion, it rotates around the solvent droplet so when it is being pulled it shuffle around in different directions. This happens even if I changed the pull-geometry to direction using a vector. Also, I tried to change the pull -dimension from YYY to NNY and the movement was not in z- direction. The only way I can get rid of that rotation or shuffling around is do my pulling without solvent and I need to see the pulling effect in the presence of solvent. I wonder if there is a possibility that my reference group is moving in the solvent and so the spring moves with it,although the reference is supposed to be fixed.
First of all, it does not sound very surprising to me that Na+ wiggles around your protein instead of being pulled nicely in a linear manner. The amount of it might be probably controlled by decreasing the pulling rate (pull-coord1-rate). But here's a few additional thoughts you might want to take into consideration:
1. You should probably use pull-dist instead of pull-direction because you don't have a really fixed surface so the directions are ambiguous (see http://manual.gromacs.org/online/mdp_opt.html#pull).
2. Simulating a protein and an ion without solvent is pointless and scientifically rather irrelevant. You need solvent to form hydration shell aroud Na+, not to mention the protein behavior.
3. The reason why it wasn't wiggling in the dehydrated state is because you have the angular comm-mode that is suppose to prevent the rotation around cented of mass. When you add solvent, Gromacs will be trying to "unrotate" the whole system: protein+Na+solvent which won't be as effective. One way to overcome this issue is to define two groups for COM coupling: e.g. comm-grps = Protein-Na Solvent
If you still find something very unnatural in the behavior of your system, it is because you introduce an external. artificial force that drives your system out of equilibrium...
Thanks for your response
I use.distance instead of direction and maybe I was not clear but I don't pull Na, I pull a molecule that is coordinated to Na on the surface of protein so what iam hoping to see if this molecule will be pulled with or without Na
Here is my mdp file if you can provide me with more insights
Here is a copy of my .mdp file
title = Umbrella pulling simulation
define = -DPOSRES_Protein
; 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 md run
; Single-range cutoff scheme
nstlist = 50
ns_type = simple
cutoff-scheme = verlet
comm_mode = linear
coulombtype = cut-off
vdwtype = cut-off
rlist = 333.3
rcoulomb = 333.3
rvdw = 333.3
Tcoupl = nose-hoover
ref_t = 370
tau_t = 0.5
tc-grps = system
nsttcouple = 1
;Pressure coupling
Pcoupl = NO
;Generate velocities is off
gen_vel = no
; Periodic boundary conditions are on in all directions
pbc = xyz
; Pull code
pull = umbrella
pull_geometry = distance ; simple distance increase
pull_dim = Y Y Y
pull_start = yes ; define initial COM distance > 0
pull_ngroups = 2
pull_ncoords = 1
pull_group1_name = X ;its is called like that in index file
pull_group2_name = Backbone
pull_coord1_groups = 1 2
pull_coord1_rate = 0.02 ;
pull_coord1_k = 1000 ; kJ mol^-1 nm^-2
It will try suggestion 3 about comm mode for sure but I provided my mdp file just if there is something else in there
The cutoff radius for list, coulomb, and vdw is a crazy value, 333.3 nm !? This will seriously slow down your simulation, 1-2 orders of magnitude, while basically doing nothing. A typical value is 1.0 nm. Also, make sure you have long range Coulomb interactions turned on.
well I tried every changes in the.mdp file but it seems that I can not get a linear pull any further suggestions?
You could try to post your current .mdp and maybe snapshots from the trajectory, otherwise it's difficult to say.
Hi,
What is the difference between pull_geometry= distance and direction( vec= 0 0 1)? i want to pull molecule along sheet( freeze in all dimensions), pull_dim= N N Y?
Does it make any difference in PMF ?
Dear all.
I want to perform an unconventional pulling work, in which the pulling work must be stopped when the displacement of ligand's COM reaches a cut-off value.
Can we do that ???
Duc Toan Truong a good reading here: http://manual.gromacs.org/documentation/5.1/user-guide/mdp-options.html#com-pulling
Coming back to your question, combining the parameter "pull_coord1_rate" and the length of simulation gives the total displacement of the pulling.
For example, if you have:
pull_coord1_rate = 0.01 ; 0.01 nm per ps = 10 nm per ns # as in the tutorial http://www.mdtutorials.com/gmx/umbrella/05_pull.html
and a total time of simulation of 500ps then your total pulling distance will be: 0.01nm/ps* 500ps = 5 nm.
So to define a cut-off you can play with those two parameters.
Bakary Ntji Diallo
Sorry Bakary.
I think the displacement (5nm) you mention is belonged to the cantilever's movement. We can't determine that ligand'll reach the same distance.
- Badges
- Science topic