Coverage for gpaw/test/xc/test_beef.py: 100%

1import warnings 

2import pytest 

3import numpy as np 

4from ase.build import bulk 

5from ase.dft.bee import BEEFEnsemble, readbee 

6from gpaw import GPAW, Mixer, PW 

7from gpaw.test import gen 

8from gpaw.mpi import world 

9import gpaw.cgpaw as cgpaw 

10 

11 

12@pytest.mark.mgga 

13@pytest.mark.libxc 

14@pytest.mark.slow 

15@pytest.mark.parametrize('xc', ['mBEEF', 'BEEF-vdW', 'mBEEF-vdW']) 

16def test_beef(in_tmp_dir, xc, gpaw_new): 

17 if xc[0] == 'm': 

18 assert cgpaw.lxcXCFuncNum('MGGA_X_MBEEF') is not None 

19 

20 results = {'mBEEF': (5.449, 0.056), 

21 'BEEF-vdW': (5.484, 0.071), 

22 'mBEEF-vdW': (5.426, 0.025)} 

23 

24 kwargs = dict() 

25 if xc == 'mBEEF-vdW': 

26 kwargs['setups'] = dict(Si=gen('Si', xcname='PBEsol')) 

27 

28 E = [] 

29 V = [] 

30 for a in np.linspace(5.4, 5.5, 3): 

31 si = bulk('Si', a=a) 

32 si.calc = GPAW(txt='Si-' + xc + '.txt', 

33 mixer=Mixer(0.8, 7, 50.0), 

34 xc=xc, 

35 kpts=[2, 2, 2], 

36 mode=PW(200), 

37 **kwargs) 

38 with warnings.catch_warnings(): 

39 warnings.filterwarnings('ignore', module='gpaw.xc.libxc') 

40 E.append(si.get_potential_energy()) 

41 ens = BEEFEnsemble(si.calc, verbose=False) 

42 ens.get_ensemble_energies(200) 

43 ens.write('Si-{}-{:.3f}'.format(xc, a)) 

44 V.append(si.get_volume()) 

45 p = np.polyfit(V, E, 2) 

46 v0 = np.roots(np.polyder(p))[0] 

47 a = (v0 * 4)**(1 / 3) 

48 

49 a0, da0 = results[xc] 

50 

51 assert abs(a - a0) < 0.002, (xc, a, a0) 

52 

53 if world.rank == 0: 

54 E = [] 

55 for a in np.linspace(5.4, 5.5, 3): 

56 e = readbee('Si-{}-{:.3f}'.format(xc, a)) 

57 E.append(e) 

58 

59 A = [] 

60 for energies in np.array(E).T: 

61 p = np.polyfit(V, energies, 2) 

62 assert p[0] > 0, (V, E, p) 

63 v0 = np.roots(np.polyder(p))[0] 

64 A.append((v0 * 4)**(1 / 3)) 

65 

66 A = np.array(A) 

67 a = A.mean() 

68 da = A.std() 

69 

70 print('a(ref) = {:.3f} +- {:.3f}'.format(a0, da0)) 

71 print('a = {:.3f} +- {:.3f}'.format(a, da)) 

72 assert abs(a - a0) < 0.01 

73 assert abs(da - da0) < 0.01