Coverage for gpaw/test/lcao/test_force.py: 81%

1# This tests calculates the force on the atoms of a small molecule. 

2# 

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

4# difference check. 

5 

6import numpy as np 

7import pytest 

8from ase.build import molecule 

9 

10from gpaw import GPAW 

11from gpaw.atom.basis import BasisMaker 

12 

13 

14@pytest.mark.old_gpaw_only # basis set cutoff? 

15def test_lcao_force(): 

16 obasis = BasisMaker.from_symbol('O').generate(2, 1, energysplit=0.3, 

17 tailnorm=0.03**.5) 

18 hbasis = BasisMaker.from_symbol('H').generate(2, 1, energysplit=0.3, 

19 tailnorm=0.03**.5) 

20 basis = {'O': obasis, 'H': hbasis} 

21 

22 system = molecule('H2O') 

23 system.center(vacuum=1.5) 

24 system.rattle(stdev=.2, seed=42) 

25 system.set_pbc(1) 

26 

27 calc = GPAW(h=0.2, 

28 mode='lcao', 

29 basis=basis, 

30 kpts=[(0., 0., 0.), (.3, .1, .4)], 

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

32 

33 system.calc = calc 

34 

35 # Previous FD result, generated by disabled code below 

36 F_ac_ref = np.array([[1.05022478, 1.63103681, -5.00612007], 

37 [-0.69739179, -0.89624274, 3.03203147], 

38 [-0.34438181, -0.7042696, 1.77209023]]) 

39 

40 for use_rho in [0, 1]: 

41 if use_rho: 

42 for kpt in calc.wfs.kpt_u: 

43 kpt.rho_MM = calc.wfs.calculate_density_matrix(kpt.f_n, 

44 kpt.C_nM) 

45 

46 F_ac = system.get_forces() 

47 system.calc.results.pop('forces') 

48 

49 err_ac = np.abs(F_ac - F_ac_ref) 

50 err = err_ac.max() 

51 

52 print('Force') 

53 print(F_ac) 

54 print() 

55 print('Reference result') 

56 print(F_ac_ref) 

57 print() 

58 print('Error') 

59 print(err_ac) 

60 print() 

61 print('Max error') 

62 print(err) 

63 

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

65 # error should be around 2e-3. 

66 # In fact 4e-3 would probably be acceptable 

67 assert err < 3e-3 

68 

69 # Set boolean to run new FD check 

70 fd = False 

71 

72 if fd: 

73 from gpaw.test import calculate_numerical_forces 

74 F_ac_fd = calculate_numerical_forces(system, 0.001) 

75 print('Self-consistent forces') 

76 print(F_ac) 

77 print('FD') 

78 print(F_ac_fd) 

79 print(repr(F_ac_fd)) 

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

81 

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