Coverage for gpaw/analyse/vdwradii.py: 100%

1from ase.data import atomic_numbers, chemical_symbols 

2from ase.units import Bohr 

3 

4from gpaw.setup import Setups 

5from gpaw.xc import XC 

6from gpaw.mpi import world 

7 

8Bondi64jpc_vdWradii = { # units Anstrom 

9 'He': 1.40, 

10 'Ne': 1.54, 

11 'Ar': 1.88, 

12 'Kr': 2.02, 

13 'Xe': 2.16 

14} 

15# Van der Waals Radii after 

16# Pekka Pyykko, Chem. Rev. 97 (1997) 597-636 

17# Table 2 

18Pyykko97cr_vdWradii = { # units Anstrom 

19 'Ne': 1.55, 

20 'Ar': 1.88, 

21 'Kr': 2.00, 

22 'Xe': 2.18, 

23 'Rn': 2.24 

24} 

25collected_vdWradii = Bondi64jpc_vdWradii 

26collected_vdWradii['Rn'] = Pyykko97cr_vdWradii['Rn'] 

27 

28 

29def vdWradii(symbols, xc): 

30 """Find the elements van der Waals radius. 

31 

32 Method proposed in: 

33 Tkatchenko and Scheffler PRL 102 (2009) 073005 

34 

35 The returned radii are given in Angstroms. 

36 """ 

37 Z_rare_gas = [atomic_numbers[symbol] for symbol in Bondi64jpc_vdWradii] 

38 Z_rare_gas.append(atomic_numbers['Rn']) 

39 Z_rare_gas.sort() 

40 

41 if isinstance(xc, str): 

42 xc = XC(xc) 

43 

44 def get_density(Z): 

45 """Return density and radial grid from setup.""" 

46 # load setup 

47 setups = Setups([Z], 'paw', {}, xc, world=world) 

48 setup = setups[0].data 

49 # create density 

50 n_g = setup.nc_g.copy() 

51 for f, phi_g in zip(setup.f_j, setup.phi_jg): 

52 n_g += f * phi_g ** 2 

53 return n_g, setup.rgd.r_g 

54 

55 radii = [] 

56 radius = {} 

57 for symbol in symbols: 

58 Z = atomic_numbers[symbol] 

59 if symbol not in radius: 

60 # find the rare gas of the elements row 

61 Zrg = None 

62 for Zr in Z_rare_gas: 

63 if Zrg is None and Z <= Zr: 

64 Zrg = Zr 

65 

66 n_g, r_g = get_density(Zrg) 

67 # find density at R 

68 R = collected_vdWradii[chemical_symbols[Zrg]] / Bohr 

69 n = 0 

70 while r_g[n] < R: 

71 n += 1 

72 # linear interpolation 

73 ncut = (n_g[n - 1] + 

74 (n_g[n] - n_g[n - 1]) * (R - r_g[n - 1]) / 

75 (r_g[n] - r_g[n - 1])) 

76# print "Z, Zrg, ncut", Z, Zrg, ncut 

77 

78 # find own R at this density 

79 n_g, r_g = get_density(Z) 

80 n = 0 

81 while n_g[n] > ncut: 

82 n += 1 

83 # linear interpolation 

84 R = (r_g[n - 1] + 

85 (r_g[n] - r_g[n - 1]) * (ncut - n_g[n - 1]) / 

86 (n_g[n] - n_g[n - 1])) 

87 radius[symbol] = R * Bohr 

88 

89 radii.append(radius[symbol]) 

90 

91 return radii