Merge pull request #196 from timcallow/fix_hartree_error

Fix an error in the Hartree potential

atoMEC/staticKS.py

@@ -930,7 +930,8 @@ def calc_v_ha(density, xgrid, grid_type):
             \exp(3x') -\int_x^{\log(r_s)} \mathrm{d}x' n(x') \exp(2x') \Big\}
         """
         # initialize v_ha
-        v_ha = np.zeros_like(xgrid)
+        v_ha_u = np.zeros_like(xgrid)
+        v_ha_l = np.zeros_like(xgrid)
 
         # construct the total (sum over spins) density
         rho = np.sum(density, axis=0)
@@ -948,14 +949,16 @@ def calc_v_ha(density, xgrid, grid_type):
 
             # now compute the hartree potential
             if grid_type == "log":
-                int_l = exp(-x0) * np.trapz(rho_l * np.exp(3.0 * x_l), x_l)
-                int_u = np.trapz(rho_u * np.exp(2 * x_u), x_u)
+                v_ha_l[i] = exp(-x0) * np.trapz(rho_l * np.exp(3.0 * x_l), x_l)
+                v_ha_u[i] = np.trapz(rho_u * np.exp(2 * x_u), x_u)
             else:
-                int_l = (1 / x0**2) * np.trapz(rho_l * 2 * x_l**5, x_l)
-                int_u = np.trapz(rho_u * 2 * x_u**3, x_u)
+                v_ha_l[i] = (1 / x0**2) * np.trapz(rho_l * 2 * x_l**5, x_l)
+                v_ha_u[i] = np.trapz(rho_u * 2 * x_u**3, x_u)
 
-            # total hartree potential is sum over integrals
-            v_ha[i] = 4.0 * pi * (int_l + int_u)
+        # total hartree potential is sum over integrals
+        # have to shift upper integral by one index to match
+        v_ha_u[1:] = v_ha_u[:-1]
+        v_ha = 4.0 * pi * (v_ha_l + v_ha_u)
 
         return v_ha
 

tests/boundary_conditions_test.py

@@ -12,9 +12,10 @@
 
 
 # expected values and tolerance
-dirichlet_expected = -14.080034
-neumann_expected = -15.161942
-bands_expected = -14.731577
+dirichlet_expected = -14.02657778662206
+neumann_expected = -15.112195227107968
+bands_expected = -14.68311960864459
+
 accuracy = 1e-3
 
 
@@ -80,6 +81,9 @@ def _run(bc):
 
 if __name__ == "__main__":
     config.numcores = -1
-    print("dirichlet = ", TestBcs._run("dirichlet"))
-    print("neumann = ", TestBcs._run("neumann"))
-    print("bands = ", TestBcs._run("bands"))
+    dirichlet = TestBcs._run("dirichlet")
+    neumann = TestBcs._run("neumann")
+    bands = TestBcs._run("bands")
+    print("dirichlet_expected =", dirichlet)
+    print("neumann_expected =", neumann)
+    print("bands_expected =", bands)

tests/conductivity_test.py

@@ -9,21 +9,22 @@
 
 
 # expected values and tolerance
-N_cc_expected_1_0 = 1.742515706139
-N_cc_expected_3_1 = 4.492209106738
-N_tt_expected_1_0 = 4.59790383269
-N_tt_expected_3_1 = 4.59790383269
+N_cc_expected_1_0 = 2.454442757747575
+N_cc_expected_3_1 = 6.727831527604185
+N_tt_expected_1_0 = 6.833387102730367
+N_tt_expected_3_1 = 6.833387102730367
 N_vv_expected_1_0 = 0.0
 N_vv_expected_3_1 = 0.0
-N_cv_expected_1_0 = 2.60323659985
-N_cv_expected_3_1 = 0.006759
-expected_integral_4 = 1.341640786499
-expected_integral_2 = 0.447213595499
-expected_sum_rule = 0.2940819621
+N_cv_expected_1_0 = 4.1257136059769985
+N_cv_expected_3_1 = 0.00749276366862972
+expected_integral_4 = 1.3416407864998738
+expected_integral_2 = 0.447213595499958
+expected_sum_rule = 0.29294247872666707
 expected_prop_sig_vv = 0.0
-expected_prop_sig_cv = 0.0013982278
-expected_prop_N_tot = 6.09241213503
-expected_prop_N_free = 4.5709301608
+expected_prop_sig_cv = 0.0014275605940770114
+expected_prop_N_tot = 6.08906292713134
+expected_prop_N_free = 4.5527425486281805
+
 accuracy = 0.001
 
 
@@ -238,18 +239,18 @@ def _run_sum_rule(input_SCF):
 if __name__ == "__main__":
     config.numcores = -1
     SCF_out = TestConductivity._run_SCF()
-    print(TestConductivity._run_cond_tot(SCF_out, "cc", (1, 0)))
-    print(TestConductivity._run_cond_tot(SCF_out, "cc", (3, 1)))
-    print(TestConductivity._run_cond_tot(SCF_out, "tt", (1, 0)))
-    print(TestConductivity._run_cond_tot(SCF_out, "tt", (3, 1)))
-    print(TestConductivity._run_cond_tot(SCF_out, "vv", (1, 0)))
-    print(TestConductivity._run_cond_tot(SCF_out, "vv", (3, 1)))
-    print(TestConductivity._run_cond_tot(SCF_out, "cv", (1, 0)))
-    print(TestConductivity._run_cond_tot(SCF_out, "cv", (3, 1)))
-    print(TestConductivity._run_int_calc(4))
-    print(TestConductivity._run_int_calc(2))
-    print(TestConductivity._run_sum_rule(SCF_out))
-    print(TestConductivity._run_prop(SCF_out, "sig_vv"))
-    print(TestConductivity._run_prop(SCF_out, "sig_cv"))
-    print(TestConductivity._run_prop(SCF_out, "N_tot"))
-    print(TestConductivity._run_prop(SCF_out, "N_free"))
+    print("N_cc_expected_1_0 =", TestConductivity._run_cond_tot(SCF_out, "cc", (1, 0)))
+    print("N_cc_expected_3_1 =", TestConductivity._run_cond_tot(SCF_out, "cc", (3, 1)))
+    print("N_tt_expected_1_0 =", TestConductivity._run_cond_tot(SCF_out, "tt", (1, 0)))
+    print("N_tt_expected_3_1 =", TestConductivity._run_cond_tot(SCF_out, "tt", (3, 1)))
+    print("N_vv_expected_1_0 =", TestConductivity._run_cond_tot(SCF_out, "vv", (1, 0)))
+    print("N_vv_expected_3_1 =", TestConductivity._run_cond_tot(SCF_out, "vv", (3, 1)))
+    print("N_cv_expected_1_0 =", TestConductivity._run_cond_tot(SCF_out, "cv", (1, 0)))
+    print("N_cv_expected_3_1 =", TestConductivity._run_cond_tot(SCF_out, "cv", (3, 1)))
+    print("expected_integral_4 =", TestConductivity._run_int_calc(4))
+    print("expected_integral_2 =", TestConductivity._run_int_calc(2))
+    print("expected_sum_rule =", TestConductivity._run_sum_rule(SCF_out))
+    print("expected_prop_sig_vv =", TestConductivity._run_prop(SCF_out, "sig_vv"))
+    print("expected_prop_sig_cv =", TestConductivity._run_prop(SCF_out, "sig_cv"))
+    print("expected_prop_N_tot =", TestConductivity._run_prop(SCF_out, "N_tot"))
+    print("expected_prop_N_free =", TestConductivity._run_prop(SCF_out, "N_free"))

tests/energy_alt_test.py

@@ -12,8 +12,8 @@
 
 
 # expected values and tolerance
-ideal_expected = -214.953809
-quantum_expected = -208.999350
+ideal_expected = -214.33339222837958
+quantum_expected = -207.98637746900067
 accuracy = 1e-3
 
 
@@ -86,5 +86,7 @@ def _run(unbound):
 
 if __name__ == "__main__":
     config.numcores = -1
-    print("ideal energy = ", TestEnergyAlt._run("ideal"))
-    print("quantum energy = ", TestEnergyAlt._run("quantum"))
+    ideal = TestEnergyAlt._run("ideal")
+    quantum = TestEnergyAlt._run("quantum")
+    print("ideal_expected =", ideal)
+    print("quantum_expected =", quantum)

tests/functionals_test.py

@@ -12,11 +12,11 @@
 
 
 # expected values and tolerance
-lda_expected = -159.210841
-gdsmfb_expected = -159.154166
-no_xc_expected = -145.019411
-no_hxc_expected = -249.507728
-gga_expected = -159.92294751868772
+lda_expected = -158.66453567079228
+gdsmfb_expected = -158.6079090856573
+no_xc_expected = -144.49795725080637
+no_hxc_expected = -249.50772961499766
+gga_expected = -159.37528906276566
 accuracy = 1e-3
 
 
@@ -100,8 +100,13 @@ def _run(func):
 
 if __name__ == "__main__":
     config.numcores = -1
-    print("lda=", TestFuncs()._run("lda"))
-    print("gdsmfb=", TestFuncs()._run("gdsmfb"))
-    print("no_xc=", TestFuncs()._run("no_xc"))
-    print("no_hxc=", TestFuncs()._run("no_hxc"))
-    print("pbe=", TestFuncs()._run("gga"))
+    lda = TestFuncs()._run("lda")
+    gdsmfb = TestFuncs()._run("gdsmfb")
+    no_xc = TestFuncs()._run("no_xc")
+    no_hxc = TestFuncs()._run("no_hxc")
+    pbe = TestFuncs()._run("gga")
+    print("lda_expected =", lda)
+    print("gdsmfb_expected =", gdsmfb)
+    print("no_xc_expected =", no_xc)
+    print("no_hxc_expected =", no_hxc)
+    print("gga_expected =", pbe)

tests/gramschmidt_test.py

@@ -13,7 +13,7 @@
 
 # expected values and tolerance
 self_overlap_expected = 1.00000
-overlap_expected = 0.041998602727
+overlap_expected = 0.041051392274703176
 accuracy = 0.0001
 
 
@@ -97,5 +97,5 @@ def _run_overlap(input_SCF, case):
 
 if __name__ == "__main__":
     SCF_out = TestGS._run_SCF()
-    print(TestGS._run_overlap(SCF_out, "self"))
-    print(TestGS._run_overlap(SCF_out, "other"))
+    print("self_overlap_expected =", TestGS._run_overlap(SCF_out, "self"))
+    print("overlap_expected =", TestGS._run_overlap(SCF_out, "other"))

tests/localization_test.py

@@ -13,11 +13,11 @@
 
 
 # expected values and tolerance
-epdc_orbs_expected = 9.2269
-epdc_dens_expected = 3.5304
-orbitals_expected = 1.0889
-density_expected = 2.1496
-IPR_expected = 78.2107
+epdc_orbs_expected = 9.22831052983569
+epdc_dens_expected = 3.527486593490338
+density_expected = 2.1469328620811075
+orbitals_expected = 1.100028932668539
+IPR_expected = 78.19952984379802
 accuracy = 0.01
 
 
@@ -221,8 +221,19 @@ def _run_IPR(input_SCF):
     config.numcores = -1
     SCF_spin_out = TestLocalization._run_SCF(True)
     SCF_no_spin_out = TestLocalization._run_SCF(False)
-    print(TestLocalization._run_epdc(SCF_spin_out, "orbitals", True))
-    print(TestLocalization._run_epdc(SCF_spin_out, "density", True))
-    print(TestLocalization._run_ELF(SCF_no_spin_out, "density", False))
-    print(TestLocalization._run_ELF(SCF_spin_out, "orbitals", True))
-    print(TestLocalization._run_IPR(SCF_no_spin_out))
+    print(
+        "epdc_orbs_expected =",
+        TestLocalization._run_epdc(SCF_spin_out, "orbitals", True),
+    )
+    print(
+        "epdc_dens_expected =",
+        TestLocalization._run_epdc(SCF_spin_out, "density", True),
+    )
+    print(
+        "density_expected =",
+        TestLocalization._run_ELF(SCF_no_spin_out, "density", False),
+    )
+    print(
+        "orbitals_expected =", TestLocalization._run_ELF(SCF_spin_out, "orbitals", True)
+    )
+    print("IPR_expected =", TestLocalization._run_IPR(SCF_no_spin_out))

tests/pressure_test_log.py

@@ -15,14 +15,15 @@
 
 
 # expected values and tolerance
-finite_diff_expected_A = 170.95
-finite_diff_expected_B = 239.70
-stress_tensor_expected_rr = 146.06
-stress_tensor_expected_tr = 109.76
-virial_expected_corr = 140.75
-virial_expected_nocorr = 175.89
-ideal_expected = 100.06
-ion_expected = 165.19
+finite_diff_expected_A = 172.56183194182995
+finite_diff_expected_B = 243.88427030813003
+stress_tensor_expected_rr = 149.58033991478263
+stress_tensor_expected_tr = 111.04150046097219
+virial_expected_corr = 143.06375646995878
+virial_expected_nocorr = 178.2430500359269
+ideal_expected = 103.0141806222132
+ion_expected = 165.19118614722603
+
 accuracy = 0.1
 
 
@@ -245,17 +246,19 @@ def _run_ion(Atom):
 if __name__ == "__main__":
     config.numcores = -1
     SCF_out = TestPressure._run_SCF()
-    print("Finite diff pressure A: ", TestPressure._run_finite_diff(SCF_out, "A"))
-    print("Finite diff pressure B: ", TestPressure._run_finite_diff(SCF_out, "B"))
+    fd_a = TestPressure._run_finite_diff(SCF_out, "A")
+    fd_b = TestPressure._run_finite_diff(SCF_out, "B")
+    print("finite_diff_expected_A =", fd_a)
+    print("finite_diff_expected_B =", fd_b)
     print(
-        "Stress tensor pressure rr: ",
+        "stress_tensor_expected_rr =",
         TestPressure._run_stress_tensor(SCF_out, True),
     )
     print(
-        "Stress tensor pressure tr: ",
+        "stress_tensor_expected_tr =",
         TestPressure._run_stress_tensor(SCF_out, False),
     )
-    print("Virial pressure corr: ", TestPressure._run_virial(SCF_out, True))
-    print("Virial pressure no corr: ", TestPressure._run_virial(SCF_out, False))
-    print("Ideal electron: ", TestPressure._run_ideal(SCF_out))
-    print("Ion pressure: ", TestPressure._run_ion(SCF_out["Atom"]))
+    print("virial_expected_corr =", TestPressure._run_virial(SCF_out, True))
+    print("virial_expected_nocorr =", TestPressure._run_virial(SCF_out, False))
+    print("ideal_expected =", TestPressure._run_ideal(SCF_out))
+    print("ion_expected =", TestPressure._run_ion(SCF_out["Atom"]))

tests/pressure_test_sqrt.py

@@ -15,13 +15,14 @@
 
 
 # expected values and tolerance
-finite_diff_expected_A = 10.7309
-finite_diff_expected_B = 24.2916
-stress_tensor_expected_rr = 9.8064
-stress_tensor_expected_tr = 7.0468
-virial_expected_corr = 10.9138
-virial_expected_nocorr = 11.7581
-ideal_expected = 8.4891
+finite_diff_expected_A = 10.6938851494863
+finite_diff_expected_B = 24.390856939007094
+stress_tensor_expected_rr = 9.831719531837031
+stress_tensor_expected_tr = 7.059453507314628
+virial_expected_corr = 10.89986054089238
+virial_expected_nocorr = 11.744866168684108
+ideal_expected = 8.514959799814081
+
 accuracy = 0.1
 
 
@@ -220,16 +221,18 @@ def _run_ideal(scf_input):
 if __name__ == "__main__":
     config.numcores = -1
     SCF_out = TestPressure._run_SCF()
-    print("Finite diff pressure A: ", TestPressure._run_finite_diff(SCF_out, "A"))
-    print("Finite diff pressure B: ", TestPressure._run_finite_diff(SCF_out, "B"))
+    fd_a = TestPressure._run_finite_diff(SCF_out, "A")
+    fd_b = TestPressure._run_finite_diff(SCF_out, "B")
+    print("finite_diff_expected_A =", fd_a)
+    print("finite_diff_expected_B =", fd_b)
     print(
-        "Stress tensor pressure rr: ",
+        "stress_tensor_expected_rr =",
         TestPressure._run_stress_tensor(SCF_out, True),
     )
     print(
-        "Stress tensor pressure tr: ",
+        "stress_tensor_expected_tr =",
         TestPressure._run_stress_tensor(SCF_out, False),
     )
-    print("Virial pressure corr: ", TestPressure._run_virial(SCF_out, True))
-    print("Virial pressure no corr: ", TestPressure._run_virial(SCF_out, False))
-    print("Ideal electron: ", TestPressure._run_ideal(SCF_out))
+    print("virial_expected_corr =", TestPressure._run_virial(SCF_out, True))
+    print("virial_expected_nocorr =", TestPressure._run_virial(SCF_out, False))
+    print("ideal_expected =", TestPressure._run_ideal(SCF_out))

tests/serial_test.py

@@ -7,8 +7,8 @@
 
 
 # expected values and tolerance
-dense_expected = -0.562437
-coarse_expected = -0.567817
+dense_expected = -0.5620194349606303
+coarse_expected = -0.5625664849484403
 accuracy = 1e-3
 
 
@@ -68,4 +68,8 @@ def _run(ngrid):
 
 
 if __name__ == "__main__":
-    print(TestSerial._run())
+    config.numcores = 0
+    dense = TestSerial._run(5001)
+    coarse = TestSerial._run(400)
+    print("dense_expected =", dense)
+    print("coarse_expected =", coarse)

tests/spin_test.py

@@ -7,8 +7,9 @@
 
 
 # expected values and tolerance
-singlet_expected = -38.0024097674
-triplet_expected = -37.987537447
+singlet_expected = -37.89807041319449
+triplet_expected = -37.88174335076866
+
 accuracy = 1e-4
 
 
@@ -76,5 +77,7 @@ def _run(spinmag):
 
 if __name__ == "__main__":
     config.numcores = -1
-    print("singlet energy = ", TestSpin._run(0))
-    print("triplet energy = ", TestSpin._run(2))
+    singlet = TestSpin._run(0)
+    triplet = TestSpin._run(2)
+    print("singlet_expected =", singlet)
+    print("triplet_expected =", triplet)