[gmx-users] Expanded ensemble simulation died with fatal error: Something wrong in choosing new lambda state with a Gibbs move

Michael Shirts mrshirts at gmail.com
Wed Jul 31 23:17:04 CEST 2013 

Hi Dejun-

The basic problem is that for this particular configuration, the
current state is the only state with nonzero weight. Note that the
state with the second highest weight has weight 10^-7.  When it tries
to compare weights in single precision, it has a numerical overflow
and fails.

A few things:

1. This really should be more robust, so that it will realize it's
supposed to stay in the most likely state, since that's the only state
with nonzero weight.  I have a fix that I've been working on for
exactly this problem, but it's not quite ready yet.  Hopefully in the
next couple of days.

2. This problem is very unlikely to occur in double precision, if you
can afford the performance hit in the meantime.

3. If this is a typical average difference, exchanges will be very
unlikely.  You should probably choose your lambda intervals to be a
bit closer together at the end range.

Hopefully this will give you enough information to move forward for
the time being until a better fix is implemented.


On Wed, Jul 31, 2013 at 2:15 PM, Dejun Lin <dejun.lin at gmail.com> wrote:
> Hi all,
>
> I'm running an expanded ensemble simulation using gromacs 4.6.3 and it
> crashed with the error:
>
> Fatal error:
> Something wrong in choosing new lambda state with a Gibbs move -- probably
> underflow in weight determination.
> Denominator is:   0 1.0000002384e+00
>   i                dE        numerator          weights
>   0 -9.1451739502e+02 0.0000000000e+00 0.0000000000e+00
>   1 -9.0000128174e+02 0.0000000000e+00 0.0000000000e+00
>   2 -8.8548516846e+02 0.0000000000e+00 0.0000000000e+00
>   3 -8.7096899414e+02 0.0000000000e+00 0.0000000000e+00
>   4 -8.5645288086e+02 0.0000000000e+00 0.0000000000e+00
>   5 -8.4193676758e+02 0.0000000000e+00 0.0000000000e+00
>   6 -8.2742059326e+02 0.0000000000e+00 0.0000000000e+00
>   7 -8.1290447998e+02 0.0000000000e+00 0.0000000000e+00
>   8 -7.9838836670e+02 0.0000000000e+00 0.0000000000e+00
>   9 -7.8387219238e+02 0.0000000000e+00 0.0000000000e+00
>  10 -7.6935607910e+02 0.0000000000e+00 0.0000000000e+00
>  11 -7.5483990479e+02 0.0000000000e+00 0.0000000000e+00
>  12 -7.4032379150e+02 0.0000000000e+00 0.0000000000e+00
>  13 -7.2580767822e+02 0.0000000000e+00 0.0000000000e+00
>  14 -7.1129150391e+02 0.0000000000e+00 0.0000000000e+00
>  15 -6.9677539062e+02 0.0000000000e+00 0.0000000000e+00
>  16 -6.8225927734e+02 0.0000000000e+00 0.0000000000e+00
>  17 -6.6774316406e+02 0.0000000000e+00 0.0000000000e+00
>  18 -6.5322698975e+02 0.0000000000e+00 0.0000000000e+00
>  19 -6.3871087646e+02 0.0000000000e+00 0.0000000000e+00
>  20 -6.2419470215e+02 0.0000000000e+00 0.0000000000e+00
>  21 -6.0967858887e+02 0.0000000000e+00 0.0000000000e+00
>  22 -5.9516247559e+02 0.0000000000e+00 0.0000000000e+00
>  23 -5.8064630127e+02 0.0000000000e+00 0.0000000000e+00
>  24 -5.6613018799e+02 0.0000000000e+00 0.0000000000e+00
>  25 -5.5161407471e+02 0.0000000000e+00 0.0000000000e+00
>  26 -5.3709790039e+02 0.0000000000e+00 0.0000000000e+00
>  27 -5.2258178711e+02 0.0000000000e+00 0.0000000000e+00
>  28 -5.0806564331e+02 0.0000000000e+00 0.0000000000e+00
>  29 -4.9354953003e+02 0.0000000000e+00 0.0000000000e+00
>  30 -4.7903335571e+02 0.0000000000e+00 0.0000000000e+00
>  31 -4.6451724243e+02 0.0000000000e+00 0.0000000000e+00
>  32 -4.5000018311e+02 0.0000000000e+00 0.0000000000e+00
>  33 -4.3548400879e+02 0.0000000000e+00 0.0000000000e+00
>  34 -4.2096792603e+02 0.0000000000e+00 0.0000000000e+00
>  35 -4.0645178223e+02 0.0000000000e+00 0.0000000000e+00
>  36 -3.9193563843e+02 0.0000000000e+00 0.0000000000e+00
>  37 -3.8107025146e+02 0.0000000000e+00 0.0000000000e+00
>  38 -3.6290338135e+02 0.0000000000e+00 0.0000000000e+00
>  39 -3.4838723755e+02 0.0000000000e+00 0.0000000000e+00
>  40 -3.3387109375e+02 0.0000000000e+00 0.0000000000e+00
>  41 -3.1935494995e+02 0.0000000000e+00 0.0000000000e+00
>  42 -3.0483883667e+02 0.0000000000e+00 0.0000000000e+00
>  43 -2.9032269287e+02 0.0000000000e+00 0.0000000000e+00
>  44 -2.7580654907e+02 0.0000000000e+00 0.0000000000e+00
>  45 -2.6129040527e+02 0.0000000000e+00 0.0000000000e+00
>  46 -2.4677430725e+02 0.0000000000e+00 0.0000000000e+00
>  47 -2.3225816345e+02 0.0000000000e+00 0.0000000000e+00
>  48 -2.1774200439e+02 0.0000000000e+00 0.0000000000e+00
>  49 -2.0322586060e+02 0.0000000000e+00 0.0000000000e+00
>  50 -1.8970976257e+02 0.0000000000e+00-1.0000000000e+00
>  51 -1.7419361877e+02 0.0000000000e+00 0.0000000000e+00
>  52 -1.5967747498e+02 0.0000000000e+00 0.0000000000e+00
>  53 -1.4516131592e+02 0.0000000000e+00 0.0000000000e+00
>  54 -1.3064523315e+02 0.0000000000e+00 0.0000000000e+00
>  55 -1.1612908173e+02 0.0000000000e+00 0.0000000000e+00
>  56 -1.0161293030e+02 7.0064923216e-45 0.0000000000e+00
>  57 -8.7096786499e+01 1.4939846888e-38 0.0000000000e+00
>  58 -7.2580688477e+01 3.0102835162e-32 0.0000000000e+00
>  59 -5.8064540863e+01 6.0658294505e-26 0.0000000000e+00
>  60 -4.3548393250e+01 1.2222865341e-19 0.0000000000e+00
>  61 -2.9032243729e+01 2.4629560544e-13 0.0000000000e+00
>  62 -1.5516148567e+01 1.8256696421e-07-1.0000000000e+00
>  63  0.0000000000e+00 9.9999976158e-01-1.0000000000e+00
>
> The mdp options for the free energy and expanded ensemble stuff are:
>
> free-energy              = expanded
>
> ; no need to mess with these for now
> ;--------
> sc-alpha                 = 0
> sc-power                 = 0
> sc-r-power               = 6
> sc-coul                  = no
> -------
>
> ; Which intermediate state are we simulating?
> -------
> init-lambda-state        = 0
>
> ; What are the values of lambda at the intermediate states?
> ;-------
> ; fep-lambdas             = 0.0 0.06667 0.1333 0.2 0.2667 0.3333 0.4 0.4667
> 0.5333 0.6 0.6667 0.7333 0.8 0.8667 0.9333 1.0
> bonded-lambdas           =  0.000000 0.000000 0.000000 0.000000 0.000000
> 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
> 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
> 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
> 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
> 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
> 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
> 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
> 0.000000 0.000000 0.000000
> vdw-lambdas              =  0.000000 0.015873 0.031746 0.047619 0.063492
> 0.079365 0.095238 0.111111 0.126984 0.142857 0.158730 0.174603 0.190476
> 0.206349 0.222222 0.238095 0.253968 0.269841 0.285714 0.301587 0.317460
> 0.333333 0.349206 0.365079 0.380952 0.396825 0.412698 0.428571 0.444444
> 0.460317 0.476190 0.492063 0.507937 0.523810 0.539683 0.555556 0.571429
> 0.58331 0.603175 0.619048 0.634921 0.650794 0.666667 0.682540 0.698413
> 0.714286 0.730159 0.746032 0.761905 0.777778 0.793651 0.809524 0.825397
> 0.841270 0.857143 0.873016 0.888889 0.904762 0.920635 0.936508 0.952381
> 0.968254 0.984127 1.000000
> mass-lambdas             =  0.000000 0.000000 0.000000 0.000000 0.000000
> 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
> 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
> 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
> 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
> 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
> 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
> 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
> 0.000000 0.000000 0.000000
> coul-lambdas             =  0.000000 0.000000 0.000000 0.000000 0.000000
> 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
> 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
> 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
> 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
> 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
> 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
> 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
> 0.000000 0.000000 0.000000
> restraint-lambdas                =  0.000000 0.000000 0.000000 0.000000
> 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
> 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
> 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
> 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
> 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
> 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
> 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
> 0.000000 0.000000 0.000000 0.000000
>
> ; This makes sure we print out the differences in Hamiltonians between all
> states, and not just the neighboring states
> ;--------
> calc-lambda-neighbors = -1
>
> ; the frequency the free energy information is calculated.  This
> ; frequency (every 0.2 ps) is pretty good for small molecule solvation.
> ;-------
> nstdhdl                  = 1000
> ; not required, but useful if you are doing any temperature reweighting.
> Without
> ; temperature reweighting, you don't need the total energy -- differences
> are enough
> dhdl-print-energy        = yes
>
> ; We are doing free energies with the LIPO_MUT molecule alone
> couple-moltype           = LIPO_MUT
> ; we are mutating on type of molecule into another.  In the initial state,
> both are on
> couple-lambda0           = vdw-q
> ; in the final state, both are on.
> couple-lambda1           = vdw-q
> ; let the intramolecular interaction be coupled too
> couple-intramol          = yes
>
> ; expanded ensemble stuff
> nstexpanded              = 100
> ; Wang-Landau algorithm to determine the free energies 'weights' of the
> states
> lmc-stats                = wang-landau
> ; Metropolized gibbs algorithm to move between states
> lmc-move                 = metropolized-gibbs
> ; we stop equilibrating when the wang-landau scaling term gets as low as
> 0.0001
> lmc-seed                 = 7890
> lmc-weights-equil        = wl-delta
> weight-equil-wl-delta    = 0.0001
>
> ; Seed for Monte Carlo in lambda space
> ; We scale our wang landau weight by 0.7, whenever the smallest state
> ; and largest state have ratio of 0.8.  The initial wang-landau weight
> ; increment delta is 1 kbT, and when this delta<1/N, where N is the
> ; number of attempted switches in state space, we use 1/N as the delta,
> ; which is less prone to saturation (stopping at the wrong value because
> ; the weight schedule lowered too quickly).
> wl-scale                 = 0.7
> wl-ratio                 = 0.8
> init-wl-delta            = 1 ; this is 1*kB*T
> wl-oneovert              = yes
>
> ; frequency to output transition matrix
> nst-transition-matrix    = 10000000
>
> Any idea?
>
> Thanks,
> Dejun
> --
> gmx-users mailing list    gmx-users at gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-request at gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

More information about the gromacs.org_gmx-users mailing list