[gmx-users] tables don't work, more diagnostics ...

m b mic0404 at yahoo.com
Sat Dec 13 23:32:00 CET 2003


dead gxm-users,

I don't manage to get gromacs to work with
pre-calculated tables for the VdW and Coulomb potentials.
Below I describe some test-calculations I made and I'd
be very grateful if somebody from the developers group
or anybody else who's familiar with the Gromacs intestines
could have a look at it.

thanks in advance for any help!
cheers,
Michael


What I see is very strange indeed ...

I've now compared two runs with both the same starting configurations
and topologies, the ONLY difference being that in one case I use
tables (the table6-12.xvg file in the dist) for both LJ and Coulomb 
and in the other case not. (md.mdp attached below)
In both cases I cut-off both LJ and Coulomb interactions
at 1.2 nm (no shift, no lattice summation). My md.mdp file and 
further diagnostic output is attached below.

If I use tables the simulation gets stuck at t=0, it gets as far
as to print out the (huge) energy of the starting config as given
below:

           Step           Time         Lambda      Annealing
              0        0.00000        0.00000        1.00000
There are 2200 atoms in your xtc output selection
           Bond          Angle    Proper Dih.          LJ-14     Coulomb-14
    3.41280e+00    6.29621e+00    3.11374e+01   -7.47061e+00   -4.78495e+01
        LJ (SR)   Coulomb (SR)      Potential    Kinetic En.   Total Energy
    4.80048e+13    5.85068e+15    5.89869e+15    1.63775e+31    1.63775e+31
    Temperature Pressure (bar)
    8.20731e+31    3.16787e+30

Why is the temperature at t=0 so high ?? (I used gen_vel=yes and 
gen_temp=0.1 in both cases!)
If I don't use tables the simulation starts and then runs
smoothly over 1000 time-steps, the energies at t=0 are:

           Step           Time         Lambda      Annealing
              0        0.00000        0.00000        1.00000
There are 2200 atoms in your xtc output selection
   Energies (kJ/mol)
           Bond          Angle    Proper Dih.          LJ-14     Coulomb-14
    3.41280e+00    6.29621e+00    3.11374e+01   -7.47061e+00   -4.78496e+01
        LJ (SR)   Coulomb (SR)      Potential    Kinetic En.   Total Energy
   -3.82688e+02    1.96422e+01   -3.77519e+02    9.69220e-01   -3.76550e+02
    Temperature Pressure (bar)
    4.85707e+00   -7.38772e+01

to make sure that the input of the two runs was the same
I diff-ed the two resulting md.log files and the ONLY
difference (apart from times and quotes) found before the 
energy output starts is:

prompt> diff no_tables-md.log tables-md.log
76c76
<    coulombtype          = Cut-off
---
>    coulombtype          = User
79c79
<    vdwtype              = Cut-off
---
>    vdwtype              = User
152,153c152,153
< Table routines are used for coulomb: FALSE
< Table routines are used for vdw:     FALSE
---
> Table routines are used for coulomb: TRUE
> Table routines are used for vdw:     TRUE
155,159c155
< Generated table with 500 data points for COUL.
< Tabscale = 500 points/nm
< Generated table with 500 data points for LJ6.
< Tabscale = 500 points/nm
< Generated table with 500 data points for LJ12.
---
> Read user tables from table6-12.xvg with 1001 data points.

Using the -debug option to compare the ctab.xvg, etc
files I find that they are practically identical in both cases.

Seeing the huge energies in the output I first thought
that possibly the energygrp_excl directive doesn't work with
tables, since there are no bonds defined for the frozen 
graphite crystal slab in this system and hence all the bonded 
carbons would interact via their non-bonded LJ and Coulomb 
interactions ...
but if I run the same input without tables and
without the energygrp_excl directive I get energies that are also 
large but decisively different to the case where I use both
tables and energygroup exclusions and the simulation does
not (at least not immediately) blow up ...

I also compared the two mdrun.log files obtained using the
-debug flag. I did find some differences but the content of these
files is rather cryptic and hard to understand ...
I'd be grateful if somebody could help me there:
There is a number of lines where it says:

"Initiating neighbourlist type X for General with 2200 SR, 0 LR atoms"

where X is a number, these numbers differ (generally they are
larger in the case using tables)

further below there's couple of lines saying, e.g.,

"fshift after SR[   13]={-1.40713e+01, -8.79709e+01, -5.82557e+02}"

the vectors on the RHS differ and are also (considerably) larger
in the case using tables. I have no idea how to interpret this output.

a final strange thing:

I compared the energy output of the run using tables
(table6-12.xvg) with another run, also using tables, but in the  
second case the tables for the Coulomb interactions are different.
In both cases the energies (including Coulomb SR) reported by 
Gromacs at t=0 just before it crashes are EXACTLY the same 
(down to the last digit) although both md.log files contain 
the lines:
"Table routines are used for coulomb: TRUE
 Table routines are used for vdw:     TRUE
 Read user tables from table6-12.xvg with 1001 data points."
but nothing starting with "Generated table ..."
The values for the Coulomb interation given in the two
tables are QUITE different ... in other words gromacs
claims to read AND use the tables but the energy it calculates
has apparently nothing to do whith what's in these tables ...

and here comes my md.mdp:

title           = MD
cpp             = /usr/bin/cpp
integrator      = md
dt              = 0.002
nsteps          = 1000
comm_mode       = linear
nstcomm         = 0
comm_grps       = System
nstxtcout       = 1
xtc_precision   = 1000.0
xtc_grps        = System
nstxout         = 0
nstvout         = 0
nstfout         = 0
nstlog          = 1
nstenergy       = 1
energygrps      = gra hxd
freezegrps      = gra
freezedim       = Y Y Y
energygrp_excl  = gra gra
nstlist         = 1
ns_type         = grid
pbc             = xyz
rlist           = 1.2
vdwtype         = User    ; (or cut-off)
rvdw            = 1.2
DispCorr        = no
coulombtype     = User    ; (or cut-off)
rcoulomb        = 1.2
tcoupl          = Berendsen
tc_grps         = System
tau_t           = 0.1
ref_t           = 300.0
pcoupl          = no
constraints     = none
gen_vel         = yes
gen_temp        = 0.1
gen_seed        = 3234


__________________________________
Do you Yahoo!?
Protect your identity with Yahoo! Mail AddressGuard
http://antispam.yahoo.com/whatsnewfree



More information about the gromacs.org_gmx-users mailing list