Coverage for gpaw/test/pw/test_si_stress_mgga.py: 96%

1import numpy as np 

2import pytest 

3from ase.io.ulm import ulmopen 

4from ase.parallel import parprint 

5 

6from gpaw import GPAW 

7from gpaw.test import calculate_numerical_forces 

8 

9 

10@pytest.mark.mgga 

11def test_pw_si_stress_mgga(gpw_files, gpaw_new): 

12 if gpaw_new and ulmopen(gpw_files['si_pw_distorted']).version < 4: 

13 pytest.skip('Unsupported new-GPAW + old gpw-file combo') 

14 

15 si = GPAW(gpw_files['si_pw_distorted']).get_atoms() 

16 

17 # Trigger nasty bug (fixed in !486): 

18 if not gpaw_new: 

19 si.calc.wfs.pt.blocksize = si.calc.wfs.pd.maxmyng - 1 

20 

21 s_analytical = si.get_stress() 

22 # Calculated numerical stress once, store here to speed up test 

23 s_numerical = np.array([-0.01140242, -0.04084746, -0.0401058, 

24 -0.02119629, 0.13584242, 0.00911572]) 

25 if 0: 

26 s_numerical = si.calc.calculate_numerical_stress(si, 1e-5) 

27 

28 s_err = s_numerical - s_analytical 

29 

30 parprint('Analytical stress:\n', s_analytical) 

31 parprint('Numerical stress:\n', s_numerical) 

32 parprint('Error in stress:\n', s_err) 

33 assert np.all(abs(s_err) < 1e-5) 

34 

35 # Check y-component of second atom: 

36 f = si.get_forces()[1, 1] 

37 fref = -2.066952082010687 

38 if 0: 

39 fref = calculate_numerical_forces(si, 0.001, [1], [1])[0, 0] 

40 print(f, fref, f - fref) 

41 assert abs(f - fref) < 0.0005