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
_______________________________________________
gmx-users mailing list gmx-users@gromacs.org
http://www.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at http://www.gromacs.org/search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to [EMAIL PROTECTED]
Can't post? Read http://www.gromacs.org/mailing_lists/users.php
