Coverage for gpaw/elph/displacements.py: 92%

1r"""Finite difference calculation of changes in potential. 

2 

3The implementation is based on finite-difference calculations of the the atomic 

4gradients of the effective potential expressed on a real-space grid. The el-ph 

5couplings are obtained from LCAO representations of the atomic gradients of the 

6effective potential and the electronic states. 

7 

8In PAW the matrix elements of the derivative of the effective potential is 

9given by the sum of the following contributions:: 

10 

11 d d 

12 < i | -- V | j > = < i | -- V | j> 

13 du eff du 

14 

15 _ 

16 \ ~a d . ~a 

17 + ) < i | p > -- /_\H < p | j > 

18 /_ i du ij j 

19 a,ij 

20 

21 _ 

22 \ d ~a . ~a 

23 + ) < i | -- p > /_\H < p | j > 

24 /_ du i ij j 

25 a,ij 

26 

27 _ 

28 \ ~a . d ~a 

29 + ) < i | p > /_\H < -- p | j > 

30 /_ i ij du j 

31 a,ij 

32 

33where the first term is the derivative of the potential (Hartree + XC) and the 

34last three terms originate from the PAW (pseudopotential) part of the effective 

35DFT Hamiltonian. 

36 

37Note: It is a bit difficult to find good references for how spin-polarisation 

38is supposed to be handled. Here we just handle the spin channels separately. 

39Use with care. 

40 

41""" 

42from __future__ import annotations 

43import numpy as np 

44 

45from ase import Atoms 

46from ase.phonons import Displacement 

47 

48from gpaw.calculator import GPAW as OldGPAW 

49from gpaw.new.ase_interface import ASECalculator 

50from gpaw.utilities import pack_hermitian 

51 

52dr_version = 1 

53# v1: saves natom, supercell, delta 

54 

55 

56class DisplacementRunner(Displacement): 

57 """Class for calculating the changes in effective potential. 

58 

59 The derivative of the effective potential wrt atomic displacements is 

60 obtained from a finite difference approximation to the derivative by doing 

61 a self-consistent calculation for atomic displacements in the +/- 

62 directions. These calculations are carried out in the ``run`` member 

63 function. 

64 """ 

65 

66 def __init__(self, 

67 atoms: Atoms, 

68 calc: OldGPAW | ASECalculator, 

69 supercell: tuple = (1, 1, 1), 

70 name: str = 'elph', 

71 delta: float = 0.01, 

72 calculate_forces: bool = True) -> None: 

73 """Initialize with base class args and kwargs. 

74 

75 Parameters 

76 ---------- 

77 atoms: Atoms 

78 The atoms to work on. Primitive cell. 

79 calc: GPAW 

80 Calculator for the supercell finite displacement calculation. 

81 supercell: tuple, list 

82 Size of supercell given by the number of repetitions (l, m, n) of 

83 the small unit cell in each direction. 

84 name: str 

85 Name to use for files (default: 'elph'). 

86 delta: float 

87 Magnitude of displacements. (default: 0.01 A) 

88 calculate_forces: bool 

89 If true, also calculate and store the dynamical matrix. 

90 """ 

91 

92 # Init base class and make the center cell in the supercell the 

93 # reference cell 

94 Displacement.__init__(self, atoms, calc=calc, supercell=supercell, 

95 name=name, delta=delta, center_refcell=True) 

96 self.calculate_forces = calculate_forces 

97 

98 def calculate(self, atoms_N: Atoms, disp): 

99 return self(atoms_N) 

100 

101 def __call__(self, atoms_N: Atoms) -> dict: 

102 """Extract effective potential and projector coefficients.""" 

103 

104 # Do calculation 

105 atoms_N.get_potential_energy() 

106 

107 # Calculate forces if desired 

108 if self.calculate_forces: 

109 forces = atoms_N.get_forces() 

110 else: 

111 forces = None 

112 

113 # Get calculator 

114 calc = atoms_N.calc 

115 if not isinstance(calc, (OldGPAW, ASECalculator)): 

116 calc = calc.dft # unwrap DFTD3 wrapper 

117 

118 # Effective potential (in Hartree) and projector coefficients 

119 # Note: Need to use coarse grid, because we project onto basis later 

120 if isinstance(calc, ASECalculator): 

121 potential = calc.dft.potential 

122 Vt_sG = potential.vt_sR.gather(broadcast=True).data 

123 dH_all_asp = {a: pack_hermitian(dH_ii) 

124 for a, dH_ii 

125 in potential.dH_asii.gather(broadcast=True).items()} 

126 else: 

127 Vt_sG = calc.hamiltonian.vt_sG 

128 Vt_sG = calc.wfs.gd.collect(Vt_sG, broadcast=True) 

129 dH_asp = calc.hamiltonian.dH_asp 

130 

131 setups = calc.wfs.setups 

132 nspins = calc.wfs.nspins 

133 

134 dH_all_asp = {} 

135 for a, setup in enumerate(setups): 

136 ni = setup.ni 

137 nii = ni * (ni + 1) // 2 

138 dH_tmp_sp = np.zeros((nspins, nii)) 

139 if a in dH_asp: 

140 dH_tmp_sp[:] = dH_asp[a] 

141 calc.wfs.gd.comm.sum(dH_tmp_sp) 

142 dH_all_asp[a] = dH_tmp_sp 

143 

144 output = {'Vt_sG': Vt_sG, 'dH_all_asp': dH_all_asp} 

145 if forces is not None: 

146 output['forces'] = forces 

147 return output 

148 

149 def save_info(self) -> None: 

150 with self.cache.lock('info') as handle: 

151 if handle is not None: 

152 info = {'natom': len(self.atoms), 'supercell': self.supercell, 

153 'delta': self.delta, 'dr_version': dr_version} 

154 handle.save(info) 

155 

156 def run(self) -> None: 

157 """Run the calculations for the required displacements.""" 

158 # Save some information about this run 

159 self.save_info() 

160 Displacement.run(self)