[gmx-users] Diagnosing an explosion

[Steven Kirk]

Hello,

I am having a problem with my BD simulation either crashing with a Range 
error or locking up the mdrun process completely.

I am running Gromacs 3.3.1 (must use this version ATM because HPC 
provider has not upgraded yet) on a quad core Intel machine, w/ Fedora 
Core 7.

I have 1000 'particles' in a cubic box (configuration here:

http://datavet.hv.se/personal/SK/publicfiles/colloid.pdb

where the particles are placed randomly with no two particle centres 
closer than 4 nm.

The particles are defined in the topology here:

; Complete topology file
;
; particles.itp

[ defaults ]
1  1  no  1.0  1.0

[ atomtypes ]
AR   28699.81134   0.00   A   1.0   1.0

[ bondtypes ]

; molecule.itp

[ nonbond_params ]
AR AR 1 1.0 1.0

[ moleculetype ]
AR 0

[ atoms ]
1 AR 1 AR AR 1 0 28699.81134

; colloid.top
; #include "particles.itp"
; #include "molecule.itp"

[ system ]
Aggregation simulation

[ molecules ]
AR 1000

The particles interact via the tabulated non-bonded potential held here:

http://datavet.hv.se/personal/SK/publicfiles/table.xvg

and the mdrun is controlled by the mdp file settings:

;       Input file
;
title               =  Aggregation
cpp                 =  /lib/cpp
constraints         =  none
integrator          =  bd
dt                  =  0.002    ; ps !
nsteps              =  1000000  ;
nstcomm             =  1
nstxout             =  10
nstvout             =  10
nstfout             =  0
nstlog              =  1
nstenergy           =  1
nstlist             =  5
ns_type             =  grid
coulombtype         = User
vdwtype             = User
rlist               =  5.036
rcoulomb            =  5.036
rvdw                =  5.036
table-extension     = 0.3
comm-mode           = None
bd_fric             = 5000
; Berendsen temperature coupling is on in two groups
Tcoupl              =  berendsen
tc-grps             =  System
tau_t               =  0.1
ref_t               =  300
; Energy monitoring
energygrps          =  AR
; Isotropic pressure coupling is now off
Pcoupl              =  no
Pcoupltype          = isotropic
tau_p               =  0.5
compressibility     =  4.5e-5
ref_p               =  1.0
; Generate velocites is on at 300 K.
gen_vel             =  yes
gen_temp            =  300.0
gen_seed            =  -1

Two strange things happen in this simulation - between the first and 
second step, the temperature approximately doubles to 600K, and remains 
approximately at this value for almost as long as the simulation runs. 
Why does the thermostat stabilize at almost exactly 2x the requested 
temperature?

The second strange event is the explosion (I have looked at the output 
from g_energy up to this point, and the total and potential energies 
drop nicely up to the point where the simulation crashes, either with a 
range error or a frozen process.

 From md.log:

            Step           Time         Lambda
           59840      119.68001        0.00000

    Energies (kJ/mol)
         LJ (SR)   Coulomb (SR)      Potential    Kinetic En.   Total Energy
    -2.18045e+04    0.00000e+00   -2.18045e+04    7.47528e+03   -1.43292e+04
     Temperature Pressure (bar)
     5.99376e+02    1.67395e-01

            Step           Time         Lambda
           59841      119.68201        0.00000

    Energies (kJ/mol)
         LJ (SR)   Coulomb (SR)      Potential    Kinetic En.   Total Energy
    -2.85992e+29    0.00000e+00   -2.85992e+29    9.18128e+25   -2.85900e+29
     Temperature Pressure (bar)
     7.36165e+24    1.50577e+21

            Step           Time         Lambda
           59842      119.68401        0.00000

    Energies (kJ/mol)
         LJ (SR)   Coulomb (SR)      Potential    Kinetic En.   Total Energy
            -inf    0.00000e+00           -inf            inf            nan
     Temperature Pressure (bar)
             inf            inf

            Step           Time         Lambda
           59843      119.68600        0.00000

    Energies (kJ/mol)
         LJ (SR)   Coulomb (SR)      Potential    Kinetic En.   Total Energy
             nan    0.00000e+00            nan            nan            nan
     Temperature Pressure (bar)
             nan            nan

and then comes the range error.

Clearly something very nasty happens between step 59840 and 59841.

Both reducing the timestep and increasing bd_fric delay the crash, but 
do not prevent it. I have inspected the trajectory visually using VMD 
and it does not look unusual. Doing an EM minimization step on the 
original coordinates before running the full BD run does not cure the 
problem either.

I would be very grateful for suggestions as to what is causing this (my 
main suspect at the moment is the tabulated potential) and what I can do 
to either fix it or diagnose the true cause.

Many thanks in advance,
Steve Kirk

-- 
Dr. Steven R. Kirk           <Steven.Kirk at hv.se, S.R.Kirk at physics.org>
Dept. of Technology, Mathematics & Computer Science  (P)+46 520 223215
University West                                      (F)+46 520 223299
P.O. Box 957 Trollhattan 461 29 SWEDEN       http://taconet.webhop.org