[gmx-users] Diagnosing an explosion

[David van der Spoel]

Steven Kirk wrote:
> 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
Looks like you are using dimensionless eps and sigma while using 
otherwise MD units (e.g. T). Check chapter 2 of the manual. You want to 
multiply your T by boltz (0.008314 or something like that) probably.


> 
> [ 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
> 


-- 
David van der Spoel, Ph.D.
Molec. Biophys. group, Dept. of Cell & Molec. Biol., Uppsala University.
Box 596, 75124 Uppsala, Sweden. Phone:	+46184714205. Fax: +4618511755.
spoel at xray.bmc.uu.se	spoel at gromacs.org   http://folding.bmc.uu.se