Coverage for gpaw/test/generic/test_relax.py: 100%

1from time import time 

2 

3import pytest 

4from ase import Atoms 

5from ase.io import read 

6from ase.optimize import BFGS 

7 

8from gpaw import GPAW 

9from gpaw.mpi import world 

10 

11 

12def timer(func, *args, **kwargs): 

13 t0 = time() 

14 ret = func(*args, **kwargs) 

15 return ret, time() - t0 

16 

17 

18@pytest.mark.old_gpaw_only 

19def test_generic_relax(in_tmp_dir): 

20 a = 4.0 # Size of unit cell (Angstrom) 

21 c = a / 2 

22 d = 0.74 # Experimental bond length 

23 molecule = Atoms('H2', 

24 [(c - d / 2, c, c), 

25 (c + d / 2, c, c)], 

26 cell=(a, a, a), 

27 pbc=False) 

28 calc = GPAW(mode='fd', 

29 h=0.2, 

30 nbands=1, 

31 xc={'name': 'PBE', 'stencil': 1}, 

32 txt=None, 

33 poissonsolver={'name': 'fd'}) 

34 molecule.calc = calc 

35 e1 = molecule.get_potential_energy() 

36 calc.write('H2.gpw') 

37 calc.write('H2a.gpw', mode='all') 

38 molecule.get_forces() 

39 calc.write('H2f.gpw') 

40 calc.write('H2fa.gpw', mode='all') 

41 

42 molecule = GPAW('H2.gpw', txt=None).get_atoms() 

43 f1, t1 = timer(molecule.get_forces) 

44 molecule = GPAW('H2a.gpw', txt=None).get_atoms() 

45 f2, t2 = timer(molecule.get_forces) 

46 molecule = GPAW('H2f.gpw', txt=None).get_atoms() 

47 f3, t3 = timer(molecule.get_forces) 

48 molecule = GPAW('H2fa.gpw', txt=None).get_atoms() 

49 f4, t4 = timer(molecule.get_forces) 

50 print('timing:', t1, t2, t3, t4) 

51 assert t2 < 0.6 * t1 

52 assert t3 < 0.5 

53 assert t4 < 0.5 

54 print(f1) 

55 print(f2) 

56 print(f3) 

57 print(f4) 

58 assert sum((f1 - f4).ravel()**2) < 1e-6 

59 assert sum((f2 - f4).ravel()**2) < 1e-6 

60 assert sum((f3 - f4).ravel()**2) < 1e-6 

61 

62 positions = molecule.get_positions() 

63 # x-coordinate x-coordinate 

64 # v v 

65 d0 = positions[1, 0] - positions[0, 0] 

66 # ^ ^ 

67 # second atom first atom 

68 

69 print('experimental bond length:') 

70 print('hydrogen molecule energy: %7.3f eV' % e1) 

71 print('bondlength : %7.3f Ang' % d0) 

72 

73 # Find the theoretical bond length: 

74 relax = BFGS(molecule) 

75 relax.run(fmax=0.05) 

76 

77 e2 = molecule.get_potential_energy() 

78 

79 positions = molecule.get_positions() 

80 # x-coordinate x-coordinate 

81 # v v 

82 d0 = positions[1, 0] - positions[0, 0] 

83 # ^ ^ 

84 # second atom first atom 

85 

86 print('PBE energy minimum:') 

87 print('hydrogen molecule energy: %7.3f eV' % e2) 

88 print('bondlength : %7.3f Ang' % d0) 

89 

90 molecule = GPAW('H2fa.gpw', txt='H2.txt').get_atoms() 

91 relax = BFGS(molecule) 

92 relax.run(fmax=0.05) 

93 e2q = molecule.get_potential_energy() 

94 positions = molecule.get_positions() 

95 d0q = positions[1, 0] - positions[0, 0] 

96 assert abs(e2 - e2q) < 2e-6 

97 assert abs(d0q - d0) < 4e-4 

98 

99 f0 = molecule.get_forces() 

100 del relax, molecule 

101 

102 world.barrier() # syncronize before reading text output file 

103 f = read('H2.txt').get_forces() 

104 assert abs(f - f0).max() < 5e-6 # 5 digits in txt file 

105 

106 energy_tolerance = 0.002 

107 assert e1 == pytest.approx(-6.287873, abs=energy_tolerance) 

108 assert e2 == pytest.approx(-6.290744, abs=energy_tolerance) 

109 assert e2q == pytest.approx(-6.290744, abs=energy_tolerance)