Coverage for gpaw/test/test_AA_enthalpy.py: 32%

1import pytest 

2from gpaw.mpi import world 

3from ase import Atoms, Atom 

4from ase.build import molecule 

5from ase.units import Hartree, mol, kcal 

6from gpaw import GPAW 

7from gpaw.mixer import Mixer, MixerSum 

8from gpaw.occupations import FermiDirac 

9from gpaw.test import gen 

10 

11 

12pytestmark = pytest.mark.skipif(world.size < 4, 

13 reason='world.size < 4') 

14 

15 

16def _xc(name): 

17 return {'name': name, 'stencil': 1} 

18 

19 

20data = {} 

21 

22 

23# data (from tables.pdf of 10.1063/1.1626543) 

24data['N'] = { 

25 # intermolecular distance (A), 

26 # formation enthalpy(298) (kcal/mol) on B3LYP geometry 

27 'exp': (1.098, 0.0, 'none', 'none'), 

28 'PBE': (1.103, -15.1, 'PBE', 'gga'), 

29 'BLYP': (1.103, -11.7, 'BLYP', 'gga'), 

30 'BP86': (1.104, -15.6, 'BP86', 'gga'), 

31 'BPW91': (1.103, -8.5, 'BPW91', 'gga'), 

32 'B3LYP': (1.092, -1.03, 'BLYP', 'hyb_gga'), 

33 'B3PW91': (1.091, 2.8, 'PW91', 'hyb_gga'), 

34 'PBE0': (1.090, 3.1, 'PBE', 'hyb_gga'), 

35 'PBEH': (1.090, 3.1, 'PBE', 'hyb_gga'), 

36 'magmom': 3.0, 

37 # tables.pdf: 

38 # https://aip.scitation.org/doi/suppl/10.1063/1.1626543/suppl_file/tables.pdf 

39 'R_AA_B3LYP': 1.092, # (from tables.pdf of 10.1063/1.1626543) (Ang) 

40 'ZPE_AA_B3LYP': 0.005457 * Hartree, # (from benchmarks.txt of 

41 # 10.1063/1.1626543) (eV) 

42 'H_298_H_0_AA_B3LYP': 0.003304 * Hartree, # (from benchmarks.txt of 

43 # 10.1063/1.1626543) (eV) 

44 'H_298_H_0_A': 1.04 / (mol / kcal), # (from 10.1063/1.473182) (eV) 

45 'dHf_0_A': 112.53 / (mol / kcal)} # (from 10.1063/1.473182) (eV) 

46 

47 

48data['O'] = { 

49 # intermolecular distance (A), 

50 # formation enthalpy(298) (kcal/mol) on B3LYP geometry 

51 'exp': (1.208, 0.0, 'none', 'none'), 

52 'PBE': (1.218, -23.6, 'PBE', 'gga'), 

53 'BLYP': (1.229, -15.4, 'BLYP', 'gga'), 

54 'BP86': (1.220, -21.9, 'BP86', 'gga'), 

55 'BPW91': (1.219, -17.9, 'BPW91', 'gga'), 

56 'B3LYP': (1.204, -3.7, 'BLYP', 'hyb_gga'), 

57 'B3PW91': (1.197, -5.1, 'PW91', 'hyb_gga'), 

58 'PBE0': (1.192, -4.3, 'PBE', 'hyb_gga'), 

59 'PBEH': (1.192, -4.3, 'PBE', 'hyb_gga'), 

60 'magmom': 2.0, 

61 # tables.pdf: 

62 # https://aip.scitation.org/doi/suppl/10.1063/1.1626543/suppl_file/tables.pdf 

63 'R_AA_B3LYP': 1.204, # (from tables.pdf of 10.1063/1.1626543) (Ang) 

64 'ZPE_AA_B3LYP': 0.003736 * Hartree, # (from benchmarks.txt of 

65 # 10.1063/1.1626543) (eV) 

66 'H_298_H_0_AA_B3LYP': 0.003307 * Hartree, 

67 # (from benchmarks.txt of 10.1063/1.1626543) (eV) 

68 'H_298_H_0_A': 1.04 / (mol / kcal), # (from 10.1063/1.473182) (eV) 

69 'dHf_0_A': 58.99 / (mol / kcal)} # (from 10.1063/1.473182) (eV) 

70 

71 

72data['H'] = { 

73 # intermolecular distance (A), 

74 # formation enthalpy(298) (kcal/mol) on B3LYP geometry 

75 'exp': (0.741, 0.0, 'none', 'none'), 

76 'PBE': (0.750, 5.1, 'PBE', 'gga'), 

77 'BLYP': (0.746, 0.3, 'BLYP', 'gga'), 

78 'BP86': (0.750, -1.8, 'BP86', 'gga'), 

79 'BPW91': (0.748, 4.0, 'BPW91', 'gga'), 

80 'B3LYP': (0.742, -0.5, 'BLYP', 'hyb_gga'), 

81 'B3PW91': (0.744, 2.4, 'PW91', 'hyb_gga'), 

82 'PBE0': (0.745, 5.3, 'PBE', 'hyb_gga'), 

83 'PBEH': (0.745, 5.3, 'PBE', 'hyb_gga'), 

84 'magmom': 1.0, 

85 # tables.pdf: 

86 # https://aip.scitation.org/doi/suppl/10.1063/1.1626543/suppl_file/tables.pdf 

87 'R_AA_B3LYP': 0.742, # (from tables.pdf of 10.1063/1.1626543) (Ang) 

88 'ZPE_AA_B3LYP': 0.010025 * Hartree, # (from benchmarks.txt of 

89 # 10.1063/1.1626543) (eV) 

90 'H_298_H_0_AA_B3LYP': 0.003305 * Hartree, # (from benchmarks.txt of 

91 # 10.1063/1.1626543) (eV) 

92 'H_298_H_0_A': 1.01 / (mol / kcal), # (from 10.1063/1.473182) (eV) 

93 'dHf_0_A': 51.63 / (mol / kcal)} # (from 10.1063/1.473182) (eV) 

