1# XCFun, an arbitrary order exchange-correlation library
2# Copyright (C) 2020 Ulf Ekström and contributors.
3#
4# This file is part of XCFun.
5#
6# This Source Code Form is subject to the terms of the Mozilla Public
7# License, v. 2.0. If a copy of the MPL was not distributed with this
8# file, You can obtain one at http://mozilla.org/MPL/2.0/.
9#
10# For information on the complete list of contributors to the
11# XCFun library, see: <https://xcfun.readthedocs.io/>
12
13import pytest
14
15import numpy
16from numpy.testing import assert_allclose
17
18import xcfun
19
20
21def test_splash():
22
23    version = xcfun.xcfun_version()
24    splash = xcfun.xcfun_splash()
25
26    print("\n Sample use of python interface to xcfun")
27
28    print("\n\nXCFun version ", version)
29    print(splash)
30
31
32def test_lda_old():
33    fun_xc_name = 'LDA'
34    fun_xc_weight = 1.0
35    fun_xc = xcfun.Functional({fun_xc_name: fun_xc_weight})
36
37    n = numpy.array([1.0])
38    result = fun_xc.eval_potential_n(n)
39
40    energy = result[0, 0]
41    potential = result[0, 1]
42
43    print("Functional used: ", fun_xc_name, " weight ", fun_xc_weight)
44    print("        density: ", n)
45    print("         energy: ", energy)
46    print("      potential: ", potential)
47
48    assert_allclose(energy, -0.8101513)
49    assert_allclose(potential, -1.06468341)
50
51
52@pytest.fixture
53def lda_fun():
54    return xcfun.Functional({'LDA': 1.0})
55
56
57@pytest.fixture
58def pbe_fun():
59    return xcfun.Functional({'PBE': 1.0})
60
61
62def test_type(lda_fun, pbe_fun):
63    assert lda_fun.type == 0
64    assert pbe_fun.type == 1
65
66
67@pytest.fixture
68def dens():
69    return numpy.arange(0.1, 5.0, 0.5)
70
71
72@pytest.fixture
73def spindens(dens):
74    return 0.1 * numpy.sin(dens)
75
76
77@pytest.fixture
78def densgrad(dens):
79    dg = numpy.zeros((dens.shape[0], 3))
80    dg[:, 0] = 0.2 * numpy.cos(dens)
81    dg[:, 1] = 0.3 * numpy.sin(dens)
82    dg[:, 2] = 0.1 * numpy.sin(dens)**2
83    return dg
84
85
86@pytest.fixture
87def denshess(dens):
88    dh = numpy.zeros((dens.shape[0], 6))
89    dh[:, 0] = 0.5
90    dh[:, 1] = 0.1 * dens
91    dh[:, 2] = numpy.cos(dens)**2
92    dh[:, 3] = numpy.sin(dens)
93    dh[:, 4] = 0.6 * numpy.cos(dens)**2
94    dh[:, 5] = 0.3 * numpy.sin(dens)
95    return dh
96
97
98@pytest.fixture
99def refen_lda():
100    return numpy.array([
101        -0.03962059, -0.41418026, -0.91826807, -1.50300217, -2.14968647, -2.84774973, -3.59023064, -4.37213299,
102        -5.18966423, -6.03982917
103    ])
104
105
106@pytest.fixture
107def refpot_lda():
108    return numpy.array([
109        -0.51789018, -0.90610372, -1.09730105, -1.23582171, -1.34751238, -1.44246858, -1.52581304, -1.60054497,
110        -1.66858921, -1.73126329
111    ])
112
113
114def test_lda_energy_restr(lda_fun, dens, refen_lda):
115    out = lda_fun.eval_energy_n(dens)
116    assert_allclose(out, refen_lda)
117
118
119def test_lda_potential_restr(lda_fun, dens, refen_lda, refpot_lda):
120    out = lda_fun.eval_potential_n(dens)
121    assert_allclose(out[:, 0], refen_lda)
122    assert_allclose(out[:, 1], refpot_lda)
123
124
125def test_lda_energy_unrestr_closedshell(lda_fun, dens, refen_lda):
126    rho = numpy.zeros((dens.shape[0], 2))
127    rho[:, 0] = 0.5 * dens[:]
128    rho[:, 1] = 0.5 * dens[:]
129
130    out = lda_fun.eval_energy_ab(rho)
131    assert_allclose(out, refen_lda)
132
133
134def test_lda_potential_unrestr_closedshell(lda_fun, dens, refen_lda, refpot_lda):
135    rho = numpy.zeros((dens.shape[0], 2))
136    rho[:, 0] = 0.5 * dens[:]
137    rho[:, 1] = 0.5 * dens[:]
138
139    out = lda_fun.eval_potential_ab(rho)
140    assert_allclose(out[:, 0], refen_lda)
141    assert_allclose(out[:, 1], refpot_lda)
142    assert_allclose(out[:, 2], refpot_lda)
143
144
145def test_lda_energy_unrestr_openshell(lda_fun, dens, spindens):
146    rho = numpy.zeros((dens.shape[0], 2))
147    rho[:, 0] = 0.5 * dens[:] + spindens[:]
148    rho[:, 1] = 0.5 * dens[:] - spindens[:]
149
150    out = lda_fun.eval_energy_ab(rho)
151
152    refen_lda_openshell = numpy.array([
153        -0.03985427, -0.41666061, -0.92248897, -1.50718908, -2.15231233, -2.84856692, -3.59023539, -4.3726223,
154        -5.19120401, -6.04193876
155    ])
156
157    assert_allclose(out, refen_lda_openshell)
158
159
160def test_lda_potential_unrestr_openshell(lda_fun, dens, spindens):
161    rho = numpy.zeros((dens.shape[0], 2))
162    rho[:, 0] = 0.5 * dens[:] + spindens[:]
163    rho[:, 1] = 0.5 * dens[:] - spindens[:]
164
165    out = lda_fun.eval_potential_ab(rho)
166
167    refpot_lda_alpha = numpy.array([
168        -0.53994997, -0.94756428, -1.14234107, -1.27610516, -1.37715554, -1.4581231, -1.52695501, -1.58940021,
169        -1.64952872, -1.70973604
170    ])
171    refpot_lda_beta = numpy.array([
172        -0.49296826, -0.85942489, -1.04738939, -1.1922106, -1.31627816, -1.42641378, -1.52466911, -1.61151607,
173        -1.68716914, -1.75220302
174    ])
175
176    assert_allclose(out[:, 1], refpot_lda_alpha)
177    assert_allclose(out[:, 2], refpot_lda_beta)
178
179    outen = lda_fun.eval_energy_ab(rho)
180    assert_allclose(out[:, 0], outen)
181
182
183@pytest.fixture
184def refen_pbe():
185    return numpy.array([
186        -0.04051725, -0.41397904, -0.91782841, -1.50231227, -2.14873903, -2.84654021, -3.58875595, -4.37039097,
187        -5.18765335, -6.03754832
188    ])
189
190
191@pytest.fixture
192def refpot_pbe():
193    return numpy.array([
194        -0.48524858, -0.90474127, -1.09660442, -1.23524903, -1.34697602, -1.44193844, -1.52527933, -1.60000832,
195        -1.66805094, -1.73072351
196    ])
197
198
199def test_pbe_energy_restr(pbe_fun, dens, densgrad, refen_pbe):
200    out = pbe_fun.eval_energy_n(dens, densgrad)
201    assert_allclose(out, refen_pbe)
202
203
204def test_pbe_potential_restr(pbe_fun, dens, densgrad, denshess, refen_pbe, refpot_pbe):
205    out = pbe_fun.eval_potential_n(dens, densgrad, denshess)
206    assert_allclose(out[:, 0], refen_pbe)
207    assert_allclose(out[:, 1], refpot_pbe)
208
209
210def test_pbe_energy_unrestr_closedshell(pbe_fun, dens, densgrad, refen_pbe):
211    rho = numpy.zeros((dens.shape[0], 2))
212    rho[:, 0] = 0.5 * dens[:]
213    rho[:, 1] = 0.5 * dens[:]
214
215    rhograd = numpy.zeros((dens.shape[0], 3, 2))
216    rhograd[:, 0:3, 0] = 0.5 * densgrad[:, 0:3]
217    rhograd[:, 0:3, 1] = 0.5 * densgrad[:, 0:3]
218
219    out = pbe_fun.eval_energy_ab(rho, rhograd)
220    assert_allclose(out, refen_pbe)
221
222
223def test_pbe_potential_unrestr_closedshell(pbe_fun, dens, densgrad, denshess, refen_pbe, refpot_pbe):
224    rho = numpy.zeros((dens.shape[0], 2))
225    rho[:, 0] = 0.5 * dens[:]
226    rho[:, 1] = 0.5 * dens[:]
227
228    rhograd = numpy.zeros((dens.shape[0], 3, 2))
229    rhograd[:, 0:3, 0] = 0.5 * densgrad[:, 0:3]
230    rhograd[:, 0:3, 1] = 0.5 * densgrad[:, 0:3]
231
232    rhohess = numpy.zeros((dens.shape[0], 6, 2))
233    rhohess[:, 0:6, 0] = 0.5 * denshess[:, 0:6]
234    rhohess[:, 0:6, 1] = 0.5 * denshess[:, 0:6]
235
236    out = pbe_fun.eval_potential_ab(rho, rhograd, rhohess)
237    assert_allclose(out[:, 0], refen_pbe)
238    assert_allclose(out[:, 1], refpot_pbe)
239    assert_allclose(out[:, 2], refpot_pbe)
240
241
242def test_pbe_energy_unrestr_openshell(pbe_fun, dens, spindens, densgrad, refen_pbe):
243    rho = numpy.zeros((dens.shape[0], 2))
244    rho[:, 0] = dens[:] + spindens[:]
245    rho[:, 1] = dens[:] - spindens[:]
246
247    rhograd = numpy.zeros((dens.shape[0], 3, 2))
248    rhograd[:, 0:3, 0] = 0.4 * densgrad[:, 0:3]
249    rhograd[:, 0:3, 1] = 0.6 * densgrad[:, 0:3]
250
251    out = pbe_fun.eval_energy_ab(rho, rhograd)
252
253    refen_pbe_openshell = numpy.array([
254        -0.09848243, -1.03067727, -2.28717791, -3.74476386, -5.35679716, -7.09762348, -8.95059771, -10.90338775,
255        -12.94601448, -15.07029572
256    ])
257
258    assert_allclose(out, refen_pbe_openshell)
259
260
261def test_pbe_potential_unrestr_closedshell(pbe_fun, dens, spindens, densgrad, denshess, refen_pbe, refpot_pbe):
262    rho = numpy.zeros((dens.shape[0], 2))
263    rho[:, 0] = dens[:] + spindens[:]
264    rho[:, 1] = dens[:] - spindens[:]
265
266    rhograd = numpy.zeros((dens.shape[0], 3, 2))
267    rhograd[:, 0:3, 0] = 0.4 * densgrad[:, 0:3]
268    rhograd[:, 0:3, 1] = 0.6 * densgrad[:, 0:3]
269
270    rhohess = numpy.zeros((dens.shape[0], 6, 2))
271    rhohess[:, 0:6, 0] = 0.7 * denshess[:, 0:6]
272    rhohess[:, 0:6, 1] = 0.2 * denshess[:, 0:6]
273
274    out = pbe_fun.eval_potential_ab(rho, rhograd, rhohess)
275
276    refpot_pbe_alpha = numpy.array([
277        -0.62768401, -1.15201996, -1.39499424, -1.56595725, -1.69939803, -1.80975816, -1.90530719, -1.99151265,
278        -2.07206904, -2.1492557
279    ])
280    refpot_pbe_beta = numpy.array([
281        -0.64069084, -1.10200908, -1.33808862, -1.51488446, -1.66214897, -1.79045119, -1.90425337, -2.00565175,
282        -2.09589401, -2.17609771
283    ])
284    assert_allclose(out[:, 1], refpot_pbe_alpha)
285    assert_allclose(out[:, 2], refpot_pbe_beta)
286
287    outen = pbe_fun.eval_energy_ab(rho, rhograd)
288    assert_allclose(out[:, 0], outen)
289
290
291def test_raises_eval_energy_n(lda_fun, pbe_fun):
292    # correct shape for all arguments (LDA)
293    lda_fun.eval_energy_n(numpy.zeros((10, )))
294
295    with pytest.raises(xcfun.XCFunException, match="Wrong shape of density argument"):
296        lda_fun.eval_energy_n(numpy.zeros((10, 2)))
297
298    # correct shape for all arguments (PBE)
299    pbe_fun.eval_energy_n(numpy.zeros((10, )), numpy.zeros((10, 3)))
300
301    with pytest.raises(xcfun.XCFunException, match="Density gradient required"):
302        pbe_fun.eval_energy_n(numpy.zeros((10, )))
303    with pytest.raises(xcfun.XCFunException, match="Wrong shape of densgrad argument"):
304        pbe_fun.eval_energy_n(numpy.zeros((10, )), numpy.zeros((10, )))
305    with pytest.raises(xcfun.XCFunException, match="Wrong shape of densgrad argument"):
306        pbe_fun.eval_energy_n(numpy.zeros((10, )), numpy.zeros((10, 4)))
307    with pytest.raises(xcfun.XCFunException, match="Wrong shape of densgrad argument"):
308        pbe_fun.eval_energy_n(numpy.zeros((10, )), numpy.zeros((20, 3)))
309
310
311def test_raises_eval_potential_n(lda_fun, pbe_fun):
312    # correct shape for all arguments (LDA)
313    lda_fun.eval_potential_n(numpy.zeros((10, )))
314
315    with pytest.raises(xcfun.XCFunException, match="Wrong shape of density argument"):
316        lda_fun.eval_potential_n(numpy.zeros((10, 2)))
317
318    # correct shape for all arguments (PBE)
319    pbe_fun.eval_potential_n(numpy.zeros((10, )), numpy.zeros((10, 3)), numpy.zeros((10, 6)))
320
321    with pytest.raises(xcfun.XCFunException, match="Density gradient and Hessian required"):
322        pbe_fun.eval_potential_n(numpy.zeros((10, )))
323    with pytest.raises(xcfun.XCFunException, match="Density gradient and Hessian required"):
324        pbe_fun.eval_potential_n(numpy.zeros((10, )), numpy.zeros((10, 3)))
325    with pytest.raises(xcfun.XCFunException, match="Wrong shape of densgrad argument"):
326        pbe_fun.eval_potential_n(numpy.zeros((10, )), numpy.zeros((10, )), numpy.zeros((10, 6)))
327    with pytest.raises(xcfun.XCFunException, match="Wrong shape of densgrad argument"):
328        pbe_fun.eval_potential_n(numpy.zeros((10, )), numpy.zeros((10, 4)), numpy.zeros((10, 6)))
329    with pytest.raises(xcfun.XCFunException, match="Wrong shape of densgrad argument"):
330        pbe_fun.eval_potential_n(numpy.zeros((10, )), numpy.zeros((20, 3)), numpy.zeros((10, 6)))
331    with pytest.raises(xcfun.XCFunException, match="Wrong shape of denshess argument"):
332        pbe_fun.eval_potential_n(numpy.zeros((10, )), numpy.zeros((10, 3)), numpy.zeros((10, )))
333    with pytest.raises(xcfun.XCFunException, match="Wrong shape of denshess argument"):
334        pbe_fun.eval_potential_n(numpy.zeros((10, )), numpy.zeros((10, 3)), numpy.zeros((10, 3)))
335    with pytest.raises(xcfun.XCFunException, match="Wrong shape of denshess argument"):
336        pbe_fun.eval_potential_n(numpy.zeros((10, )), numpy.zeros((10, 3)), numpy.zeros((20, 6)))
337
338
339def test_raises_eval_energy_ab(lda_fun, pbe_fun):
340    # correct shape for all arguments (LDA)
341    lda_fun.eval_energy_ab(numpy.zeros((10, 2)))
342
343    with pytest.raises(xcfun.XCFunException, match="Wrong shape of density argument"):
344        lda_fun.eval_energy_ab(numpy.zeros((10, )))
345    with pytest.raises(xcfun.XCFunException, match="Wrong shape of density argument"):
346        lda_fun.eval_energy_ab(numpy.zeros((10, 4)))
347
348    # correct shape for all arguments (PBE)
349    pbe_fun.eval_energy_ab(numpy.zeros((10, 2)), numpy.zeros((10, 3, 2)))
350
351    with pytest.raises(xcfun.XCFunException, match="Density gradient required"):
352        pbe_fun.eval_energy_ab(numpy.zeros((10, 2)))
353    with pytest.raises(xcfun.XCFunException, match="Wrong shape of densgrad argument"):
354        pbe_fun.eval_energy_ab(numpy.zeros((10, 2)), numpy.zeros((10, 3)))
355    with pytest.raises(xcfun.XCFunException, match="Wrong shape of densgrad argument"):
356        pbe_fun.eval_energy_ab(numpy.zeros((10, 2)), numpy.zeros((10, 4, 2)))
357    with pytest.raises(xcfun.XCFunException, match="Wrong shape of densgrad argument"):
358        pbe_fun.eval_energy_ab(numpy.zeros((10, 2)), numpy.zeros((20, 3, 2)))
359
360
361def test_raises_eval_potential_ab(lda_fun, pbe_fun):
362    # correct shape for all arguments (LDA)
363    lda_fun.eval_potential_ab(numpy.zeros((10, 2)))
364
365    with pytest.raises(xcfun.XCFunException, match="Wrong shape of density argument"):
366        lda_fun.eval_potential_ab(numpy.zeros((10, )))
367    with pytest.raises(xcfun.XCFunException, match="Wrong shape of density argument"):
368        lda_fun.eval_potential_ab(numpy.zeros((10, 4)))
369
370    # correct shape for all arguments (PBE)
371    pbe_fun.eval_potential_ab(numpy.zeros((10, 2)), numpy.zeros((10, 3, 2)), numpy.zeros((10, 6, 2)))
372
373    with pytest.raises(xcfun.XCFunException, match="Density gradient and Hessian required"):
374        pbe_fun.eval_potential_ab(numpy.zeros((10, 2)))
375    with pytest.raises(xcfun.XCFunException, match="Density gradient and Hessian required"):
376        pbe_fun.eval_potential_ab(numpy.zeros((10, 2)), numpy.zeros((10, 3, 2)))
377    with pytest.raises(xcfun.XCFunException, match="Wrong shape of densgrad argument"):
378        pbe_fun.eval_potential_ab(numpy.zeros((10, 2)), numpy.zeros((10, 3)), numpy.zeros((10, 6, 2)))
379    with pytest.raises(xcfun.XCFunException, match="Wrong shape of densgrad argument"):
380        pbe_fun.eval_potential_ab(numpy.zeros((10, 2)), numpy.zeros((10, 4, 2)), numpy.zeros((10, 6, 2)))
381    with pytest.raises(xcfun.XCFunException, match="Wrong shape of densgrad argument"):
382        pbe_fun.eval_potential_ab(numpy.zeros((10, 2)), numpy.zeros((20, 3, 2)), numpy.zeros((10, 6, 2)))
383    with pytest.raises(xcfun.XCFunException, match="Wrong shape of denshess argument"):
384        pbe_fun.eval_potential_ab(numpy.zeros((10, 2)), numpy.zeros((10, 3, 2)), numpy.zeros((10, 6)))
385    with pytest.raises(xcfun.XCFunException, match="Wrong shape of denshess argument"):
386        pbe_fun.eval_potential_ab(numpy.zeros((10, 2)), numpy.zeros((10, 3, 2)), numpy.zeros((10, 3, 2)))
387    with pytest.raises(xcfun.XCFunException, match="Wrong shape of denshess argument"):
388        pbe_fun.eval_potential_ab(numpy.zeros((10, 2)), numpy.zeros((10, 3, 2)), numpy.zeros((20, 6, 2)))
389