Coverage for gpaw/test/directopt/test_mom_directopt_pw.py: 84%

1import pytest 

2 

3from gpaw import GPAW, restart 

4from gpaw.mom import prepare_mom_calculation 

5from gpaw.directmin.tools import excite 

6import numpy as np 

7 

8 

9@pytest.mark.do 

10def test_mom_directopt_pw(in_tmp_dir, gpaw_new, gpw_files): 

11 calc = GPAW(gpw_files['h2o_mom_do_pw']) 

12 atoms = calc.atoms 

13 atoms.calc = calc 

14 calc.write('h2o.gpw', mode='all') 

15 

16 for canonical in [True, False]: 

17 # Mixed-spin excited state calculation 

18 atoms, calc = restart('h2o.gpw', txt='-') 

19 mom_after_canonical = False # Test MOM after canonical only 

20 if mom_after_canonical: 

21 momevery = np.inf 

22 else: 

23 momevery = 3 

24 calc.set(eigensolver=dict(name='etdm-fdpw', 

25 excited_state=True, 

26 momevery=momevery, 

27 restart_canonical=canonical, 

28 printinnerloop=True)) 

29 f_sn = excite(calc, 0, 0, (0, 0)) 

30 prepare_mom_calculation(calc, atoms, f_sn) 

31 

32 def rotate_homo_lumo(ctx=None, calc=calc): 

33 a = 70 * np.pi / 180.0 

34 if gpaw_new: 

35 iters = ctx.niter 

36 else: 

37 iters = calc.get_number_of_iterations() 

38 if iters == 3: 

39 psit_nG_old = calc.wfs.kpt_u[0].psit_nG.copy() 

40 calc.wfs.kpt_u[0].psit_nG[3] = \ 

41 np.cos(a) * psit_nG_old[3] + np.sin(a) * psit_nG_old[4] 

42 calc.wfs.kpt_u[0].psit_nG[4] = \ 

43 np.cos(a) * psit_nG_old[4] - np.sin(a) * psit_nG_old[3] 

44 for kpt in calc.wfs.kpt_u: 

45 calc.wfs.pt.integrate(kpt.psit_nG, kpt.P_ani, kpt.q) 

46 

47 if gpaw_new: 

48 calc.hooks['scf_step'] = rotate_homo_lumo 

49 else: 

50 calc.attach(rotate_homo_lumo, 1) 

51 e = atoms.get_potential_energy() 

52 assert e == pytest.approx(0.027152, abs=1.0e-3) 

53 

54 f = atoms.get_forces() 

55 

56 # Numeric forces, generated by disabled code below 

57 f2 = np.array([[-4.07053122, -5.46395110, -2.54491307e-04], 

58 [5.57197214, -1.00373986e-01, 1.87598849e-04], 

59 [-1.52873845, 5.38474614, 2.07965790e-04]]) 

60 assert f2 == pytest.approx(f, abs=0.1) 

61 

62 numeric = False 

63 if numeric: 

64 calc.observers = [] 

65 from gpaw.test import calculate_numerical_forces 

66 f_num = calculate_numerical_forces(atoms, 0.001) 

67 print('Numerical forces') 

68 print(f_num) 

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

70 

71 calc.write('h2o.gpw', mode='all') 

72 

73 # Test restart and fixed occupations 

74 atoms, calc = restart('h2o.gpw', txt='-') 

75 if gpaw_new: 

76 return 

77 atoms.calc.results.pop('energy') 

78 atoms.calc.scf.converged = False 

79 e2 = atoms.get_potential_energy() 

80 niter = calc.get_number_of_iterations() 

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

82 assert e == pytest.approx(e2, abs=1.0e-3) 

83 

84 prepare_mom_calculation(calc, atoms, f_sn, use_fixed_occupations='True') 

85 e2 = atoms.get_potential_energy() 

86 for spin in range(calc.get_number_of_spins()): 

87 f_n = calc.get_occupation_numbers(spin=spin) 

88 assert (np.allclose(f_sn[spin], f_n)) 

89 assert (np.allclose(f_sn[spin], calc.wfs.occupations.numbers[spin])) 

90 niter = calc.get_number_of_iterations() 

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

92 assert e == pytest.approx(e2, abs=1.0e-3) 

93 

94 

95if __name__ == '__main__': 

96 test_mom_directopt_pw(1, 1)