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