Coverage for gpaw/test/lcao/test_largecellforce.py: 71%

1# This test calculates the force on the atoms in a hydrogen chain, comparing 

2# the results to finite-difference values. 

3# 

4# The reference result is obtained from an FD calculation, which can be rerun 

5# by setting the fd boolean below. 

6# 

7# The purpose is to test domain decomposition with large cells. The basis 

8# functions of one atom are defined to not overlap the rightmost domain 

9# for z-decompositions of two or more slices. This will change the set 

10# of atomic indices stored in BasisFunctions objects and other things 

11# 

12# This test also ensures that lcao forces are tested with non-pbc. 

13import numpy as np 

14from numpy import array 

15from ase import Atoms 

16from gpaw import GPAW 

17from gpaw.atom.basis import BasisMaker 

18 

19 

20def test_lcao_largecellforce(gpaw_new): 

21 hbasis = BasisMaker.from_symbol('H').generate(1, 0, energysplit=1.8, 

22 tailnorm=0.03**.5) 

23 basis = {'H': hbasis} 

24 

25 atom = Atoms('H') 

26 atom.center(vacuum=.8) 

27 system = atom.repeat((1, 1, 4)) 

28 

29 system.center(vacuum=2.5) 

30 

31 calc = GPAW(h=0.23, 

32 mode='lcao', 

33 basis=basis, 

34 convergence={'density': 1e-4, 'energy': 1e-7}) 

35 

36 system.calc = calc 

37 

38 F_ac = system.get_forces() 

39 

40 # Check that rightmost domain is in fact outside range of basis functions 

41 from gpaw.mpi import rank, size 

42 if rank == 0 and size > 1: 

43 if gpaw_new: 

44 basis = calc.dft.scf_loop.hamiltonian.basis 

45 else: 

46 basis = calc.wfs.basis_functions 

47 assert len(basis.atom_indices) < len(system) 

48 

49 fd = 0 

50 

51 # Values taken from FD calculation below 

52 # (Symmetry means only z-component may be nonzero) 

53 ref = array([[0, 0, 4.616841597363752], 

54 [0, 0, -2.7315136482540803], 

55 [0, 0, 2.7315116638237935], 

56 [0, 0, -4.616840606709416]]) 

57 

58 if fd: 

59 from gpaw.test import calculate_numerical_forces 

60 ref = calculate_numerical_forces(system, 0.002, icarts=[2]) 

61 print('Calced') 

62 print(F_ac) 

63 print('FD') 

64 print(ref) 

65 print(repr(ref)) 

66 

67 err = np.abs(F_ac - ref).max() 

68 print('Ref') 

69 print(ref) 

70 print('Calculated') 

71 print(F_ac) 

72 print('Max error', err) 

73 assert err < 6e-4