Coverage for gpaw/test/noncollinear/test_complex_densmat.py: 100%

1import numpy as np 

2import pytest 

3 

4from gpaw.new.pot_calc import calculate_non_local_potential1 

5from gpaw.new.xc import create_functional 

6from gpaw.setup import create_setup 

7from gpaw.core import UGDesc 

8from gpaw.xc import XC 

9 

10 

11@pytest.mark.soc 

12def test_energy_from_complex_densmat(): 

13 

14 # Set up objects and matrices 

15 

16 setup = create_setup('Ga') 

17 grid = UGDesc(cell=[1, 1, 1], size=[9, 9, 9]) 

18 xc = create_functional(XC('LDA', collinear=False), grid) 

19 soc = True 

20 err = 1.0e-6 

21 D_sii = np.zeros((4, setup.ni, setup.ni), complex) 

22 D_sii[0, 0, 0] = 2 

23 

24 # Calculate non-local potential energy for some random state 

25 

26 P_sm = np.array([[0.2 + 0.3j, 0.1 + 0.2j, 0.3 + 0.4j], 

27 [0.4 - 0.5j, 0.2 + 0.3j, 0.6 - 0.7j]]) 

28 

29 D_ssmm = np.einsum('si, zj -> szij', P_sm.conj(), P_sm) 

30 D_sii[:, 1:4, 1:4] = [D_ssmm[0, 0] + D_ssmm[1, 1], 

31 D_ssmm[0, 1] + D_ssmm[1, 0], 

32 -1j * (D_ssmm[0, 1] - D_ssmm[1, 0]), 

33 D_ssmm[0, 0] - D_ssmm[1, 1]] 

34 

35 def calc_energies(D_sii): 

36 _, energies = calculate_non_local_potential1( 

37 setup, xc, D_sii, np.zeros(1), soc, [], 0) 

38 return energies 

39 

40 energies1 = calc_energies(D_sii) 

41 

42 assert energies1['kinetic_correction'] == pytest.approx( 

43 0.04340694003, abs=err) 

44 assert energies1['coulomb'] == pytest.approx(-5.5575386716, abs=err) 

45 assert energies1['zero'] == pytest.approx(-2.432694074696, abs=err) 

46 assert energies1['xc'] == pytest.approx(1.5938337327, abs=err) 

47 

48 # Rotate the state 90 degrees around the z-axis (x -> y, y -> -x, z -> z). 

49 # Assert that this does not change the energies. 

50 

51 # First rotate the spins 

52 P_sm = [P_sm[0, :] * (1 - 1j), P_sm[1, :] * (1 + 1j)] / np.sqrt(2) 

53 # Then rotate the density 

54 P_sm = np.matmul([[0, 0, 1], [0, 1, 0], [-1, 0, 0]], P_sm.T).T 

55 

56 D_ssmm = np.einsum('si, zj -> szij', P_sm.conj(), P_sm) 

57 D_sii[:, 1:4, 1:4] = [D_ssmm[0, 0] + D_ssmm[1, 1], 

58 D_ssmm[0, 1] + D_ssmm[1, 0], 

59 -1j * (D_ssmm[0, 1] - D_ssmm[1, 0]), 

60 D_ssmm[0, 0] - D_ssmm[1, 1]] 

61 

62 energies2 = calc_energies(D_sii) 

63 

64 assert energies2['kinetic_correction'] == pytest.approx( 

65 energies1['kinetic_correction'], abs=err) 

66 assert energies2['coulomb'] == pytest.approx(energies1['coulomb'], abs=err) 

67 assert energies2['zero'] == pytest.approx(energies1['zero'], abs=err) 

68 assert energies2['xc'] == pytest.approx(energies1['xc'], abs=err) 

69 

70 # Assert that only the kinetic energy changes when the density 

71 # matrix is forced to be real 

72 

73 energies3 = calc_energies(D_sii.real) 

74 

75 assert energies3['kinetic_correction'] == pytest.approx( 

76 0.0446930609623, abs=err) 

77 assert energies3['coulomb'] == pytest.approx(energies1['coulomb'], abs=err) 

78 assert energies3['zero'] == pytest.approx(energies1['zero'], abs=err) 

79 assert energies3['xc'] == pytest.approx(energies1['xc'], abs=err)