Coverage for gpaw/test/sic/test_pwsic.py: 90%

1import io 

2import pytest 

3 

4from gpaw import GPAW, restart 

5import numpy as np 

6from ase.dft.bandgap import bandgap 

7from ase.units import Ha 

8 

9import numpy.testing as npt 

10from gpaw.io.logger import GPAWLogger 

11from gpaw.wavefunctions.base import eigenvalue_string 

12from gpaw.test.sic._utils import (mk_arr_from_str, 

13 extract_lagrange_section, 

14 MockWorld) 

15from gpaw.mpi import rank 

16 

17 

18@pytest.mark.old_gpaw_only 

19@pytest.mark.sic 

20def test_pwsic(in_tmp_dir, gpw_files): 

21 """ 

22 test Perdew-Zunger Self-Interaction 

23 Correction in PW mode using DirectMin 

24 :param in_tmp_dir: 

25 :return: 

26 """ 

27 calc = GPAW(gpw_files["h2o_pwsic"]) 

28 H2O = calc.atoms 

29 H2O.calc = calc 

30 e = H2O.get_potential_energy() 

31 f = H2O.get_forces() 

32 efermi = calc.wfs.fermi_levels[0] * Ha 

33 gap = bandgap(calc, efermi=efermi)[0] 

34 

35 # Numeric forces, generated by disabled code below 

36 f2 = np.array( 

37 [ 

38 [0.22161312, -0.98564396, -0.00204214], 

39 [-0.34986867, 0.17494903, 0.00029861], 

40 [0.01085528, 0.56112341, 0.00129632], 

41 ] 

42 ) 

43 assert f2 == pytest.approx(f, abs=3e-2) 

44 assert e == pytest.approx(-9.98919, abs=1e-3) 

45 assert gap == pytest.approx(9.555, abs=1e-2) 

46 

47 numeric = False 

48 if numeric: 

49 from gpaw.test import calculate_numerical_forces 

50 

51 f_num = calculate_numerical_forces(H2O, 0.001) 

52 print("Numerical forces") 

53 print(f_num) 

54 print(f - f_num, np.abs(f - f_num).max()) 

55 

56 # 

57 calc.write("h2o.gpw", mode="all") 

58 

59 H2O, calc = restart("h2o.gpw", txt="-") 

60 H2O.positions += 1.0e-6 

61 f3 = H2O.get_forces() 

62 niter = calc.get_number_of_iterations() 

63 assert niter == pytest.approx(4, abs=3) 

64 assert f2 == pytest.approx(f3, abs=3e-2) 

65 

66 if rank == 0: 

67 logger = GPAWLogger(MockWorld(rank=0)) 

68 string_io = io.StringIO() 

69 logger.fd = string_io 

70 calc.wfs.summary_func(logger) 

71 lstr = extract_lagrange_section(string_io.getvalue()) 

72 

73 expect_lagrange_str = """\ 

74 Band L_ii Occupancy Band L_ii Occupancy 

75 0 -21.28876 1.00000 0 -21.29052 1.00000 

76 1 -21.02682 1.00000 1 -21.03380 1.00000 

77 2 -13.96149 1.00000 2 -13.94529 1.00000 

78 3 -13.91558 1.00000 3 -13.92438 1.00000 

79 4 -0.94710 0.00000 4 -0.94712 0.00000 

80 5 0.76805 0.00000 5 0.76806 0.00000 

81 """ 

82 expect_eigen_str = """\ 

83 Band Eigenvalues Occupancy Eigenvalues Occupancy 

84 0 -30.18943 1.00000 -30.19005 1.00000 

85 1 -16.73073 1.00000 -16.73110 1.00000 

86 2 -12.77049 1.00000 -12.77136 1.00000 

87 3 -10.50200 1.00000 -10.50147 1.00000 

88 4 -0.94720 0.00000 -0.94720 0.00000 

89 5 0.76815 0.00000 0.76815 0.00000 

90 """ 

91 

92 npt.assert_allclose( 

93 mk_arr_from_str(expect_lagrange_str, 6), 

94 mk_arr_from_str(lstr, 6), 

95 atol=0.3, 

96 ) 

97 

98 npt.assert_allclose( 

99 mk_arr_from_str(expect_eigen_str, 5), 

100 mk_arr_from_str(eigenvalue_string(calc.wfs), 5, skip_rows=1), 

101 atol=0.3, 

102 )