[gmx-users] Problem in velocity plot while doing MD simulation of DNA protein with E field

[ARNAB MUKHERJEE]

Hi,

I am trying to study the motion of protein along a DNA with application of
electric field. I have position restrained the DNA, and the DNA axis is
along Z, and I am applying an electric field also along +Z. Now since the
protein is strongly positively charged it moves along the DNA in the
direction of the field (along +Z).

Now when I extract the z coordinate of the central C_alpha atom of the
protein, as expected I see for each cycle (since my box is periodic in all
directions) the z coordinate increases with time. But the problem is when I
extract and plot the z component of velocity of the central C_alpha atom of
protein, it fluctuates about 0. I don't understand how can the z component
of the velocity be even negative, when z coordinate shows positive slope
all the time.

In order to extract the position and velocities, I used the following
commands :

gmx_mpi traj -f traj_comp.xtc -s md_run.tpr -n index1.ndx -ox coor-CB.xv

gmx_mpi traj -f traj.trr -s md_run.tpr -n index1.ndx -ov vel-CB.xvg

I hope that the commands are correct.

Below I have pasted my .mdp file

title           = NVT equilibration with position restraint on all solute
(topology modified)
; Run parameters
integrator      = sd            ; leap-frog integrator
nsteps          = 2500000       ; 1 * 500000 = 500 ps
;nsteps          = 5000
dt              = 0.02          ; 1 fs
; Output control
nstxout         = 1000  ; save coordinates every 10 ps
nstvout         = 1000  ; save velocities every 10 ps
nstcalcenergy   = 50
nstenergy       = 1000  ; save energies every 1 ps
nstxtcout       = 2500
;nstxout-compressed  = 5000   ; save compressed coordinates every 1.0 ps
                             ; nstxout-compressed replaces nstxtcout
;compressed-x-grps  = System  ; replaces xtc-grps
nstlog          = 1000  ; update log file every 1 ps
; Bond parameters
continuation    = no   ; first dynamics run
constraint_algorithm = lincs ; holonomic constraints
constraints     = none          ; all bonds (even heavy atom-H bonds)
constrained
;lincs_iter     = 2                         ; accuracy of LINCS
lincs_order     = 4                         ; also related to accuracy
epsilon_r       = 15
; Neighborsearching
cutoff-scheme   = Verlet
ns_type         = grid          ; search neighboring grid cels
nstlist         = 10            ; 20 fs
rvdw_switch     = 1.0
rlist           = 1.2           ; short-range neighborlist cutoff (in nm)
rcoulomb        = 1.2           ; short-range electrostatic cutoff (in nm)
rvdw            = 1.2           ; short-range van der Waals cutoff (in nm)
vdwtype         = Cut-off       ; Twin range cut-offs rvdw >= rlist
;vdw-modifier    = Force-switch
;Electrostatics
coulombtype     = PME           ; Particle Mesh Ewald for long-range
electrostatics
pme_order       = 4                 ; cubic interpolation
fourierspacing  = 0.12          ; grid spacing for FFT

ld-seed = -1

; Temperature coupling is on
tcoupl          = v-rescale
tc_grps         = System
tau_t           = 2.0
ref_t           = 300

;energygrps = DNA Protein W ION

;freezegrps = Frozen-DNA-atoms
;freezedim = Y Y Y

E-x = 1 0 1             ;
E-y = 1 0 1             ;
E-z = 1 0.28 1         ;

; Pressure coupling is off
pcoupl          = no            ; no pressure coupling in NVT
; Periodic boundary conditions
pbc             = xyz           ; 3-D PBC
; Dispersion correctiion
DispCorr        = no            ; account for cut-off vdW scheme
; Velocity generation
gen_vel         = yes           ; assign velocities from Maxwell
distribution
gen_temp        = 300           ; temperature for Maxwell distribution
gen_seed        = -1            ; generate a random seed
; COM motion removal
; These options remove motion of the protein/bilayer relative to the
solvent/ions
nstcomm         = 50
comm-mode       = Linear
comm-grps       = DNA_W
;
refcoord_scaling = all

Earlier I had "coom-grps = System". So I was thinking that because it was
removing the velocity of COM of the system, hence the velocity of protein
that I was getting was only due to the thermal fluctuation, since it is
removing the velocity due to the electric field. Hence now I used
"coom-grps = DNA_W" where "DNA_W" is the group containing all the DNA and
water atoms, but still when I plot the z component of velocity of the
central C_alpha atom of protein, it again shows fluctuation about 0, I
don't understand why.

Can anyone please help me understand the problem here?

Thank you in advance,

Regards,

Arnab Mukherjee