Coverage for gpaw/lcao/el_ph.py: 8%
239 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-08 00:17 +0000
1import pickle
2import numpy as np
3from os.path import isfile
5from ase.parallel import parprint
6from ase.units import Bohr, Ha
8from gpaw.mpi import world
9from gpaw.utilities import unpack_hermitian
10from gpaw.lcao.projected_wannier import dots
11from gpaw.utilities.tools import tri2full
12from gpaw.lfc import LocalizedFunctionsCollection as LFC
15"""This module is used to calculate the electron-phonon coupling matrix,
16 expressed in terms of GPAW LCAO orbitals.
18 This module is not maintained and possibly broken.
19 Use gpaw/elph/* instead.
20 """
23class ElectronPhononCouplingMatrix:
24 r"""Class for calculating the electron-phonon coupling matrix, defined
25 by the electron phonon interaction.
27 ::
29 __ _____
30 \ l cc / h cc
31 H = ) M c c /------ ( b + b ),
32 el-ph /_ ij i j \/ 2 W l l
33 l,ij l
35 where the electron phonon coupling matrix is given by::
37 l ___
38 M = < i | \ / V * v |j>
39 ij 'u eff l
41 """
43 def __init__(self, atoms, indices=None, name='v', delta=0.005, nfree=2,
44 derivativemethod='tci'):
45 assert nfree in [2, 4]
46 self.nfree = nfree
47 self.delta = delta
49 if indices is None:
50 indices = range(len(self.atoms))
51 self.calc = atoms.calc
52 self.atoms = atoms
53 self.indices = np.asarray(indices)
54 self.name = name
55 self.p0 = self.atoms.positions.copy()
57 self.derivativemethod = derivativemethod
58 if derivativemethod == 'grid':
59 self.get_dP_aMix = get_grid_dP_aMix
60 elif derivativemethod == 'grid2':
61 self.get_dP_aMix = get_grid2_dP_aMix
62 elif derivativemethod == 'tci':
63 self.get_dP_aMix = get_tci_dP_aMix
64 else:
65 raise ValueError('derivativemethod must be grid, grid2, or tci')
67 def run(self):
68 if not isfile(self.name + '.eq.pckl'):
70 self.calc.calculate(self.atoms)
71 Vt_G = self.calc.get_effective_potential(pad=False)
72 Vt_G = self.calc.wfs.gd.collect(Vt_G, broadcast=True) / Ha
73 dH_asp = self.calc.hamiltonian.dH_asp
74 setups = self.calc.wfs.setups
75 nspins = self.calc.wfs.nspins
76 gd_comm = self.calc.wfs.gd.comm
78 alldH_asp = {}
79 for a, setup in enumerate(setups):
80 ni = setup.ni
81 nii = ni * (ni + 1) // 2
82 tmpdH_sp = np.zeros((nspins, nii))
83 if a in dH_asp:
84 tmpdH_sp[:] = dH_asp[a]
85 gd_comm.sum(tmpdH_sp)
86 alldH_asp[a] = tmpdH_sp
88 forces = self.atoms.get_forces()
89 self.calc.write('eq.gpw')
91 world.barrier()
92 if world.rank == 0:
93 vd = open(self.name + '.eq.pckl', 'wb')
94 fd = open('vib.eq.pckl', 'wb')
95 pickle.dump((Vt_G, alldH_asp), vd, 2)
96 pickle.dump(forces, fd)
97 vd.close()
98 fd.close()
99 world.barrier()
101 p = self.atoms.positions.copy()
102 for a in self.indices:
103 for j in range(3):
104 for sign in [-1, 1]:
105 for ndis in range(1, self.nfree / 2 + 1):
106 name = '.%d%s%s.pckl' % (a, 'xyz'[j],
107 ndis * ' +-'[sign])
108 if isfile(self.name + name):
109 continue
110 self.atoms.positions[a, j] = (p[a, j] +
111 sign * ndis * self.delta)
112 self.calc.calculate(self.atoms)
113 Vt_G = self.calc.get_effective_potential(pad=False)
114 Vt_G = self.calc.wfs.gd.collect(Vt_G,
115 broadcast=True) / Ha
116 dH_asp = self.calc.hamiltonian.dH_asp
118 alldH_asp = {}
119 for a2, setup in enumerate(setups):
120 ni = setup.ni
121 nii = ni * (ni + 1) // 2
122 tmpdH_sp = np.zeros((nspins, nii))
123 if a2 in dH_asp:
124 tmpdH_sp[:] = dH_asp[a2]
125 gd_comm.sum(tmpdH_sp)
126 alldH_asp[a2] = tmpdH_sp
128 forces = self.atoms.get_forces()
129 world.barrier()
130 if world.rank == 0:
131 vd = open(self.name + name, 'wb')
132 fd = open('vib' + name, 'wb')
133 pickle.dump((Vt_G, alldH_asp), vd)
134 pickle.dump(forces, fd)
135 vd.close()
136 fd.close()
137 world.barrier()
138 self.atoms.positions[a, j] = p[a, j]
139 self.atoms.set_positions(p)
141 def get_gradient(self):
142 """Calculates gradient"""
143 nx = len(self.indices) * 3
144 veqt_G, dHeq_asp = pickle.load(open(self.name + '.eq.pckl', 'rb'))
145 gpts = veqt_G.shape
146 dvt_Gx = np.zeros(gpts + (nx, ))
147 ddH_aspx = {}
148 for a, dH_sp in dHeq_asp.items():
149 ddH_aspx[a] = np.empty(dH_sp.shape + (nx,))
151 x = 0
152 for a in self.indices:
153 for i in range(3):
154 name = '%s.%d%s' % (self.name, a, 'xyz'[i])
155 vtm_G, dHm_asp = pickle.load(open(name + '-.pckl', 'rb'))
156 vtp_G, dHp_asp = pickle.load(open(name + '+.pckl', 'rb'))
158 if self.nfree == 4:
159 vtmm_G, dHmm_asp = pickle.load(
160 open(name + '--.pckl', 'rb'))
161 vtpp_G, dHpp_asp = pickle.load(
162 open(name + '++.pckl', 'rb'))
163 dvtdx_G = (-vtpp_G + 8.0 * vtp_G -
164 8 * vtm_G + vtmm_G) / (12 * self.delta / Bohr)
165 dvt_Gx[:, :, :, x] = dvtdx_G
166 for atom, ddH_spx in ddH_aspx.items():
167 ddH_aspx[atom][:, :, x] = (
168 -dHpp_asp[atom]
169 + 8.0 * dHp_asp[atom]
170 - 8.0 * dHm_asp[atom]
171 + dHmm_asp[atom]) / (12 * self.delta / Bohr)
172 else: # nfree = 2
173 dvtdx_G = (vtp_G - vtm_G) / (2 * self.delta / Bohr)
174 dvt_Gx[..., x] = dvtdx_G
175 for atom, ddH_spx in ddH_aspx.items():
176 ddH_aspx[atom][:, :, x] = (
177 dHp_asp[atom]
178 - dHm_asp[atom]) / (2 * self.delta / Bohr)
179 x += 1
180 return dvt_Gx, ddH_aspx
182 def get_M(self, modes, q=0):
183 r"""Calculate el-ph coupling matrix for given modes(s).
185 XXX:
186 kwarg "q=0" is an ugly hack for k-points.
187 There shuold be loops over q!
189 Note that modes must be given as a dictionary with mode
190 frequencies in eV and corresponding mode vectors in units
191 of 1/sqrt(amu), where amu = 1.6605402e-27 Kg is an atomic mass unit.
192 In short frequencies and mode vectors must be given in ase units.
194 ::
196 d d ~
197 < w | -- v | w' > = < w | -- v | w'>
198 dP dP
200 _
201 \ ~a d . ~a
202 + ) < w | p > -- /_\H < p | w' >
203 /_ i dP ij j
204 a,ij
206 _
207 \ d ~a . ~a
208 + ) < w | -- p > /_\H < p | w' >
209 /_ dP i ij j
210 a,ij
212 _
213 \ ~a . d ~a
214 + ) < w | p > /_\H < -- p | w' >
215 /_ i ij dP j
216 a,ij
218 """
220 modes1 = modes.copy()
221 # convert to atomic units
222 amu = 1.6605402e-27 # atomic unit mass [Kg]
223 me = 9.1093897e-31 # electron mass [Kg]
224 modes = {}
225 for k in modes1.keys():
226 modes[k / Ha] = modes1[k] / np.sqrt(amu / me)
228 dvt_Gx, ddH_aspx = self.get_gradient()
230 from gpaw import restart
231 atoms, calc = restart('eq.gpw', txt=None)
232 if calc.wfs.S_qMM is None:
233 calc.initialize(atoms)
234 calc.initialize_positions(atoms)
236 wfs = calc.wfs
237 nao = wfs.setups.nao
238 bfs = wfs.basis_functions
239 dtype = wfs.dtype
240 spin = 0 # XXX
242 M_lii = {}
243 parprint('Starting gradient of pseudo part')
244 for f, mode in modes.items():
245 mo = []
246 M_ii = np.zeros((nao, nao), dtype)
247 for a in self.indices:
248 mo.append(mode[a])
249 mode = np.asarray(mo).flatten()
250 dvtdP_G = np.dot(dvt_Gx, mode)
251 bfs.calculate_potential_matrix(dvtdP_G, M_ii, q=q)
252 tri2full(M_ii, 'L')
253 M_lii[f] = M_ii
254 parprint('Finished gradient of pseudo part')
256 P_aqMi = calc.wfs.P_aqMi
257 # Add the term
258 # _
259 # \ ~a d . ~a
260 # ) < w | p > -- /_\H < p | w' >
261 # /_ i dP ij j
262 # a,ij
264 Ma_lii = {}
265 for f, mode in modes.items():
266 Ma_lii[f] = np.zeros_like(M_lii.values()[0])
268 parprint('Starting gradient of dH^a part')
269 for f, mode in modes.items():
270 mo = []
271 for a in self.indices:
272 mo.append(mode[a])
273 mode = np.asarray(mo).flatten()
275 for a, ddH_spx in ddH_aspx.items():
276 ddHdP_sp = np.dot(ddH_spx, mode)
277 ddHdP_ii = unpack_hermitian(ddHdP_sp[spin])
278 Ma_lii[f] += dots(P_aqMi[a][q], ddHdP_ii, P_aqMi[a][q].T)
279 parprint('Finished gradient of dH^a part')
281 parprint('Starting gradient of projectors part')
282 dP_aMix = self.get_dP_aMix(calc.spos_ac, wfs, q)
283 parprint('Finished gradient of projectors part')
285 dH_asp = pickle.load(open('v.eq.pckl', 'rb'))[1]
287 Mb_lii = {}
288 for f, mode in modes.items():
289 Mb_lii[f] = np.zeros_like(M_lii.values()[0])
291 for f, mode in modes.items():
292 for a, dP_Mix in dP_aMix.items():
293 dPdP_Mi = np.dot(dP_Mix, mode[a])
294 dH_ii = unpack_hermitian(dH_asp[a][spin])
295 dPdP_MM = dots(dPdP_Mi, dH_ii, P_aqMi[a][q].T)
296 Mb_lii[f] -= dPdP_MM + dPdP_MM.T
297 # XXX The minus sign here is quite subtle.
298 # It is related to how the derivative of projector
299 # functions in GPAW is calculated.
300 # More thorough explanations, anyone...?
302 # Units of M_lii are Hartree/(Bohr * sqrt(m_e))
303 for mode in M_lii.keys():
304 M_lii[mode] += Ma_lii[mode] + Mb_lii[mode]
306 # conversion to eV. The prefactor 1 / sqrt(hb^2 / 2 * hb * f)
307 # has units Bohr * sqrt(me)
308 M_lii_1 = M_lii.copy()
309 M_lii = {}
311 for f in M_lii_1.keys():
312 M_lii[f * Ha] = M_lii_1[f] * Ha / np.sqrt(2 * f)
314 return M_lii
317#####################################################
318# XXX grid and grid 2 sometimes gives random numbers,
319# XXX sometimes even nan!
320#####################################################
322def get_grid_dP_aMix(spos_ac, wfs, q): # XXXXXX q
323 nao = wfs.setups.nao
324 C_MM = np.identity(nao, dtype=wfs.dtype)
325 # XXX In the future use the New Two-Center integrals
326 # to evaluate this
327 dP_aMix = {}
328 for a, setup in enumerate(wfs.setups):
329 ni = 0
330 dP_Mix = np.zeros((nao, setup.ni, 3))
331 pt = LFC(wfs.gd, [setup.pt_j],
332 wfs.kd.comm, dtype=wfs.dtype, forces=True)
333 spos1_ac = [spos_ac[a]]
334 pt.set_k_points(wfs.ibzk_qc)
335 pt.set_positions(spos1_ac)
336 for b, setup_b in enumerate(wfs.setups):
337 nao = setup_b.nao
338 phi_MG = wfs.gd.zeros(nao, wfs.dtype)
339 phi_MG = wfs.gd.collect(phi_MG, broadcast=False)
340 wfs.basis_functions.lcao_to_grid(C_MM[ni:ni + nao], phi_MG, q)
341 dP_bMix = pt.dict(len(phi_MG), derivative=True)
342 pt.derivative(phi_MG, dP_bMix, q=q)
343 dP_Mix[ni:ni + nao] = dP_bMix[0]
344 ni += nao
345 parprint(f'projector grad. doing atoms ({a}, {b}) ')
347 dP_aMix[a] = dP_Mix
348 return dP_aMix
351def get_grid2_dP_aMix(spos_ac, wfs, q, *args, **kwargs): # XXXXXX q
352 nao = wfs.setups.nao
353 C_MM = np.identity(nao, dtype=wfs.dtype)
354 bfs = wfs.basis_functions
355 phi_MG = wfs.gd.zeros(nao, wfs.dtype)
356 bfs.lcao_to_grid(C_MM, phi_MG, q)
357 setups = wfs.setups
358 pt = LFC(wfs.gd, [setup.pt_j for setup in setups],
359 wfs.kd.comm, dtype=wfs.dtype, forces=True)
360 pt.set_k_points(wfs.ibzk_qc)
361 pt.set_positions(spos_ac)
362 dP_aMix = pt.dict(len(phi_MG), derivative=True)
363 pt.derivative(phi_MG, dP_aMix, q=q)
364 return dP_aMix
367def get_tci_dP_aMix(spos_ac, wfs, q, *args, **kwargs):
368 # container for spline expansions of basis function-projector pairs
369 # (note to self: remember to conjugate/negate because of that)
370 from gpaw.lcao.overlap import ManySiteDictionaryWrapper, \
371 TwoCenterIntegralCalculator, NewTwoCenterIntegrals
373 if not isinstance(wfs.tci, NewTwoCenterIntegrals):
374 raise RuntimeError('Please remember --gpaw=usenewtci=True')
376 dP_aqxMi = {}
377 nao = wfs.setups.nao
378 nq = len(wfs.ibzk_qc)
379 for a, setup in enumerate(wfs.setups):
380 dP_aqxMi[a] = np.zeros((nq, 3, nao, setup.ni), wfs.dtype)
382 calc = TwoCenterIntegralCalculator(wfs.ibzk_qc, derivative=True)
383 expansions = ManySiteDictionaryWrapper(wfs.tci.P_expansions, dP_aqxMi)
384 calc.calculate(wfs.tci.atompairs, [expansions], [dP_aqxMi])
386 dP_aMix = {}
387 for a in dP_aqxMi:
388 dP_aMix[a] = dP_aqxMi[a].transpose(0, 2, 3, 1).copy()[q] # XXX q
389 return dP_aMix