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

[Nadir Kaplan]

Dear Mark,


>>
>> 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.
>
> This doesn't serve a purpose that I can see.
>
>> (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!
>
> Keeping them really doesn't serve a purpose - you've radically changed 
> the system, so those old velocities can't be from a meaningful 
> ensemble any more, even if they were to start with.



I see.



>
>> (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.
>
> With gen_vel = no then this should be expected.
>
>> (5) I do an energy minimisation step, I check the output (both the
>> *.trr file and *.gro file) again there is no velocity info.
>
> Since EM is mechanical, not dynamical, there are no velocities.
>
>> (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). 
>
> You need 0 < *nstvout* to see velocity output. You won't see 
> velocities in the .xtc file because that's not the purpose of the .xtc 
> file. See http://wiki.gromacs.org/index.php/.xtc_file. IIRC the value 
> of gen_vel is immaterial here, since the dynamics will start from zero 
> velocities if you don't generate them. grompp will probably give you a 
> warning about this, so make sure you check for them! mdrun will always 
> write velocities at the final timestep, regardless of nstvout too. You 
> will also need to change the integrator to md.



OK, I see. But still there is no velocity output in the trr files. I
checked the grompp output, it is as follows:

checking input for internal consistency...
calling /usr/bin/cpp...
processing topology...
Generated 165 of the 1596 non-bonded parameter combinations
Excluding 3 bonded neighbours for CNT 1
Excluding 2 bonded neighbours for SOL 8962
processing coordinates...
double-checking input for internal consistency...
renumbering atomtypes...
converting bonded parameters...
#      BONDS:   17924
#     ANGLES:   8962
initialising group options...
processing index file...
Analysing residue names:
Opening library file /usr/share/gromacs/top/aminoacids.dat
There are:  8963      OTHER residues
There are:     0    PROTEIN residues
There are:     0        DNA residues
Analysing Other...
Making dummy/rest group for Acceleration containing 1280 elements
Making dummy/rest group for Freeze containing 26886 elements
Making dummy/rest group for VCM containing 28166 elements
Number of degrees of freedom in T-Coupling group CNT is 0.00
Number of degrees of freedom in T-Coupling group SOL is 80655.00
Making dummy/rest group for User1 containing 28166 elements
Making dummy/rest group for User2 containing 28166 elements
Making dummy/rest group for Or. Res. Fit containing 28166 elements
Making dummy/rest group for QMMM containing 28166 elements
T-Coupling       has 2 element(s): CNT SOL
Energy Mon.      has 2 element(s): CNT SOL
Acceleration     has 2 element(s): SOL rest
Freeze           has 2 element(s): CNT rest
User1            has 1 element(s): rest
User2            has 1 element(s): rest
VCM              has 1 element(s): rest
XTC              has 1 element(s): CNT
Or. Res. Fit     has 1 element(s): rest
QMMM             has 1 element(s): rest
WARNING 1 [file aminoacids.dat, line 1]:
   Can not exclude the lattice Coulomb energy between energy groups
Checking consistency between energy and charge groups...
Net Acceleration in X direction, will  be corrected
Calculating fourier grid dimensions for X Y Z
Using a fourier grid of 70x50x50, spacing 0.114 0.120 0.120
writing run input file...
There was 1 warning


>
> > I also tried to put to
> > input *.gro file velocity columns by hand, but it didn't work.
>
> That doesn't make sense in the context.
>
>> 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.
>
> Turn off the weird stuff (freeze groups and accelerate groups) until 
> you've got your basic issues sorted out. MD simulations in general and 
> CNT simulations in particular are not easy things to do right. Learn 
> to walk before you try to run! :-)



Freezing is for simplicity. One needs velocity to run.

Below I paste another mdp file which I used for the MD part of the same
simulation. That MD simulation was succesfully over and I was able to
reproduce some results in the literature by using position info. Again,
in the .trr file of that simulation, there are no velocities.

LINES STARTING WITH ';' ARE COMMENTS
title        = Nanotube in Water (MD)    ; Title of run

; The following lines tell the program the standard locations where to
find certain files
cpp        = /usr/bin/cpp    ; Preprocessor
define          =-DFLEXIBLE
; Parameters describing what to do, when to stop and what to save
integrator    = md        ; Algorithm (md = molecular dynamics)
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


; 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