The DFT solver, as written
ge_it_check.viv
Generated by RUNNING the independent numpy reference at the deck's settings and pinning what it printed; then the published LDA values; then the figure.
# Ge: what the calculation must reproduce, and the picture of it. # # dft/run.sh ge it # # GENERATED by dft/ref/make_deck.py, which RAN dft/ref/pwdft_ref.py -- an independent numpy implementation that diagonalizes H densely -- # at exactly this deck's settings and pinned what it printed. # A disagreement here is a bug in one of the two implementations. GAMMA_P = NSEG LAST = NPTS - 1.0 test "SCF converged: density, energy, and every occupied band's residual" { assert_lt(DRHO, 1.0e-8) assert_lt(EIGRES, EIG_TOL) assert_lt(ITER, MAXIT) } test "every band-path point converged" { assert_lt(BANDRES, EIG_TOL) } test "total energy agrees with the numpy reference" { assert_near(ETOT, 0.090394639792, 1.0e-7) } test "band energies at L, Gamma and X agree with the numpy reference (eV)" { assert_near(EVB[0.0, 0.0], -10.712318984, 1.0e-5) # L, band 1 assert_near(EVB[0.0, 1.0], -7.622551135, 1.0e-5) # L, band 2 assert_near(EVB[0.0, 2.0], -1.388442563, 1.0e-5) # L, band 3 assert_near(EVB[0.0, 3.0], -1.388442563, 1.0e-5) # L, band 4 assert_near(EVB[0.0, 4.0], 0.129280654, 1.0e-5) # L, band 5 assert_near(EVB[0.0, 5.0], 3.748442984, 1.0e-5) # L, band 6 assert_near(EVB[GAMMA_P, 0.0], -12.780989557, 1.0e-5) # Gamma, band 1 assert_near(EVB[GAMMA_P, 1.0], -0.015418954, 1.0e-5) # Gamma, band 2 assert_near(EVB[GAMMA_P, 2.0], -0.000000000, 1.0e-5) # Gamma, band 3 assert_near(EVB[GAMMA_P, 3.0], 0.000000000, 1.0e-5) # Gamma, band 4 assert_near(EVB[GAMMA_P, 4.0], 0.000001415, 1.0e-5) # Gamma, band 5 assert_near(EVB[GAMMA_P, 5.0], 2.580416300, 1.0e-5) # Gamma, band 6 assert_near(EVB[LAST, 0.0], -8.918309515, 1.0e-5) # X, band 1 assert_near(EVB[LAST, 1.0], -8.918309283, 1.0e-5) # X, band 2 assert_near(EVB[LAST, 2.0], -3.052852834, 1.0e-5) # X, band 3 assert_near(EVB[LAST, 3.0], -3.052852331, 1.0e-5) # X, band 4 assert_near(EVB[LAST, 4.0], 0.682617968, 1.0e-5) # X, band 5 assert_near(EVB[LAST, 5.0], 0.682618117, 1.0e-5) # X, band 6 } #: REFPATH the reference's lowest 6 bands at all 61 path points, eV REFPATH = [[-10.7123190, -7.6225511, -1.3884426, -1.3884426, 0.1292807, 3.7484430], [-10.7246885, -7.6027757, -1.3859669, -1.3859669, 0.1318768, 3.7498587], [-10.7611510, -7.5453779, -1.3788391, -1.3788391, 0.1390682, 3.7530968], [-10.8189922, -7.4526601, -1.3669690, -1.3669690, 0.1510419, 3.7582248], [-10.8945919, -7.3284312, -1.3503838, -1.3503838, 0.1677693, 3.7648262], [-10.9839473, -7.1764360, -1.3290923, -1.3290923, 0.1892190, 3.7723144], [-11.0832327, -7.0005838, -1.3029115, -1.3029115, 0.2153397, 3.7799459], [-11.1893584, -6.8042939, -1.2723530, -1.2723530, 0.2460551, 3.7867146], [-11.2995017, -6.5903512, -1.2372392, -1.2372392, 0.2812770, 3.7915185], [-11.4114609, -6.3612528, -1.1976786, -1.1976786, 0.3206571, 3.7929759], [-11.5234826, -6.1185067, -1.1536985, -1.1536985, 0.3644845, 3.7898278], [-11.6341863, -5.8638274, -1.1054462, -1.1054462, 0.4123394, 3.7804907], [-11.7424810, -5.5984769, -1.0530628, -1.0530628, 0.4639595, 3.7634752], [-11.8475041, -5.3235151, -0.9967217, -0.9967217, 0.5190045, 3.7374198], [-11.9485698, -5.0398534, -0.9361524, -0.9361524, 0.5770772, 3.7013315], [-12.0451317, -4.7483265, -0.8728442, -0.8728442, 0.6374000, 3.6542898], [-12.1367486, -4.4497813, -0.8062747, -0.8062747, 0.6994645, 3.5962345], [-12.2229577, -4.1451630, -0.7363378, -0.7363378, 0.7620738, 3.5275349], [-12.3036825, -3.8342616, -0.6645410, -0.6645410, 0.8239388, 3.4492656], [-12.3785813, -3.5182948, -0.5910173, -0.5910173, 0.8831355, 3.3628938], [-12.4474598, -3.1982813, -0.5165164, -0.5165164, 0.9371063, 3.2703715], [-12.5101417, -2.8742981, -0.4414176, -0.4414176, 0.9824906, 3.1742440], [-12.5665300, -2.5476922, -0.3676221, -0.3676221, 1.0139670, 3.0765339], [-12.6164246, -2.2193389, -0.2953584, -0.2953584, 1.0250304, 2.9803149], [-12.6599277, -1.8901481, -0.2272997, -0.2272997, 1.0074675, 2.8878054], [-12.6967979, -1.5618667, -0.1646695, -0.1646695, 0.9513289, 2.8028954], [-12.7270124, -1.2357695, -0.1093675, -0.1093675, 0.8490177, 2.7279574], [-12.7505865, -0.9143992, -0.0634384, -0.0634384, 0.6962990, 2.6658768], [-12.7674865, -0.6001149, -0.0288569, -0.0288569, 0.4971249, 2.6191233], [-12.7776127, -0.2963852, -0.0073184, -0.0073184, 0.2595115, 2.5902191], [-12.7809896, -0.0154190, -0.0000000, 0.0000000, 0.0000014, 2.5804163], [-12.7764867, -0.3308468, -0.0230582, -0.0230573, 0.3140514, 2.5637405], [-12.7629551, -0.6543759, -0.0892045, -0.0892036, 0.6331626, 2.5157274], [-12.7404527, -0.9778316, -0.1908586, -0.1908577, 0.9482695, 2.4393288], [-12.7089646, -1.3008296, -0.3188353, -0.3188345, 1.2569646, 2.3404025], [-12.6685621, -1.6228939, -0.4650497, -0.4650489, 1.5566493, 2.2248025], [-12.6191620, -1.9439592, -0.6224488, -0.6224481, 1.8429584, 2.0982188], [-12.5608942, -2.2633852, -0.7862848, -0.7862841, 1.9648697, 2.1122420], [-12.4936865, -2.5812407, -0.9522076, -0.9522069, 1.8286174, 2.3570524], [-12.4176282, -2.8974353, -1.1179247, -1.1179240, 1.6923102, 2.5683301], [-12.3327623, -3.2118723, -1.2814488, -1.2814481, 1.5582483, 2.7357810], [-12.2391378, -3.5244655, -1.4412943, -1.4412937, 1.4282329, 2.8492726], [-12.1368617, -3.8345528, -1.5966825, -1.5966819, 1.3036773, 2.9021857], [-12.0258901, -4.1428768, -1.7463340, -1.7463334, 1.1859988, 2.8933847], [-11.9063490, -4.4494187, -1.8893549, -1.8893544, 1.0754987, 2.8299471], [-11.7783121, -4.7537285, -2.0254504, -2.0254498, 0.9732114, 2.7232477], [-11.6418640, -5.0556738, -2.1541732, -2.1541727, 0.8797340, 2.5858449], [-11.4970934, -5.3547915, -2.2748240, -2.2748234, 0.7960928, 2.4290387], [-11.3442460, -5.6516233, -2.3877983, -2.3877978, 0.7213844, 2.2611609], [-11.1833505, -5.9454923, -2.4926342, -2.4926337, 0.6568716, 2.0895370], [-11.0142769, -6.2362703, -2.5887865, -2.5887860, 0.6028755, 1.9186700], [-10.8373464, -6.5237299, -2.6762440, -2.6762435, 0.5595482, 1.7519927], [-10.6526997, -6.8076286, -2.7547773, -2.7547768, 0.5273008, 1.5919554], [-10.4604888, -7.0877124, -2.8244129, -2.8244124, 0.5060063, 1.4403041], [-10.2608763, -7.3633925, -2.8849534, -2.8849529, 0.4960346, 1.2984584], [-10.0540366, -7.6350603, -2.9363304, -2.9363299, 0.4975610, 1.1669205], [-9.8401560, -7.9021088, -2.9784594, -2.9784589, 0.5107419, 1.0465043], [-9.6194322, -8.1642656, -3.0112822, -3.0112817, 0.5357148, 0.9376281], [-9.3920670, -8.4211352, -3.0346650, -3.0346645, 0.5726993, 0.8406309], [-9.1582657, -8.6727049, -3.0483993, -3.0483988, 0.6216152, 0.7555448], [-8.9183095, -8.9183093, -3.0528528, -3.0528523, 0.6826180, 0.6826181]] test "every band at every point of the path agrees with the numpy reference (eV)" { assert_lt(max(abs(EVB[:, 0:shape(REFPATH)[1]] - REFPATH)), 1.0e-5) } # --- the physics: published LDA values, loose on purpose ----------------------- NV = NELEC / 2.0 # At Gamma the four levels around the gap are the Gamma25' TRIPLET (the valence-band top) # and the s-like Gamma2' SINGLET. In LDA the singlet comes out on the wrong side, so find # them by their degeneracy rather than by their order. E1 = EVB[GAMMA_P, NV - 3.0] E2 = EVB[GAMMA_P, NV - 2.0] E3 = EVB[GAMMA_P, NV - 1.0] E4 = EVB[GAMMA_P, NV] #: G2P the Gamma2' singlet, from the triplet. Measured: +0.89 eV, the direct gap. G2P = E4 #: TRIP the spread of the three-fold level: zero if the density has the crystal's symmetry TRIP = E3 - E1 if E2 - E1 > E4 - E3 { G2P = E1 TRIP = E4 - E2 } #: LC, XC the lowest empty state at L and at X LC = EVB[0.0, NV] XC = EVB[LAST, NV] show G2P show LC show XC test "the valence-band top at Gamma is a TRIPLET: the density has the crystal's full symmetry" { assert_near(TRIP, 0.0, 1.0e-4) } test "valence band width is about 12.7 eV" { assert_near(EVB[GAMMA_P, 0.0], -12.75, 0.3) } test "LDA's known failure: germanium has NO gap -- Gamma2' lands on the valence-band top (measured: 0.89 eV above it)" { assert_near(G2P, 0.0, 0.15) } test "the L valley is within 0.3 eV of the valence-band top (measured gap, at L: 0.74 eV) and below X, as in real germanium" { assert_gt(LC, 0.0) assert_lt(LC, 0.3) assert_lt(LC, XC) } figure("Germanium: LDA band structure", "wavevector along the path", "energy from the valence-band maximum (eV)") for b to NBANDS { plot(PATHX, EVB[:, b], "band") } xtick("L", PATHX[0]) xtick("Γ", PATHX[GAMMA_P]) xtick("X", PATHX[LAST])