[gmx-users] cannot print velocity information to trr/xtc files (Carbon Naotube inside water)

[CIHAN NADIR KAPLAN]

Dear all,

I am doing a simulation of a carbon nanotube (CNT) inside water. First I am going to describe my problem and then give info about my simulation setup (initial conditions, parameters, etc.).

My problem is as follows:

(1) I am doing a *.gro file of a *.pdb of a (16,16) CNT, without using pdb2gmx, that is with my own code. I also add manually velocity columns (vx, vy, vz) = (0, 0, 0) to the *.gro file. It works properly, since VMD and various gromacs commands can process the *.gro file.
 
(2) Then, I use editconf  in order to center the nanotube and set up a rectangular box. There, when I give a finite velocity (just to check) to every carbon (C) atom in CNT , such as (vx, vy, vz) = (0.0001, 0.0001, 0.0001) for every C, then   editconf preserves the velocity info and includen in the output *.gro file. In fact I should fix the position of the nanotube. But when (vx, vy, vz) = (0, 0, 0), then it says "cannot find velocity information" or something like that and ignores the velocity columns.

(3) So, after that, when I use genbox in order to fill the box with water, than all the velocity information disappears, no matter whether the initial velocity for C atoms are finite or simply 0. The gro file includes the coordinates of both CNT and H2O molecules, but there is no velocity info!

(4) Then I run grompp with my output *.gro file, *.top (topology) file and *.mdp file which I will describe the settings below, in order to produce a *.tpr file. When I read the resulting *.tpr file with gmxdump, again there are just coordinates, no velocity columns.

(5) I do an energy minimisation step, I check the output (both the *.trr file and *.gro file) again there is no velocity info.

(6)After that, again I modify my *.mdp file in order to do MD, I run grompp just as before, than mdrun, and there is no velocity info both in *.trr and *.xtc files, even if nvout=positive integer and gen_vel=yes  (when it is set to be"no", then again there should be velocity info in the following time steps). I also tried to put to input *.gro file velocity columns by hand, but it didn't work.

What is going on here? How can I get velocity info? Any help will be greatly appreciated. Note that, there is an external acceleration on the water molecules, where is the position of CNT is fixed. So obviously water molecules should have finite velocities. 

Best,
Nadir

P.S. :Here is a part of my mdp file for md:

cpp		= /usr/bin/cpp	; Preprocessor
define          =-DFLEXIBLE
; Parameters describing what to do, when to stop and what to save
integrator	= md		; Algorithm (steep = steepest descent minimization)
dt              = 0.002         ; time step
nsteps		= 2500000	; 5 ns
nstenergy	= 1000		; Write energies to disk every nstenergy steps
nstxout         = 5000          ; collect data every 10 ps
nstvout         = 5000
nstxtcout	= 1000		; Write coordinates to disk every nstxtcout steps
xtc_grps	= CNT SOL 	; Which coordinate group(s) to write to disk
energygrps	= CNT SOL 	; Which energy group(s) to write to disk
energygrp_excl  = CNT CNT

gen-vel         = no
gen-temp        = 300
gen-seed        = 173529


; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
compressibility = 4.5E-5
nstlist		= 15		; Frequency to update the neighbor list and long range forces
ns_type		= grid 		; Method to determine neighbor list (simple, grid)
rlist		= 1.0		; Cut-off for making neighbor list (short range forces)
coulombtype	= PME 		; Treatment of long range electrostatic interactions
rcoulomb	= 1.0		; long range electrostatic cut-off
rvdw		= 1.0		; long range Van der Waals cut-off
fourierspacing  = 0.12          ; FFT grid spacing
pme_order        = 4
optimize_fft    = yes
freezegrps      = CNT
freezedim       = Y Y Y
acc_grps	= SOL
accelerate      = 0.002 0 0   ;
constraints     = none        ; Bond types to replace by constraints
pbc		= xyz           ; Periodic Boundary Conditions (yes/no)
Pcoupl	        = parrinello-rahman
Pcoupltype      = isotropic
ref_p           = 1.0 1.0
tau_p           = 0.5
tcoupl          = nose-hoover
tc_grps         = CNT SOL
ref_t           = 300 300
tau_t           = 0.2 0.2 neighbor list (simple, grid)
rlist		= 1.0		; Cut-off for making neighbor list (short range forces)
coulombtype