[gmx-users] Re: diagnosing an explosion
Steven Kirk
Steven.Kirk at hv.se
Mon Dec 10 16:10:01 CET 2007
Reply below
>
>> 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.
Aha! Thank you David.
My goal was to use MD units in everything, since my tabulated data file
was extracted from another GROMACS simulation and is hence already in MD
units (kJ/mol vs. nm), and I would prefer to specify temperatures in K.
My reason for specifying 1.0 for the last two parameters was that I
thought this would allow me to use, unscaled, the tabulated potential
(which is already in MD units).
If I want to do this, what alternative values should I use for the last
two parameters on the [ atomtypes ] line above (given that I am really
only using the g(r) and g''(r) columns to hold my potential data: the
h(r) and h''(r), which I don't care about, are all zeroes) ?
I have checked Chapter 2 in the manual and I am still unsure how to
handle this.
I took your advice (which matches with manual Section 2.4 on *reduced*
units) and scaled my temperature to 2.494353 = 300K * 0.008314. This
did not give a crash, but increasing the timestep and decreasing bd_fric
made the problem return.
Scaling the temperature value in this way did not fix the 'equilibration
around 2x requested temperature' problem mentioned below - g_energy
still shows the temperature settling at a numerical value around 4.988.
Many thanks for any further suggestions!
>>
>> [ 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
More information about the gromacs.org_gmx-users
mailing list