[gmx-users] looks like a bug to me ...

[m b]

dear gmx-user,

First, thanks to David for clarifying
the details about the table-file content!

Second: the problems seems to be rather involved: 

I used a shifted VdW potential (vdwtype = shift). In this case 
Gromacs uses the tabulated VdW function and NOT the built-in routine
as I would have expected. In my table file I had a value of 1.0 for all r 
for all dispersion and repulsion terms which of course gives strange 
results. (would changing the according code be easy ?)
If I use vdwtype = cut-off, the system again behaves well, both with 
the table6-12.xvg as well as with my self-made table file ... 

BUT: this still leaves the question why the system also blows up with 
the table6-12.xvg when I use (vdw-type = shift) 
actually there's more:
with the following in md.mdp :

[...]
nstlist = 1
rlist = 1.2
vdwtype = shift
rvdw = 1.2
rvdw_switch = 1.1
coulombtype = User
rcoulomb = 1.2
constraints = none
[...]

after "mdrun -s topol.tpr -table table6-12.xvg"
it says in md.log:

Table routines are used for coulomb: TRUE
Table routines are used for vdw:     TRUE
Using shifted Lennard-Jones, switch between 1.1 and 1.2 nm
Cut-off's:   NS: 1.2   Coulomb: 1.2   LJ: 1.2
Read user tables from table6-12.xvg with 1001 data points.
Tabscale = 500 points/nm
Generated table with 1000 data points for LJ6Shift.
Tabscale = 500 points/nm
Generated table with 1000 data points for LJ12Shift.
Tabscale = 500 points/nm
[...]
           Step           Time         Lambda      Annealing
              0        0.00000        0.00000        1.00000

   Energies (kJ/mol)
           Bond          Angle    Proper Dih.          LJ-14     Coulomb-14
    3.41280e+00    6.29621e+00    3.11374e+01   -7.44360e+00   -2.48470e+01
        LJ (SR)   Coulomb (SR)      Potential    Kinetic En.   Total Energy
    4.80048e+13    5.33325e+15    5.38126e+15    3.49249e+30    3.49249e+30
    Temperature Pressure (bar)
    1.75020e+31    6.75544e+29

While with "vdwtype=cut-off" the energies and tmeperature look perfectly
normal (i.e. much lower).
The strangest thing is yet to come: I get exactly the same output,
including the Coulomb (SR)-term, in md.log when I use the table6-12.xvg 
file or a self made table file with the same VdW values as in table6-12.xvg 
but with a (quite) different electrostaic potential. The ONLY difference 
I find in the two log files are the Coulomb-14 terms ... 
(again, with "vdw-type=cut-off" the Coulomb energies are quite 
different in the two cases but always look reasonable)

???

Apparently strange things are going on when one combines tables
and vdwtype=shift ... I dread going into the code myself, done so 
before, always gives me a headache ...   ;)
no offense meant, Gromacs is a marvellous program, the best MD 
package out there IMHO, and I've seen a few ...

thanks in advance for any comments/suggestions!
cheers!
Michael

ps: another (un-related) thing that looks bug-ish:
when a considerable part of the system is frozen I always
end up with ridiculously high starting temperatures (>10000 K)
when I say "gen_vel = yes" and "gem_temp = 300.0".
I get around this by simple saying something like: "gen_temp = 0.1"
but its still a bit worrying, possibly the reduced degrees of freedom 
are not accounted for consistently ...



__________________________________
Do you Yahoo!?
New Yahoo! Photos - easier uploading and sharing.
http://photos.yahoo.com/