Coverage for gpaw/test/mom/test_mom_lcao_forces.py: 100%

1import numpy as np 

2import pytest 

3from ase import Atoms 

4 

5from gpaw import GPAW, restart 

6from gpaw.mom import prepare_mom_calculation 

7 

8 

9@pytest.mark.mom 

10def test_mom_lcao_forces(in_tmp_dir): 

11 force_ref = 11.52 

12 

13 f_sn = [[1., 1., 1., 1., 0., 1., 0.], 

14 [1., 1., 1., 1., 1., 0., 0.]] 

15 L = 4.0 

16 d = 1.13 

17 delta = 0.01 

18 

19 atoms = Atoms('CO', 

20 [[0.5 * L, 0.5 * L, 0.5 * L - 0.5 * d], 

21 [0.5 * L, 0.5 * L, 0.5 * L + 0.5 * d]]) 

22 atoms.set_cell([L, L, L]) 

23 atoms.rotate(1, 'x', center=[0.5 * L, 0.5 * L, 0.5 * L]) 

24 

25 calc = GPAW(mode='lcao', 

26 basis='dzp', 

27 nbands=7, 

28 h=0.24, 

29 xc='PBE', 

30 spinpol=True, 

31 symmetry='off', 

32 convergence={'energy': 100, 

33 'density': 1e-4}) 

34 

35 atoms.calc = calc 

36 # Ground-state calculation 

37 atoms.get_potential_energy() 

38 calc.write('co_lcao_gs.gpw', 'all') 

39 

40 for mom in [False, True]: 

41 atoms, calc = restart('co_lcao_gs.gpw', txt='-') 

42 

43 occ = prepare_mom_calculation(calc, atoms, f_sn, use_projections=mom) 

44 F = atoms.get_forces() 

45 

46 # Test overlaps 

47 occ.initialize_reference_orbitals() 

48 for kpt in calc.wfs.kpt_u: 

49 f_n = calc.get_occupation_numbers(spin=kpt.s) 

50 P = occ.calculate_weights(kpt, 1.0) 

51 assert (np.allclose(P, f_n)) 

52 

53 E = [] 

54 p = atoms.positions.copy() 

55 for i in [-1, 1]: 

56 pnew = p.copy() 

57 pnew[0, 2] -= delta / 2. * i 

58 pnew[1, 2] += delta / 2. * i 

59 atoms.set_positions(pnew) 

60 

61 E.append(atoms.get_potential_energy()) 

62 

63 f = np.sqrt(((F[1, :] - F[0, :])**2).sum()) * 0.5 

64 fnum = (E[0] - E[1]) / (2. * delta) # central difference 

65 

66 print(fnum, f) 

67 assert fnum == pytest.approx(force_ref, abs=0.016) 

68 assert f == pytest.approx(fnum, abs=0.1)