Coverage for gpaw/test/generic/test_guc_force.py: 75%

1# This tests calculates the force on the atoms of a 

2# slightly distorted Silicon primitive cell. 

3# 

4# If the test fails, set the fd boolean below to enable a (costly) finite 

5# difference check. 

6 

7import numpy as np 

8import pytest 

9from ase import Atoms 

10 

11from gpaw import GPAW 

12from gpaw.atom.basis import BasisMaker 

13 

14 

15def test_generic_guc_force(): 

16 sibasis = BasisMaker.from_symbol('Si').generate( 

17 2, 1, energysplit=0.3, tailnorm=0.03**.5) 

18 basis = {'Si': sibasis} 

19 

20 a = 5.475 

21 system = Atoms(symbols='Si2', pbc=True, 

22 cell=0.5 * a * np.array([(1, 1, 0), 

23 (1, 0, 1), 

24 (0, 1, 1)]), 

25 scaled_positions=[(0.0, 0.0, 0.0), 

26 (0.23, 0.23, 0.23)]) 

27 

28 calc = GPAW(h=0.2, 

29 mode='lcao', 

30 basis=basis, 

31 kpts=(2, 2, 2), 

32 convergence={'density': 1e-5, 'energy': 1e-6} 

33 ) 

34 system.calc = calc 

35 

36 F_ac = system.get_forces() 

37 

38 # Previous FD result, generated by disabled code below 

39 F_ac_ref = np.array( 

40 [[-1.3967114867039498, -1.3967115816022613, -1.396711581510779], 

41 [1.397400325299003, 1.3974003410455182, 1.3974003410801572]]) 

42 

43 err_ac = np.abs(F_ac - F_ac_ref) 

44 err = err_ac.max() 

45 

46 print('Force') 

47 print(F_ac) 

48 print() 

49 print('Reference result') 

50 print(F_ac_ref) 

51 print() 

52 print('Error') 

53 print(err_ac) 

54 print() 

55 print('Max error') 

56 print(err) 

57 

58 # ASE uses dx = [+|-] 0.001 by default, 

59 # error should be around 2e-3. In fact 4e-3 would probably be acceptable 

60 

61 assert err == pytest.approx(0, abs=6e-3) 

62 

63 # Set boolean to run new FD check 

64 fd = False 

65 

66 if fd: 

67 from gpaw.test import calculate_numerical_forces 

68 system.calc = calc.new(symmetry='off') 

69 F_ac_fd = calculate_numerical_forces(system, 0.001) 

70 print('Self-consistent forces') 

71 print(F_ac) 

72 print('FD') 

73 print(F_ac_fd) 

74 print(repr(F_ac_fd)) 

75 print(F_ac - F_ac_fd, np.abs(F_ac - F_ac_fd).max()) 

76 

77 assert np.abs(F_ac - F_ac_fd).max() < 4e-3