Coverage for gpaw/test/vdw/test_ts09.py: 100%

1import ase.io 

2import numpy as np 

3import pytest 

4from ase.build import molecule 

5from ase.calculators.vdwcorrection import vdWTkatchenko09prl 

6from ase.parallel import barrier, parprint 

7 

8from gpaw import GPAW 

9from gpaw.analyse.hirshfeld import HirshfeldPartitioning 

10from gpaw.analyse.vdwradii import vdWradii 

11from gpaw.utilities.adjust_cell import adjust_cell 

12 

13 

14@pytest.mark.old_gpaw_only 

15def test_vdw_ts09(in_tmp_dir): 

16 h = 0.4 

17 s = molecule('LiH') 

18 adjust_cell(s, 3., h=h) 

19 

20 def print_charge_and_check(hp, q=0, label='unpolarized'): 

21 q_a = np.array(hp.get_charges()) 

22 parprint('Charges ({0})='.format(label), q_a, ', sum=', q_a.sum()) 

23 assert q_a.sum() == pytest.approx(q, abs=0.03) 

24 return q_a 

25 

26 # spin unpolarized 

27 

28 if 1: 

29 out_traj = 'LiH.traj' 

30 out_txt = 'LiH.txt' 

31 

32 cc = GPAW(mode='fd', h=h, xc='PBE', txt=out_txt) 

33 

34 # this is needed to initialize txt output 

35 cc.initialize(s) 

36 

37 hp = HirshfeldPartitioning(cc) 

38 c = vdWTkatchenko09prl(hp, 

39 vdWradii(s.get_chemical_symbols(), 'PBE')) 

40 s.calc = c 

41 E = s.get_potential_energy() 

42 F_ac = s.get_forces() 

43 s.write(out_traj) 

44 q_a = print_charge_and_check(hp) 

45 

46 barrier() 

47 

48 # test I/O, accuracy due to text output 

49 accuracy = 1.e-5 

50 for fname in [out_traj, out_txt]: 

51 s_out = ase.io.read(fname) 

52 assert s_out.get_potential_energy() == pytest.approx(E, 

53 abs=accuracy) 

54 for fi, fo in zip(F_ac, s_out.get_forces()): 

55 assert fi == pytest.approx(fo, abs=accuracy) 

56 

57 # spin polarized 

58 

59 if 0: 

60 ccs = GPAW(mode='fd', h=h, xc='PBE', spinpol=True, 

61 txt=None) 

62 hps = HirshfeldPartitioning(ccs) 

63 cs = vdWTkatchenko09prl(hps, vdWradii(s.get_chemical_symbols(), 'PBE')) 

64 s.calc = cs 

65 Es = s.get_potential_energy() 

66 Fs_ac = s.get_forces() 

67 

68 qs_a = print_charge_and_check(hps, label='spin') 

69 

70 assert q_a == pytest.approx(qs_a, abs=1.e-6) 

71 assert E == pytest.approx(Es, abs=1.e-4) 

72 assert F_ac == pytest.approx(Fs_ac, abs=1.e-4) 

73 

74 # charged 

75 

76 if 0: 

77 cc = cc.new(charge=1) 

78 hpp = HirshfeldPartitioning(cc) 

79 cp = vdWTkatchenko09prl(hpp, 

80 vdWradii(s.get_chemical_symbols(), 'PBE')) 

81 s.calc = cp 

82 E = s.get_potential_energy() 

83 

84 print_charge_and_check(hpp, 1, label='+1')