Coverage for gpaw/dscf.py: 59%
142 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-09 00:21 +0000
1# Copyright (C) 2008 CAMd
2# Please see the accompanying LICENSE file for further information.
4"""This module is used in delta self-consistent field (dSCF) calculations
6dSCF is a simple 'ad hoc' method to estimating excitation energies within
7DFT. The only difference to ordinary DFT is that one or more electrons(s)
8are forced to occupy one or more predefined orbitals. The only restriction
9on these orbitals is that they must be linear combinations of available
10Kohn-Sham orbitals.
12"""
14import sys
15import numpy as np
17from ase.parallel import paropen
18from gpaw.utilities import devnull
19from ase.units import Ha
21from gpaw.occupations import FermiDiracCalculator
24def dscf_calculation(paw, orbitals, atoms):
25 """Helper function to prepare a calculator for a dSCF calculation
27 Parameters
28 ==========
29 orbitals: list of lists
30 Orbitals which one wants to occupy. The format is
31 orbitals = [[1.0,orb1,0],[1.0,orb2,1],...], where 1.0 is the no.
32 of electrons, orb1 and orb2 are the orbitals (see MolecularOrbitals
33 below for an example of an orbital class). 0 and 1 represents the
34 spin (up and down). This number is ignored in a spin-paired
35 calculation.
37 Example
38 =======
40 ::
42 atoms.calc = calc
43 e_gs = atoms.get_potential_energy() #ground state energy
44 sigma_star=MolecularOrbitals(calc, molecule=[0,1],
45 w=[[1.,0.,0.,0.],[-1.,0.,0.,0.]])
46 dscf_calculation(calc, [[1.0,sigma_star,1]], atoms)
47 e_exc = atoms.get_potential_energy() #excitation energy
49 """
51 # If the calculator has not been initialized the occupation object
52 # is None
53 if paw.wfs is None:
54 paw.initialize(atoms)
55 assert paw.wfs.mode != 'pw'
56 occ = paw.wfs.occupations
57 if isinstance(occ, OccupationsDSCF):
58 occ.orbitals = orbitals
59 else:
60 if occ.name == 'zero-width':
61 occ = FermiDiracCalculator(1e-4, occ.parallel_layout)
62 new_occ = OccupationsDSCF(occ, paw, orbitals)
63 paw.wfs.occupations = new_occ
64 # If the calculator has already converged (for the ground state),
65 # reset self-consistency and let the density be updated right away
66 if paw.scf.converged:
67 paw.scf.niter_fixdensity = 0
68 paw.scf.reset()
71class OccupationsDSCF:
72 """Occupation class.
74 Corresponds to the ordinary FermiDirac class in occupation.py. Only
75 difference is that it forces some electrons in the supplied orbitals
76 in stead of placing all the electrons by a Fermi-Dirac distribution.
77 """
79 def __init__(self, occ, calc, orbitals):
80 self.occ = occ
81 self.extrapolate_factor = occ.extrapolate_factor
82 self.wfs = calc.wfs
83 self.orbitals = orbitals
84 self.norbitals = len(self.orbitals)
86 self.cnoe = 0.0
87 for orb in self.orbitals:
88 self.cnoe += orb[0]
90 def calculate(self,
91 nelectrons,
92 eigenvalues,
93 weights,
94 fermi_levels_guess,
95 fix_fermi_level=False):
96 assert not fix_fermi_level
97 f_qn, fermi_levels, e_entropy = self.occ.calculate(
98 nelectrons - self.cnoe,
99 eigenvalues,
100 weights,
101 fermi_levels_guess)
103 # Get the expansion coefficients c_un for each dscf-orbital
104 # and incorporate their respective occupations into kpt.ne_o
105 c_oun = []
106 for orb in self.orbitals:
107 c_oun.append(orb[1].expand([e / Ha for e in fermi_levels],
108 self.wfs,
109 self.occ.todict().get('fixmagmom')))
111 for u, kpt in enumerate(self.wfs.kpt_u):
112 kpt.ne_o = np.zeros(self.norbitals, dtype=float)
113 kpt.c_on = np.zeros((self.norbitals, len(kpt.eps_n)),
114 dtype=complex)
116 for o, orb in enumerate(self.orbitals):
117 # TODO XXX false if orb[0]<0 since abs(c_n)**2>0
118 # kpt.c_on[o,:] = abs(orb[0])**0.5 * c_oun[o][u]
120 kpt.ne_o[o] = orb[0]
121 kpt.c_on[o, :] = c_oun[o][u]
123 if self.wfs.nspins == 2:
124 assert orb[2] in range(2), 'Invalid spin index'
126 if orb[2] == kpt.s:
127 kpt.ne_o[o] *= kpt.weight
128 else:
129 kpt.ne_o[o] = 0.0
130 else:
131 kpt.ne_o[o] *= 0.5 * kpt.weight
133 return f_qn, fermi_levels, e_entropy
135 def calculate_band_energy(self, wfs):
136 de_band = 0.0
137 for kpt in wfs.kpt_u:
138 for ne, c_n in zip(kpt.ne_o, kpt.c_on):
139 de_band += ne * np.dot(np.abs(c_n)**2, kpt.eps_n)
140 return de_band
143class MolecularOrbital:
144 r"""Class defining orbitals that should be filled in a dSCF calculation.
146 An orbital is defined through a linear combination of the atomic
147 partial waves. In each self-consistent cycle the method expand
148 is called. This method take the Kohn-Sham orbitals fulfilling the
149 criteria given by Estart, Eend and nos and return the best
150 possible expansion of the orbital in this basis. The integral
151 of the Kohn-Sham all-electron wavefunction ``|u,n>`` (u being local spin
152 and kpoint index) and the partial wave ``|\phi_i^a>`` is approximated
153 by::
155 wfs.kpt_u[u].P_ani = <\tilde p_i^a|\tilde\psi_{un}>.
157 Parameters
158 ----------
159 paw: gpaw calculator instance
160 The calculator used in the dSCF calculation.
161 Estart: float
162 Kohn-Sham orbitals with an energy above Efermi+Estart are used
163 in the linear expansion.
164 Eend: float
165 Kohn-Sham orbitals with an energy below Efermi+Eend are used
166 in the linear expansion.
167 nos: int
168 The maximum Number Of States used in the linear expansion.
169 weights: dictionary
170 The keys represent atoms and have values which are lists of weights
171 of the contributing partial waves. The default value thuis corresponds
172 to an antibonding 2\sigma orbital of atoms 0 and 1.
173 Format::
175 {1. atom: [weight of 1. projector function of the 1. atom,
176 weight of 2. projector function of the 1. atom, ...],
177 2. atom: [weight of 1. projector function of the 2. atom,
178 weight of 2. projector function of the 2. atom, ...],
179 ...}
180 """
182 def __init__(self, paw, Estart=0.0, Eend=1.e6,
183 nos=None, weights={0: [1], 1: [-1]}):
185 self.w = weights
186 self.Estart = Estart
187 self.Eend = Eend
188 self.nos = nos
190 def expand(self, epsF, wfs, fixmom):
191 if wfs.nspins == 2 and not fixmom:
192 epsF = epsF * 2
194 if self.nos is None:
195 self.nos = wfs.bd.nbands
197 c_un = []
198 for u, kpt in enumerate(wfs.kpt_u):
199 Porb_n = np.zeros(wfs.bd.nbands, dtype=complex)
200 for a, P_ni in kpt.P_ani.items():
201 if a in self.w.keys():
202 for i in range(len(self.w[a])):
203 Porb_n += self.w[a][i] * P_ni[:, i]
204 wfs.gd.comm.sum(Porb_n)
206 # Starting from KS orbitals with largest overlap,
207 # fill in the expansion coeffients between Estart and Eend
208 c_n = np.zeros(wfs.bd.nbands, dtype=complex)
209 nos = 0
210 bandpriority = np.argsort(abs(Porb_n)**2)[::-1]
212 for n in bandpriority:
213 if (kpt.eps_n[n] > epsF[kpt.s] + self.Estart and
214 kpt.eps_n[n] < epsF[kpt.s] + self.Eend):
215 c_n[n] = Porb_n[n].conj()
216 nos += 1
217 if nos == self.nos:
218 break
220 # Normalize expansion coefficients
221 c_n /= np.sqrt(sum(abs(c_n)**2))
222 c_un.append(c_n)
224 return c_un
227class AEOrbital:
228 """Class defining the orbitals that should be filled in a dSCF calculation.
230 An orbital is defined through a linear combination of KS orbitals
231 which is determined by this class as follows: For each kpoint and spin
232 we calculate the quantity ``c_n = <n|a>`` where ``|n>`` is the
233 all-electron KS states in the calculation and ``|a>`` is the
234 all-electron resonant state to be kept occupied. We can then
235 write ``|a> = Sum(c_n|n>)`` and in each self-consistent cycle the
236 method expand is called. This method take the Kohn-Sham
237 orbitals fulfilling the criteria given by Estart, Eend and
238 nos (Number Of States) and return the best possible expansion of
239 the orbital in this basis.
241 Parameters
242 ----------
243 paw: gpaw calculator instance
244 The calculator used in the dSCF calculation.
245 molecule: list of integers
246 The atoms, which are a part of the molecule.
247 Estart: float
248 Kohn-Sham orbitals with an energy above Efermi+Estart are used
249 in the linear expansion.
250 Eend: float
251 Kohn-Sham orbitals with an energy below Efermi+Eend are used
252 in the linear expansion.
253 nos: int
254 The maximum Number Of States used in the linear expansion.
255 wf_u: list of wavefunction arrays
256 Wavefunction to be occupied on the kpts on this processor:
258 wf_u = [kpt.psit_nG[n] for kpt in calc_mol.wfs.kpt_u]
260 p_uai: list of dictionaries
261 Projector overlaps with the wavefunction to be occupied for each
262 kpoint. These are used when correcting to all-electron wavefunction
263 overlaps:
265 p_uai = [dict([(mol[a], P_ni[n]) for a, P_ni in kpt.P_ani.items()])
266 for kpt in paw.wfs.kpt_u]
268 where mol is a list of atoms contributing to the state and n is the
269 band. See examples in the dscf documentation on the gpaw webpage.
270 """
272 def __init__(self, paw, wf_u, p_uai, Estart=0.0, Eend=1.e6, nos=None,
273 txt='-'):
275 self.wf_u = wf_u
276 self.p_uai = p_uai
277 self.Estart = Estart
278 self.Eend = Eend
279 self.nos = nos
281 if txt is None:
282 self.txt = devnull
283 elif txt == '-':
284 self.txt = sys.stdout
285 elif isinstance(txt, str):
286 self.txt = paropen(txt, 'w')
287 else:
288 self.txt = txt
290 def expand(self, epsF, wfs, fixmom):
292 if wfs.mode == 'lcao':
293 wfs.initialize_wave_functions_from_lcao()
294 if wfs.nspins == 2 and not fixmom:
295 epsF = epsF * 2
297 if self.nos is None:
298 self.nos = wfs.bd.nbands
300 # Check dimension of lists
301 if len(self.wf_u) == len(wfs.kpt_u):
302 wf_u = self.wf_u
303 p_uai = self.p_uai
304 else:
305 raise RuntimeError('List of wavefunctions has wrong size')
307 c_un = []
308 p_u = []
309 for u, kpt in enumerate(wfs.kpt_u):
311 # Inner product of pseudowavefunctions
312 wf = np.reshape(wf_u[u], -1)
313 Wf_n = kpt.psit_nG
314 Wf_n = np.reshape(Wf_n, (len(Wf_n), -1))
315 Porb_n = np.dot(Wf_n.conj(), wf) * wfs.gd.dv
317 # Correction to obtain inner product of AE wavefunctions
318 for a, p_i in p_uai[u].items():
319 for n in range(wfs.bd.nbands):
320 for i in range(len(p_i)):
321 for j in range(len(p_i)):
322 Porb_n[n] += (kpt.P_ani[a][n][i].conj() *
323 wfs.setups[a].dO_ii[i][j] *
324 p_i[j])
325 wfs.gd.comm.sum(Porb_n)
327 # print 'Kpt:', kpt.k, ' Spin:', kpt.s, \
328 # ' Sum_n|<orb|nks>|^2:', sum(abs(Porb_n)**2)
329 p_u.append(np.array([sum(abs(Porb_n)**2)], dtype=float))
331 # Starting from KS orbitals with largest overlap,
332 # fill in the expansion coeffients
333 c_n = np.zeros(wfs.bd.nbands, dtype=complex)
334 nos = 0
335 bandpriority = np.argsort(abs(Porb_n)**2)[::-1]
337 for n in bandpriority:
338 if (kpt.eps_n[n] > epsF[kpt.s] + self.Estart and
339 kpt.eps_n[n] < epsF[kpt.s] + self.Eend):
340 c_n[n] = Porb_n[n]
341 nos += 1
342 if nos == self.nos:
343 break
345 # Normalize expansion coefficients
346 c_n /= np.sqrt(sum(abs(c_n)**2))
347 c_un.append(c_n)
349 for s in range(wfs.nspins):
350 for k in range(wfs.kd.nibzkpts):
351 p = wfs.collect_auxiliary(p_u, k, s)
352 if wfs.world.rank == 0:
353 self.txt.write('Kpt: %d, Spin: %d, '
354 'Sum_n|<orb|nks>|^2: %f\n' % (k, s, p))
356 return c_un