94 

95 

96def calculate(element, vacuum, xc, magmom): 

97 

98 atom = Atoms([Atom(element, (0, 0, 0))]) 

99 if magmom > 0.0: 

100 mms = [magmom for i in range(len(atom))] 

101 atom.set_initial_magnetic_moments(mms) 

102 

103 atom.center(vacuum=vacuum) 

104 

105 mixer = MixerSum(beta=0.4) 

106 if element == 'O': 

107 mixer = MixerSum(0.4, nmaxold=1, weight=100) 

108 atom.set_positions(atom.get_positions() + [0.0, 0.0, 0.0001]) 

109 

110 calc_atom = GPAW(mode='fd', 

111 xc=_xc(data[element][xc][2]), 

112 experimental={'niter_fixdensity': 2}, 

113 eigensolver='rmm-diis', 

114 occupations=FermiDirac(0.0, fixmagmom=True), 

115 mixer=mixer, 

116 parallel=dict(augment_grids=True), 

117 nbands=-2, 

118 txt=f'{element}.{xc}.txt') 

119 atom.calc = calc_atom 

120 

121 mixer = Mixer(beta=0.4, weight=100) 

122 compound = molecule(element + '2') 

123 if compound == 'O2': 

124 mixer = MixerSum(beta=0.4) 

125 mms = [1.0 for i in range(len(compound))] 

126 compound.set_initial_magnetic_moments(mms) 

127 

128 calc = GPAW(mode='fd', 

129 xc=_xc(data[element][xc][2]), 

130 experimental={'niter_fixdensity': 2}, 

131 eigensolver='rmm-diis', 

132 mixer=mixer, 

133 parallel=dict(augment_grids=True), 

134 txt=f'{element}2.{xc}.txt') 

135 compound.set_distance(0, 1, data[element]['R_AA_B3LYP']) 

136 compound.center(vacuum=vacuum) 

137 

138 compound.calc = calc 

139 

140 if data[element][xc][3] == 'hyb_gga': # only for hybrids 

141 e_atom = atom.get_potential_energy() 

142 e_compound = compound.get_potential_energy() 

143 

144 atom.calc = calc_atom.new(xc=_xc(xc)) 

145 compound.calc = calc.new(xc=_xc(xc)) 

146 

147 e_atom = atom.get_potential_energy() 

148 e_compound = compound.get_potential_energy() 

149 

150 dHf_0 = (e_compound - 2 * e_atom + data[element]['ZPE_AA_B3LYP'] + 

151 2 * data[element]['dHf_0_A']) 

152 dHf_298 = (dHf_0 + data[element]['H_298_H_0_AA_B3LYP'] - 

153 2 * data[element]['H_298_H_0_A']) * (mol / kcal) 

154 de = dHf_298 - data[element][xc][1] 

155 

156 print((xc, vacuum, dHf_298, data[element][xc][1], de, 

157 de / data[element][xc][1])) 

158 if element == 'H': 

159 assert dHf_298 == pytest.approx(data[element][xc][1], abs=0.25) 

160 

161 elif element == 'O': 

162 assert dHf_298 == pytest.approx(data[element][xc][1], abs=7.5) 

163 else: 

164 assert dHf_298 == pytest.approx(data[element][xc][1], abs=2.15) 

165 assert de == pytest.approx(E_ref[element][xc], abs=0.06) 

166 

167 

168E_ref = {'H': {'B3LYP': -0.11369634560501423, 

169 'PBE0': -0.21413764474738262, 

170 'PBEH': -0.14147808591211231}, 

171 'N': {'B3LYP': 0.63466589919873972, 

172 'PBE0': -0.33376468078480226, 

173 'PBEH': -0.30365500626180042}} # svnversion 5599 # -np 4 

174 

175 

176@pytest.mark.old_gpaw_only 

177@pytest.mark.slow 

178@pytest.mark.parametrize('xc', ['PBE0', 'B3LYP']) 

179def test_exx_AA_enthalpy(in_tmp_dir, add_cwd_to_setup_paths, xc): 

180 element = 'H' 

181 vacuum = 4.5 

182 

183 setup = data[element][xc][2] 

184 enable_exx = data[element][xc][3] == 'hyb_gga' # only for hybrids 

185 gen(element, exx=enable_exx, xcname=setup, write_xml=True) 

186 calculate(element, vacuum, xc, data[element]['magmom'])