1import pytest 2from gpaw.mpi import world 3from ase import Atoms, Atom 4from ase.build import molecule 5from ase.parallel import barrier 6from ase.units import Hartree, mol, kcal 7from gpaw import GPAW 8from gpaw.mixer import Mixer, MixerSum 9from gpaw.occupations import FermiDirac 10from gpaw.test import gen 11from gpaw.mpi import rank 12 13pytestmark = pytest.mark.skipif(world.size < 4, 14 reason='world.size < 4') 15 16 17@pytest.mark.slow 18def test_exx_AA_enthalpy(in_tmp_dir): 19 data = {} 20 21 def _xc(name): 22 return {'name': name, 'stencil': 1} 23 24 # data (from tables.pdf of 10.1063/1.1626543) 25 data['N'] = { 26 # intermolecular distance (A), 27 # formation enthalpy(298) (kcal/mol) on B3LYP geometry 28 'exp': (1.098, 0.0, 'none', 'none'), 29 'HCTH407': (1.097, 7.9, 'HCTH407', 'gga'), 30 'PBE': (1.103, -15.1, 'PBE', 'gga'), 31 'BLYP': (1.103, -11.7, 'BLYP', 'gga'), 32 'BP86': (1.104, -15.6, 'BP86', 'gga'), 33 'BPW91': (1.103, -8.5, 'BPW91', 'gga'), 34 'B3LYP': (1.092, -1.03, 'BLYP', 'hyb_gga'), 35 'B3PW91': (1.091, 2.8, 'PW91', 'hyb_gga'), 36 'PBE0': (1.090, 3.1, 'PBE', 'hyb_gga'), 37 'PBEH': (1.090, 3.1, 'PBE', 'hyb_gga'), 38 'magmom': 3.0, 39 # tables.pdf: 40 # http://ftp.aip.org/epaps/journ_chem_phys/E-JCPSA6-119-302348/tables.pdf 41 'R_AA_B3LYP': 1.092, # (from tables.pdf of 10.1063/1.1626543) (Ang) 42 'ZPE_AA_B3LYP': 0.005457 * Hartree, # (from benchmarks.txt of 43 # 10.1063/1.1626543) (eV) 44 'H_298_H_0_AA_B3LYP': 0.003304 * Hartree, # (from benchmarks.txt of 45 # 10.1063/1.1626543) (eV) 46 'H_298_H_0_A': 1.04 / (mol / kcal), # (from 10.1063/1.473182) (eV) 47 'dHf_0_A': 112.53 / (mol / kcal)} # (from 10.1063/1.473182) (eV) 48 49 data['O'] = { 50 # intermolecular distance (A), 51 # formation enthalpy(298) (kcal/mol) on B3LYP geometry 52 'exp': (1.208, 0.0, 'none', 'none'), 53 'HCTH407': (1.202, -14.5, 'HCTH407', 'gga'), 54 'PBE': (1.218, -23.6, 'PBE', 'gga'), 55 'BLYP': (1.229, -15.4, 'BLYP', 'gga'), 56 'BP86': (1.220, -21.9, 'BP86', 'gga'), 57 'BPW91': (1.219, -17.9, 'BPW91', 'gga'), 58 'B3LYP': (1.204, -3.7, 'BLYP', 'hyb_gga'), 59 'B3PW91': (1.197, -5.1, 'PW91', 'hyb_gga'), 60 'PBE0': (1.192, -4.3, 'PBE', 'hyb_gga'), 61 'PBEH': (1.192, -4.3, 'PBE', 'hyb_gga'), 62 'magmom': 2.0, 63 # tables.pdf: 64 # http://ftp.aip.org/epaps/journ_chem_phys/E-JCPSA6-119-302348/tables.pdf 65 'R_AA_B3LYP': 1.204, # (from tables.pdf of 10.1063/1.1626543) (Ang) 66 'ZPE_AA_B3LYP': 0.003736 * Hartree, # (from benchmarks.txt of 67 # 10.1063/1.1626543) (eV) 68 'H_298_H_0_AA_B3LYP': 0.003307 * Hartree, 69 # (from benchmarks.txt of 10.1063/1.1626543) (eV) 70 'H_298_H_0_A': 1.04 / (mol / kcal), # (from 10.1063/1.473182) (eV) 71 'dHf_0_A': 58.99 / (mol / kcal)} # (from 10.1063/1.473182) (eV) 72 73 data['H'] = { 74 # intermolecular distance (A), 75 # formation enthalpy(298) (kcal/mol) on B3LYP geometry 76 'exp': (0.741, 0.0, 'none', 'none'), 77 'HCTH407': (0.744, 1.8, 'HCTH407', 'gga'), 78 'PBE': (0.750, 5.1, 'PBE', 'gga'), 79 'BLYP': (0.746, 0.3, 'BLYP', 'gga'), 80 'BP86': (0.750, -1.8, 'BP86', 'gga'), 81 'BPW91': (0.748, 4.0, 'BPW91', 'gga'), 82 'B3LYP': (0.742, -0.5, 'BLYP', 'hyb_gga'), 83 'B3PW91': (0.744, 2.4, 'PW91', 'hyb_gga'), 84 'PBE0': (0.745, 5.3, 'PBE', 'hyb_gga'), 85 'PBEH': (0.745, 5.3, 'PBE', 'hyb_gga'), 86 'magmom': 1.0, 87 # tables.pdf: 88 # http://ftp.aip.org/epaps/journ_chem_phys/E-JCPSA6-119-302348/tables.pdf 89 'R_AA_B3LYP': 0.742, # (from tables.pdf of 10.1063/1.1626543) (Ang) 90 'ZPE_AA_B3LYP': 0.010025 * Hartree, # (from benchmarks.txt of 91 # 10.1063/1.1626543) (eV) 92 'H_298_H_0_AA_B3LYP': 0.003305 * Hartree, # (from benchmarks.txt of 93 # 10.1063/1.1626543) (eV) 94 'H_298_H_0_A': 1.01 / (mol / kcal), # (from 10.1063/1.473182) (eV) 95 'dHf_0_A': 51.63 / (mol / kcal)} # (from 10.1063/1.473182) (eV) 96 97 def calculate(element, h, vacuum, xc, magmom): 98 99 atom = Atoms([Atom(element, (0, 0, 0))]) 100 if magmom > 0.0: 101 mms = [magmom for i in range(len(atom))] 102 atom.set_initial_magnetic_moments(mms) 103 104 atom.center(vacuum=vacuum) 105 106 mixer = MixerSum(beta=0.4) 107 if element == 'O': 108 mixer = MixerSum(0.4, nmaxold=1, weight=100) 109 atom.set_positions(atom.get_positions() + [0.0, 0.0, 0.0001]) 110 111 calc_atom = GPAW(h=h, 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='%s.%s.txt' % (element, xc)) 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(h=h, xc=_xc(data[element][xc][2]), 129 experimental={'niter_fixdensity': 2}, 130 eigensolver='rmm-diis', 131 mixer=mixer, 132 parallel=dict(augment_grids=True), 133 txt='%s2.%s.txt' % (element, xc)) 134 compound.set_distance(0, 1, data[element]['R_AA_B3LYP']) 135 compound.center(vacuum=vacuum) 136 137 compound.calc = calc 138 139 if data[element][xc][3] == 'hyb_gga': # only for hybrids 140 e_atom = atom.get_potential_energy() 141 e_compound = compound.get_potential_energy() 142 143 calc_atom.set(xc=_xc(xc)) 144 calc.set(xc=_xc(xc)) 145 146 e_atom = atom.get_potential_energy() 147 e_compound = compound.get_potential_energy() 148 149 dHf_0 = (e_compound - 2 * e_atom + data[element]['ZPE_AA_B3LYP'] + 150 2 * data[element]['dHf_0_A']) 151 dHf_298 = (dHf_0 + data[element]['H_298_H_0_AA_B3LYP'] - 152 2 * data[element]['H_298_H_0_A']) * (mol / kcal) 153 de = dHf_298 - data[element][xc][1] 154 E[element][xc] = de 155 if rank == 0: 156 print((xc, h, 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 elif element == 'O': 161 assert dHf_298 == pytest.approx(data[element][xc][1], abs=7.5) 162 else: 163 assert dHf_298 == pytest.approx(data[element][xc][1], abs=2.15) 164 assert de == pytest.approx(E_ref[element][xc], abs=0.06) 165 166 E = {} 167 168 E_ref = {'H': {'HCTH407': 0.19286893273630645, 169 'B3LYP': -0.11369634560501423, 170 'PBE0': -0.21413764474738262, 171 'PBEH': -0.14147808591211231}, 172 'N': {'HCTH407': 2.1354017840869268, 173 'B3LYP': 0.63466589919873972, 174 'PBE0': -0.33376468078480226, 175 'PBEH': -0.30365500626180042}} # svnversion 5599 # -np 4 176 177 for element in ['H']: # , 'N']:#, 'O']: # oxygen atom fails to converge 178 E[element] = {} 179 for xc in ['HCTH407', 'PBE0', 'B3LYP']: 180 setup = data[element][xc][2] 181 if data[element][xc][3] == 'hyb_gga': # only for hybrids 182 exx = True 183 else: 184 exx = False 185 gen(element, exx=exx, xcname=setup, write_xml=True) 186 for h in [0.20]: 187 for vacuum in [4.5]: 188 calculate(element, h, vacuum, xc, data[element]['magmom']) 189 barrier() 190