Coverage for gpaw/xc/sic.py: 32%
750 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-12 00:18 +0000
1# Copyright (C) 2009-2012 P. Kluepfel and the CoMPaS group
2# Please see the accompanying LICENSE file for further information.
3#
4# ============================================================================
5# IMPORTANT INFORMATION
6# ============================================================================
7# The SIC module is developed by the CoMPaS Group, currently
8# affiliated with the Science Institute of the University of Iceland.
9#
10# The last updates of our contribution to the GPAW code took place
11# 01.11.2011 and a minor update on 16.01.2012. Since then there is no active
12# development in the original GPAW repositories, but of course we continued
13# work in our development version.
14#
15# If you are interested in running self-interaction corrected DFT
16# calculations we strongly recommend you to read the LICENSE file carefully
17# (especially regarding the FITNESS FOR A PARTICULAR PURPOSE of THIS version
18# of the code). In case of questions, a persistent interest in running SIC
19# calculations or if you want to support the necessary further development
20# of SIC and the theory of orbital-density dependent energy functionals
21# in GPAW you can reach the developers at:
22#
23# peter-Dot-kluepfel-At-gmail-Dot-com
24#
25# ============================================================================
26"""Perdew-Zunger SIC to DFT functionals.
28Self-consistent minimization of self-interaction corrected
29functionals (Perdew-Zunger).
30"""
32from math import pi
33from typing import Callable, cast
35import numpy as np
36from ase.units import Bohr, Hartree
37from scipy.linalg import eigh
39from gpaw.utilities.blas import mmmx
40from gpaw.xc import XC
41from gpaw.xc.functional import XCFunctional
42from gpaw.poisson import PoissonSolver
43from gpaw.transformers import Transformer
44from gpaw.typing import ArrayND, IntVector, RNG
45from gpaw.utilities import pack_density, unpack_hermitian
46from gpaw.lfc import LFC
47import gpaw.cgpaw as cgpaw
50def matrix_exponential(G_nn, dlt):
51 """Computes the matrix exponential of an antihermitian operator
53 U = exp(i*dlt*G)
55 G_nn: ndarray
56 anti-hermitian (skew-symmetric) matrix.
58 dlt: float
59 scaling factor for G_nn.
60 """
62 ndim = G_nn.shape[1]
63 w_n = np.zeros((ndim), dtype=float)
65 V_nn = np.zeros((ndim, ndim), dtype=complex)
66 O_nn = np.zeros((ndim, ndim), dtype=complex)
67 if G_nn.dtype == complex:
68 V_nn = 1j * G_nn.real - G_nn.imag
69 else:
70 V_nn = 1j * G_nn.real
72 w_n, V_nn = eigh(V_nn)
74 O_nn = np.diag(np.exp(1j * dlt * w_n))
76 if G_nn.dtype == complex:
77 U_nn = np.dot(V_nn, np.dot(O_nn, V_nn.T.conj())).copy()
78 else:
79 U_nn = np.dot(V_nn, np.dot(O_nn, V_nn.T.conj())).real.copy()
81 return U_nn
84def ortho(W_nn, maxerr=1E-10):
85 """ Orthogonalizes the column vectors of a matrix by Symmetric
86 Loewdin orthogonalization
88 W_nn: ndarray
89 unorthonormal matrix.
91 maxerr: float
92 maximum error for using explicit diagonalization.
94 """
96 ndim = np.shape(W_nn)[1]
98 # overlap matrix
99 O_nn = np.dot(W_nn, W_nn.T.conj())
101 # check error in orthonormality
102 err = np.sum(np.abs(O_nn - np.eye(ndim)))
103 if (err < maxerr):
104 # perturbative Symmetric-Loewdin
105 X_nn = 1.5 * np.eye(ndim) - 0.5 * O_nn
106 else:
107 # diagonalization
108 n_n = np.zeros(ndim, dtype=float)
109 n_n, U_nn = eigh(O_nn)
110 nsqrt_n = np.diag(1.0 / np.sqrt(n_n))
111 X_nn = np.dot(np.dot(U_nn, nsqrt_n), U_nn.T.conj())
113 # apply orthonormalizing transformation
114 O_nn = np.dot(X_nn, W_nn)
116 return O_nn
119def random_unitary_matrix(delta, n, rng: RNG = cast(RNG, np.random)):
120 """ Initializaes a (random) unitary matrix
122 delta: float
123 strength of deviation from unit-matrix.
124 > 0 random matrix
125 < 0 non-random matrix
127 n: int
128 dimensionality of matrix.
130 rng: numpy.random.Generator
131 optional RNG
132 """
134 assert n > 0
136 # special case n=1:
137 if n == 1:
138 return np.identity(1)
140 # non-random unitary matrix
141 if delta < 0:
142 W_nn = np.zeros((n, n))
143 for i in range(n):
144 for j in range(i):
145 W_nn[i, j] = 1.0 / (i + j)
146 W_nn = (W_nn - W_nn.T)
148 return matrix_exponential(W_nn, delta)
150 # random unitary matrix
151 elif delta > 0:
152 sample_unit_interval: Callable[[IntVector], ArrayND] = rng.random
153 W_nn = sample_unit_interval((n, n))
154 W_nn = (W_nn - W_nn.T)
156 return matrix_exponential(W_nn, delta)
157 else:
158 return np.identity(n)
161class SIC(XCFunctional):
163 orbital_dependent = True
164 unitary_invariant = False
166 def __init__(self, xc='LDA', finegrid=False, **parameters):
167 """Self-Interaction Corrected Functionals (PZ-SIC).
169 finegrid: boolean
170 Use fine grid for energy functional evaluations?
171 """
173 if isinstance(xc, str):
174 xc = XC(xc)
176 if xc.orbital_dependent:
177 raise ValueError('SIC does not support ' + xc.name)
179 self.xc = xc
180 XCFunctional.__init__(self, xc.name + '-PZ-SIC', xc.type)
181 self.finegrid = finegrid
182 self.parameters = parameters
184 # Proper IO of general SIC objects should work by means of something like:
185 def todict(self):
186 return {
187 'type': 'SIC',
188 'name': self.xc.name,
189 'finegrid': self.finegrid,
190 'parameters': dict(self.parameters)}
192 def initialize(self, density, hamiltonian, wfs):
193 assert wfs.kd.gamma
194 assert not wfs.gd.pbc_c.any()
195 assert wfs.bd.comm.size == 1 # band parallelization unsupported
197 self.wfs = wfs
198 self.dtype = float
199 self.xc.initialize(density, hamiltonian, wfs)
200 self.kpt_comm = wfs.kd.comm
201 self.nspins = wfs.nspins
202 self.nbands = wfs.bd.nbands
204 self.finegd = density.gd.refine() if self.finegrid else density.gd
205 self.ghat = LFC(self.finegd,
206 [setup.ghat_l for setup in density.setups],
207 integral=np.sqrt(4 * np.pi),
208 forces=True)
210 # XXX Use hamiltonian.poisson instead?
211 poissonsolver = PoissonSolver()
212 poissonsolver.set_grid_descriptor(self.finegd)
214 self.spin_s = {}
215 for kpt in wfs.kpt_u:
216 self.spin_s[kpt.s] = SICSpin(kpt, self.xc, density, hamiltonian,
217 wfs, poissonsolver, self.ghat,
218 self.finegd, **self.parameters)
220 def get_setup_name(self):
221 return self.xc.get_setup_name()
223 def calculate_paw_correction(self,
224 setup,
225 D_sp,
226 dEdD_sp=None,
227 addcoredensity=True,
228 a=None):
229 return self.xc.calculate_paw_correction(setup, D_sp, dEdD_sp,
230 addcoredensity, a)
232 def set_positions(self, spos_ac):
233 self.ghat.set_positions(spos_ac)
235 def calculate(self, gd, n_sg, v_sg=None, e_g=None):
236 # Normal XC contribution:
237 exc = self.xc.calculate(gd, n_sg, v_sg, e_g)
239 # SIC:
240 self.esic = 0.0
241 self.ekin = 0.0
242 for spin in self.spin_s.values():
243 if spin.kpt.psit_nG is not None:
244 desic, dekin = spin.calculate()
245 self.esic += desic
246 self.ekin += dekin
247 self.esic = self.kpt_comm.sum_scalar(self.esic)
248 self.ekin = self.kpt_comm.sum_scalar(self.ekin)
250 return exc + self.esic
252 def apply_orbital_dependent_hamiltonian(self,
253 kpt,
254 psit_nG,
255 Htpsit_nG=None,
256 dH_asp=None):
257 spin = self.spin_s[kpt.s]
258 if spin.W_mn is None:
259 return
260 spin.apply_orbital_dependent_hamiltonian(psit_nG)
262 def correct_hamiltonian_matrix(self, kpt, H_nn):
263 spin = self.spin_s[kpt.s]
264 if spin.W_mn is None:
265 return
266 spin.correct_hamiltonian_matrix(H_nn)
268 def add_correction(self,
269 kpt,
270 psit_xG,
271 Htpsit_xG,
272 P_axi,
273 c_axi,
274 n_x,
275 calculate_change=False):
276 spin = self.spin_s[kpt.s]
277 if spin.W_mn is None:
278 return
280 if calculate_change:
281 spin.calculate_residual_change(psit_xG, Htpsit_xG, P_axi, c_axi,
282 n_x)
283 else:
284 spin.calculate_residual(psit_xG, Htpsit_xG, P_axi, c_axi)
286 def rotate(self, kpt, U_nn):
287 self.spin_s[kpt.s].rotate(U_nn)
289 def setup_force_corrections(self, F_av):
290 self.dF_av = np.zeros_like(F_av)
291 for spin in self.spin_s.values():
292 spin.add_forces(self.dF_av)
293 self.wfs.kd.comm.sum(self.dF_av)
295 def add_forces(self, F_av):
296 F_av += self.dF_av
298 def summary(self, log):
299 for s in range(self.nspins):
300 if s in self.spin_s:
301 stabpot = self.spin_s[s].stabpot
302 spin = self.spin_s[s]
303 pos_mv = spin.get_centers()
304 exc_m = spin.exc_m
305 ecoulomb_m = spin.ecoulomb_m
306 if self.kpt_comm.rank == 1 and self.finegd.comm.rank == 0:
307 nocc = self.kpt_comm.sum(spin.nocc)
308 self.kpt_comm.send(pos_mv, 0)
309 self.kpt_comm.send(exc_m, 0)
310 self.kpt_comm.send(ecoulomb_m, 0)
311 else:
312 if self.kpt_comm.rank == 0 and self.finegd.comm.rank == 0:
313 nocc = self.kpt_comm.sum(0)
314 pos_mv = np.zeros((nocc, 3))
315 exc_m = np.zeros(nocc)
316 ecoulomb_m = np.zeros(nocc)
317 self.kpt_comm.receive(pos_mv, 1)
318 self.kpt_comm.receive(exc_m, 1)
319 self.kpt_comm.receive(ecoulomb_m, 1)
320 if self.kpt_comm.rank == 0 and self.finegd.comm.rank == 0:
321 log('\nSIC orbital centers and energies:')
322 log(' %5.2fx %5.2fx' %
323 (self.spin_s[0].xc_factor, self.spin_s[0].coulomb_factor))
324 log(' x y z XC Coulomb')
325 log('--------------------------------------------------')
326 m = 0
327 for pos_v, exc, ecoulomb in zip(pos_mv, exc_m, ecoulomb_m):
328 log('%3d (%7.3f,%7.3f,%7.3f): %8.3f %8.3f' %
329 ((m, ) + tuple(pos_v) +
330 (exc * Hartree, ecoulomb * Hartree)))
331 m += 1
332 log('--------------------------------------------------')
333 log('\nTotal SIC energy : %12.5f' % (self.esic * Hartree))
334 log('Stabilizing potential: %12.5f' % (stabpot * Hartree))
336 def read(self, reader):
337 xc_factor = reader.hamiltonian.xc.sic_xc_factor
338 coulomb_factor = reader.hamiltonian.xc.sic_coulomb_factor
340 for s in range(self.nspins):
341 W_mn = reader.hamiltonian.xc.get(
342 f'unitary_transformation{s}')
344 if s in self.spin_s:
345 self.spin_s[s].initial_W_mn = W_mn
346 self.spin_s[s].xc_factor = xc_factor
347 self.spin_s[s].coulomb_factor = coulomb_factor
349 def write(self, writer):
350 for s in self.spin_s:
351 spin = self.spin_s[s]
352 writer.write(
353 sic_xc_factor=spin.xc_factor,
354 sic_coulomb_factor=spin.coulomb_factor)
355 break
357 for s in range(self.nspins):
358 W_mn = self.get_unitary_transformation(s)
360 if W_mn is not None:
361 writer.write(f'unitary_transformation{s}', W_mn)
363 def get_unitary_transformation(self, s):
364 if s in self.spin_s.keys():
365 spin = self.spin_s[s]
367 if spin.W_mn is None or spin.finegd.rank != 0:
368 n = 0
369 else:
370 n = spin.W_mn.shape[0]
371 else:
372 n = 0
374 n = self.wfs.world.sum_scalar(n)
376 if n > 0:
377 W_mn = np.zeros((n, n), dtype=self.dtype)
378 else:
379 W_mn = None
380 return W_mn
382 if s in self.spin_s.keys():
383 spin = self.spin_s[s]
384 #
385 if spin.W_mn is None or spin.finegd.rank != 0:
386 W_mn[:] = 0.0
387 else:
388 W_mn[:] = spin.W_mn[:]
389 #
390 else:
391 W_mn[:] = 0.0
393 self.wfs.world.sum(W_mn)
394 return W_mn
397class SICSpin:
398 def __init__(self,
399 kpt,
400 xc,
401 density,
402 hamiltonian,
403 wfs,
404 poissonsolver,
405 ghat,
406 finegd,
407 dtype=float,
408 coulomb_factor=0.5,
409 xc_factor=0.5,
410 uominres=1E-1,
411 uomaxres=1E-10,
412 uorelres=1E-4,
413 uonscres=1E-10,
414 rattle=-0.1,
415 stabpot=0.0,
416 maxuoiter=10,
417 logging=2,
418 rng: RNG = cast(RNG, np.random)):
419 """Single spin SIC object.
422 coulomb_factor:
423 Scaling factor for Hartree-functional
425 xc_factor:
426 Scaling factor for xc-functional
428 uominres:
429 Minimum residual before unitary optimization starts
431 uomaxres:
432 Target accuracy for unitary optimization
433 (absolute variance)
435 uorelres:
436 Target accuracy for unitary optimization
437 (rel. to basis residual)
439 maxuoiter:
440 Maximum number of unitary optimization steps
442 rattle:
443 perturbation to the initial states
444 rng:
445 optional numpy.random.Generator instance
446 """
448 self.wfs = wfs
449 self.kpt = kpt
450 self.xc = xc
451 self.poissonsolver = poissonsolver
452 self.ghat = ghat
453 self.pt = wfs.pt
455 self.gd = wfs.gd
456 self.finegd = finegd
458 if self.finegd is self.gd:
459 self.interpolator = None
460 self.restrictor = None
461 else:
462 self.interpolator = Transformer(self.gd, self.finegd, 3)
463 self.restrictor = Transformer(self.finegd, self.gd, 3)
465 self.nspins = wfs.nspins
466 self.spin = kpt.s
467 self.timer = wfs.timer
468 self.setups = wfs.setups
470 self.dtype = dtype
471 self.coulomb_factor = coulomb_factor
472 self.xc_factor = xc_factor
474 self.nocc = None # number of occupied states
475 self.W_mn = None # unit. transf. to energy optimal states
476 self.initial_W_mn = None # initial unitary transformation
477 self.vt_mG = None # orbital dependent potentials
478 self.exc_m = None # PZ-SIC contributions (from E_xc)
479 self.ecoulomb_m = None # PZ-SIC contributions (from E_H)
481 self.rattle = rattle # perturb the initial unitary transformation
482 self.stabpot = stabpot # stabilizing constant potential to avoid
483 # occupation of unoccupied states
485 self.uominres = uominres
486 self.uomaxres = uomaxres
487 self.uorelres = uorelres
488 self.uonscres = uonscres
489 self.maxuoiter = maxuoiter
490 self.maxlsiter = 1 # maximum number of line-search steps
491 self.maxcgiter = 2 # maximum number of CG-iterations
492 self.lsinterp = True # interpolate for minimum during line search
493 self.basiserror = 1E+20
494 self.logging = logging
496 self.rng = rng
498 def initialize_orbitals(self, rattle=-0.1, localize=True):
499 if self.initial_W_mn is not None:
500 self.nocc = self.initial_W_mn.shape[0]
501 else:
502 self.nocc, x = divmod(int(self.kpt.f_n.sum()), 3 - self.nspins)
503 assert x == 0
505 if self.nocc == 0:
506 return
508 Z_mmv = np.empty((self.nocc, self.nocc, 3), complex)
509 for v in range(3):
510 G_v = np.zeros(3)
511 G_v[v] = 1
512 Z_mmv[:, :, v] = self.gd.wannier_matrix(
513 self.kpt.psit_nG, self.kpt.psit_nG, G_v, self.nocc)
514 self.finegd.comm.sum(Z_mmv)
516 if self.initial_W_mn is not None:
517 self.W_mn = self.initial_W_mn
519 elif localize:
520 W_nm = np.identity(self.nocc)
521 localization = 0.0
522 for iter in range(30):
523 loc = cgpaw.localize(Z_mmv, W_nm)
524 if loc - localization < 1e-6:
525 break
526 localization = loc
528 self.W_mn = W_nm.T.copy()
529 else:
530 self.W_mn = np.identity(self.nocc)
532 if (rattle != 0.0 and self.W_mn is not None and
533 self.initial_W_mn is None):
534 U_mm = random_unitary_matrix(rattle, self.nocc, rng=self.rng)
535 self.W_mn = np.dot(U_mm, self.W_mn)
537 if self.W_mn is not None:
538 self.finegd.comm.broadcast(self.W_mn, 0)
540 spos_mc = -np.angle(Z_mmv.diagonal()).T / (2 * pi)
541 self.initial_pos_mv = np.dot(spos_mc % 1.0, self.gd.cell_cv)
543 def localize_orbitals(self):
545 assert self.gd.orthogonal
547 # calculate wannier matrixelements
548 Z_mmv = np.empty((self.nocc, self.nocc, 3), complex)
549 for v in range(3):
550 G_v = np.zeros(3)
551 G_v[v] = 1
552 Z_mmv[:, :, v] = self.gd.wannier_matrix(
553 self.kpt.psit_nG, self.kpt.psit_nG, G_v, self.nocc)
554 self.finegd.comm.sum(Z_mmv)
556 # setup the initial configuration (identity)
557 W_nm = np.identity(self.nocc)
559 # localize the orbitals
560 localization = 0.0
561 for iter in range(30):
562 loc = cgpaw.localize(Z_mmv, W_nm)
563 if loc - localization < 1e-6:
564 break
565 localization = loc
567 # apply localizing transformation
568 self.W_mn = W_nm.T.copy()
570 def rattle_orbitals(self, rattle=-0.1):
572 # check for the trivial cases
573 if rattle == 0.0:
574 return
576 if self.W_mn is None:
577 return
579 # setup a "random" unitary matrix
580 nocc = self.W_mn.shape[0]
581 U_mm = random_unitary_matrix(rattle, nocc, rng=self.rng)
583 # apply unitary transformation
584 self.W_mn = np.dot(U_mm, self.W_mn)
586 def get_centers(self):
587 assert self.gd.orthogonal
589 # calculate energy optimal states (if necessary)
590 if self.phit_mG is None:
591 self.update_optimal_states()
593 # calculate wannier matrixelements
594 Z_mmv = np.empty((self.nocc, self.nocc, 3), complex)
595 for v in range(3):
596 G_v = np.zeros(3)
597 G_v[v] = 1
598 Z_mmv[:, :, v] = self.gd.wannier_matrix(self.phit_mG, self.phit_mG,
599 G_v, self.nocc)
600 self.finegd.comm.sum(Z_mmv)
602 # calculate positions of localized orbitals
603 spos_mc = -np.angle(Z_mmv.diagonal()).T / (2 * pi)
605 return np.dot(spos_mc % 1.0, self.gd.cell_cv) * Bohr
607 def calculate_sic_matrixelements(self):
608 # overlap of pseudo wavefunctions
609 Htphit_mG = self.vt_mG * self.phit_mG
610 V_mm = np.zeros((self.nocc, self.nocc), dtype=self.dtype)
611 # gemm(self.gd.dv, self.phit_mG, Htphit_mG, 0.0, V_mm, 't')
612 mmmx(self.gd.dv,
613 Htphit_mG, 'N',
614 self.phit_mG, 'T', 0.0, V_mm)
616 # PAW
617 for a, P_mi in self.P_ami.items():
618 for m, dH_p in enumerate(self.dH_amp[a]):
619 dH_ii = unpack_hermitian(dH_p)
620 V_mm[m, :] += np.dot(P_mi[m], np.dot(dH_ii, P_mi.T))
622 # accumulate over grid-domains
623 self.finegd.comm.sum(V_mm)
624 self.V_mm = V_mm
626 # Symmetrization of V and kappa-matrix:
627 K_mm = 0.5 * (V_mm - V_mm.T.conj())
628 V_mm = 0.5 * (V_mm + V_mm.T.conj())
630 # evaluate the kinetic correction
631 self.ekin = -np.trace(V_mm) * (3 - self.nspins)
633 return V_mm, K_mm, np.vdot(K_mm, K_mm).real
635 def update_optimal_states(self):
636 # pseudo wavefunctions
637 self.phit_mG = self.gd.zeros(self.nocc)
638 if self.nocc > 0:
639 mmmx(1.0, self.W_mn, 'N', self.kpt.psit_nG[:self.nocc], 'N', 0.0,
640 self.phit_mG)
642 # PAW
643 self.P_ami = {}
644 for a, P_ni in self.kpt.P_ani.items():
645 self.P_ami[a] = np.dot(self.W_mn, P_ni[:self.nocc])
647 def calculate_density(self, m):
648 # pseudo density
649 nt_G = self.phit_mG[m]**2
651 if self.finegd is self.gd:
652 nt_g = nt_G
653 else:
654 nt_g = self.finegd.empty()
655 self.interpolator.apply(nt_G, nt_g)
656 # normalize single-particle density
657 nt_g *= self.gd.integrate(nt_G) / self.finegd.integrate(nt_g)
659 # PAW corrections
660 Q_aL = {}
661 D_ap = {}
662 for a, P_mi in self.P_ami.items():
663 P_i = P_mi[m]
664 D_ii = np.outer(P_i, P_i.conj()).real
665 D_ap[a] = D_p = pack_density(D_ii)
666 Q_aL[a] = np.dot(D_p, self.setups[a].Delta_pL)
668 return nt_g, Q_aL, D_ap
670 def update_potentials(self, save=False, restore=False):
671 if restore:
672 self.exc_m = self.exc_save_m
673 self.ecoulomb_m = self.eco_save_m
674 self.esic = self.esic_save
675 self.vt_mG = self.vt_save_mG.copy()
676 self.dH_amp = self.dH_save_amp.copy()
677 return
679 self.timer.start('ODD-potentials')
681 nt_sg = self.finegd.empty(2)
682 nt_sg[1] = 0.0
683 vt_sg = self.finegd.empty(2)
685 # PAW
686 W_aL = self.ghat.dict()
687 zero_initial_phi = False
689 # initialize some bigger fields
690 if self.vt_mG is None or self.nocc != self.phit_mG.shape[0]:
691 self.vt_mG = self.gd.empty(self.nocc)
692 self.exc_m = np.zeros(self.nocc)
693 self.ecoulomb_m = np.zeros(self.nocc)
694 self.vHt_mg = self.finegd.zeros(self.nocc)
695 zero_initial_phi = True
696 #
697 # PAW
698 self.dH_amp = {}
699 for a, P_mi in self.P_ami.items():
700 ni = P_mi.shape[1]
701 self.dH_amp[a] = np.empty((self.nocc, ni * (ni + 1) // 2))
702 #
703 self.Q_maL = {}
704 # loop all occupied orbitals
705 for m, phit_G in enumerate(self.phit_mG):
706 #
707 # setup the single-particle density and PAW density-matrix
708 nt_sg[0], Q_aL, D_ap = self.calculate_density(m)
709 vt_sg[:] = 0.0
710 #
711 # xc-SIC
712 self.timer.start('XC')
713 if self.xc_factor != 0.0:
714 exc = self.xc.calculate(self.finegd, nt_sg, vt_sg)
715 exc /= self.finegd.comm.size
716 vt_sg[0] *= -self.xc_factor
717 #
718 # PAW
719 for a, D_p in D_ap.items():
720 setup = self.setups[a]
721 dH_p = self.dH_amp[a][m]
722 dH_sp = np.zeros((2, len(dH_p)))
723 #
724 D_sp = np.array([D_p, np.zeros_like(D_p)])
725 exc += self.xc.calculate_paw_correction(
726 setup, D_sp, dH_sp, addcoredensity=False)
727 dH_p[:] = -dH_sp[0] * self.xc_factor
729 self.exc_m[m] = -self.xc_factor * exc
730 self.timer.stop('XC')
731 #
732 # Hartree-SIC
733 self.timer.start('Hartree')
734 if self.coulomb_factor != 0.0:
735 #
736 # add compensation charges to pseudo density
737 self.ghat.add(nt_sg[0], Q_aL)
738 #
739 # solve the coulomb problem
740 self.poissonsolver.solve(
741 self.vHt_mg[m],
742 nt_sg[0],
743 zero_initial_phi=zero_initial_phi)
744 ecoulomb = 0.5 * self.finegd.integrate(nt_sg[0] *
745 self.vHt_mg[m])
746 ecoulomb /= self.finegd.comm.size
747 vt_sg[0] -= self.coulomb_factor * self.vHt_mg[m]
748 #
749 # PAW
750 self.ghat.integrate(self.vHt_mg[m], W_aL)
751 for a, D_p in D_ap.items():
752 setup = self.setups[a]
753 dH_p = self.dH_amp[a][m]
754 M_p = np.dot(setup.M_pp, D_p)
755 ecoulomb += np.dot(D_p, M_p)
756 #
757 dH_p -= self.coulomb_factor * (
758 2.0 * M_p + np.dot(setup.Delta_pL, W_aL[a]))
759 #
760 self.ecoulomb_m[m] = -self.coulomb_factor * ecoulomb
761 self.timer.stop('Hartree')
763 # restrict to course grid
764 if self.finegd is self.gd:
765 self.vt_mG[m] = vt_sg[0]
766 else:
767 self.timer.start('Restrictor')
768 self.restrictor.apply(vt_sg[0], self.vt_mG[m])
769 self.timer.stop('Restrictor')
770 self.Q_maL[m] = Q_aL
772 self.timer.stop('ODD-potentials')
774 # accumulate total xc-SIC and coulomb-SIC
775 self.finegd.comm.sum(self.exc_m)
776 self.finegd.comm.sum(self.ecoulomb_m)
778 # total sic (including spin-degeneracy)
779 self.esic = self.exc_m.sum()
780 self.esic += self.ecoulomb_m.sum()
781 self.esic *= (3 - self.nspins)
783 if save:
784 self.exc_save_m = self.exc_m
785 self.eco_save_m = self.ecoulomb_m
786 self.esic_save = self.esic
787 self.vt_save_mG = self.vt_mG.copy()
788 self.dH_save_amp = self.dH_amp.copy()
790 def apply_orbital_dependent_hamiltonian(self, psit_nG):
791 """...
793 Setup::
795 |V phi_m> and <l|Vphi_m>,
797 for occupied states m and unoccupied states l."""
799 # nocc = self.nocc
800 # nvirt = psit_nG.shape[0] - nocc
802 self.Htphit_mG = self.vt_mG * self.phit_mG
804 def correct_hamiltonian_matrix(self, H_nn):
805 """ Add contributions of the non-local part of the
806 interaction potential to the Hamiltonian matrix.
808 on entry:
809 H_nn[n,m] = <n|H_{KS}|m>
810 on exit:
811 H_nn[n,m] = <n|H_{KS} + V_{u}|m>
813 where V_u is the unified Hamiltonian
815 V_u = ...
817 """
819 nocc = self.nocc
820 nvirt = H_nn.shape[0] - nocc
822 W_mn = self.W_mn
823 # R_mk = self.R_mk
825 if self.gd.comm.rank == 0:
826 V_mm = 0.5 * (self.V_mm + self.V_mm.T)
827 H_nn[:nocc, :nocc] += np.dot(W_mn.T, np.dot(V_mm, W_mn))
828 if self.stabpot != 0.0:
829 H_nn[nocc:, nocc:] += np.eye(nvirt) * self.stabpot
831 if nvirt != 0:
832 H_nn[:nocc, nocc:] = 0.0 # R_nk
833 H_nn[nocc:, :nocc] = 0.0 # R_nk.T
834 # R_nk = np.dot(W_mn.T, R_mk) # CHECK THIS
835 # H_nn[:nocc, nocc:] += R_nk
836 # H_nn[nocc:, :nocc] += R_nk.T
838 def calculate_residual(self, psit_nG, Htpsit_nG, P_ani, c_ani):
839 """ Calculate the action of the unified Hamiltonian on an
840 arbitrary state:
842 H_u|Psi> =
843 """
845 nocc = self.nocc
846 nvirt = psit_nG.shape[0] - nocc
848 # constraint for unoccupied states
849 R_mk = np.zeros((nocc, nvirt), dtype=self.dtype)
850 if nvirt > 0:
851 mmmx(self.gd.dv, self.Htphit_mG, 'N', psit_nG[nocc:], 'T',
852 0.0, R_mk)
853 # PAW
854 for a, P_mi in self.P_ami.items():
855 P_ni = P_ani[a]
857 for m, dH_p in enumerate(self.dH_amp[a]):
858 dH_ii = unpack_hermitian(dH_p)
859 R_mk[m] += np.dot(P_mi[m], np.dot(dH_ii, P_ni[nocc:].T))
861 self.finegd.comm.sum(R_mk)
863 # self.R_mk = R_mk
864 # R_mk = self.R_mk
865 W_mn = self.W_mn
866 Htphit_mG = self.Htphit_mG
867 phit_mG = self.phit_mG
868 K_mm = 0.5 * (self.V_mm - self.V_mm.T)
869 Q_mn = np.dot(K_mm, W_mn)
871 # Action of unified Hamiltonian on occupied states:
872 if nocc > 0:
873 mmmx(1.0, W_mn.T.copy(), 'N', Htphit_mG, 'N',
874 1.0, Htpsit_nG[:nocc])
875 mmmx(1.0, Q_mn.T.copy(), 'N', phit_mG, 'N', 1.0, Htpsit_nG[:nocc])
876 if nvirt > 0:
877 mmmx(1.0, R_mk.T.copy(), 'N', phit_mG, 'N', 1.0, Htpsit_nG[nocc:])
878 if self.stabpot != 0.0:
879 Htpsit_nG[nocc:] += self.stabpot * psit_nG[nocc:]
881 # PAW
882 for a, P_mi in self.P_ami.items():
883 #
884 c_ni = c_ani[a]
885 ct_mi = P_mi.copy()
886 #
887 dO_ii = self.setups[a].dO_ii
888 c_ni[:nocc] += np.dot(Q_mn.T, np.dot(P_mi, dO_ii))
889 c_ni[nocc:] += np.dot(R_mk.T, np.dot(P_mi, dO_ii))
890 #
891 for m, dH_p in enumerate(self.dH_amp[a]):
892 dH_ii = unpack_hermitian(dH_p)
893 ct_mi[m] = np.dot(P_mi[m], dH_ii)
894 c_ni[:nocc] += np.dot(W_mn.T, ct_mi)
895 c_ni[nocc:] += self.stabpot * np.dot(P_ani[a][nocc:], dO_ii)
897 def calculate_residual_change(self, psit_xG, Htpsit_xG, P_axi, c_axi, n_x):
898 """
900 """
901 assert len(n_x) == 1
903 Htphit_mG = self.Htphit_mG
904 phit_mG = self.phit_mG
906 w_mx = np.zeros((self.nocc, 1), dtype=self.dtype)
907 v_mx = np.zeros((self.nocc, 1), dtype=self.dtype)
909 mmmx(self.gd.dv, phit_mG, 'N', psit_xG, 'T', 0.0, w_mx)
910 mmmx(self.gd.dv, Htphit_mG, 'N', psit_xG, 'T', 0.0, v_mx)
912 # PAW
913 for a, P_mi in self.P_ami.items():
914 P_xi = P_axi[a]
915 dO_ii = self.setups[a].dO_ii
916 #
917 w_mx += np.dot(P_mi, np.dot(dO_ii, P_xi.T))
919 for m, dH_p in enumerate(self.dH_amp[a]):
920 dH_ii = unpack_hermitian(dH_p)
921 v_mx[m] += np.dot(P_mi[m], np.dot(dH_ii, P_xi.T))
923 # sum over grid-domains
924 self.finegd.comm.sum(w_mx)
925 self.finegd.comm.sum(v_mx)
927 V_mm = 0.5 * (self.V_mm + self.V_mm.T)
928 q_mx = v_mx - np.dot(V_mm, w_mx)
930 if self.stabpot != 0.0:
931 q_mx -= self.stabpot * w_mx
933 mmmx(1.0, w_mx.T.copy(), 'N', Htphit_mG, 'N', 1.0, Htpsit_xG)
934 mmmx(1.0, q_mx.T.copy(), 'N', phit_mG, 'N', 1.0, Htpsit_xG)
936 # PAW
937 for a, P_mi in self.P_ami.items():
938 c_xi = c_axi[a]
939 ct_mi = P_mi.copy()
941 dO_ii = self.setups[a].dO_ii
942 c_xi += np.dot(q_mx.T, np.dot(P_mi, dO_ii))
944 for m, dH_p in enumerate(self.dH_amp[a]):
945 dH_ii = unpack_hermitian(dH_p)
946 ct_mi[m] = np.dot(P_mi[m], dH_ii)
947 c_xi += np.dot(w_mx.T, ct_mi)
949 if self.stabpot != 0.0:
950 Htphit_mG += self.stabpot * psit_xG
951 for a, P_xi in P_axi.items():
952 dO_ii = self.setups[a].dO_ii
953 c_axi[a] += self.stabpot * np.dot(P_xi, dO_ii)
955 def rotate(self, U_nn):
956 """ Unitary transformations amongst the canonic states
957 require to apply a counter-acting transformation to
958 the energy optimal states. This subroutine takes
959 care of it.
961 Reorthogonalization is required whenever unoccupied
962 states are mixed in.
963 """
964 # skip if no transformation to optimal states is set-up
965 if self.W_mn is None:
966 return
968 # compensate the transformation amongst the occupied states
969 self.W_mn = np.dot(self.W_mn, U_nn[:self.nocc, :self.nocc])
971 # reorthogonalize if unoccupied states may have been mixed in
972 if self.nocc != U_nn.shape[0]:
973 self.W_mn = ortho(self.W_mn)
974 # self.R_mk = np.dot(self.R_mk, U_nn[self.nocc:, self.nocc:].T)
976 def add_forces(self, F_av):
977 # Calculate changes in projections
978 deg = 3 - self.nspins
979 F_amiv = self.pt.dict(self.nocc, derivative=True)
980 self.pt.derivative(self.phit_mG, F_amiv)
981 for m in range(self.nocc):
982 # Force from compensation charges:
983 dF_aLv = self.ghat.dict(derivative=True)
984 self.ghat.derivative(self.vHt_mg[m], dF_aLv)
985 for a, dF_Lv in dF_aLv.items():
986 F_av[a] -= deg * self.coulomb_factor * \
987 np.dot(self.Q_maL[m][a], dF_Lv)
989 # Force from projectors
990 for a, F_miv in F_amiv.items():
991 F_vi = F_miv[m].T.conj()
992 dH_ii = unpack_hermitian(self.dH_amp[a][m])
993 P_i = self.P_ami[a][m]
994 F_v = np.dot(np.dot(F_vi, dH_ii), P_i)
995 F_av[a] += deg * 2 * F_v.real
997 def calculate(self):
998 """ Evaluate the non-unitary invariant part of the
999 energy functional and returns
1001 esic: float
1002 self-interaction energy
1004 ekin: float
1005 correction to the kinetic energy
1006 """
1007 # initialize transformation from canonic to energy
1008 # optimal states (if necessary)
1009 if self.W_mn is None:
1010 self.initialize_orbitals(rattle=self.rattle)
1012 # optimize the non-unitary invariant part of the
1013 # functional
1014 self.unitary_optimization()
1016 return self.esic, self.ekin
1018 def unitary_optimization(self):
1019 """ Optimization of the non-unitary invariant part of the
1020 energy functional.
1021 """
1023 optstep = 0.0
1024 Gold = 0.0
1025 cgiter = 0
1026 #
1027 epsstep = 0.005 # 0.005
1028 dltstep = 0.1 # 0.1
1029 prec = 1E-7
1031 # get the initial ODD potentials/energy/matrixelements
1032 self.update_optimal_states()
1033 self.update_potentials(save=True)
1034 ESI = self.esic
1035 V_mm, K_mm, norm = self.calculate_sic_matrixelements()
1037 if norm < self.uonscres and self.maxuoiter > 0:
1038 return
1040 if self.nocc <= 1:
1041 return
1043 # optimize the unitary transformation
1044 #
1045 # allocate arrays for the search direction,
1046 # i.e., the (conjugate) gradient
1047 D_old_mm = np.zeros_like(self.W_mn)
1049 for iter in range(abs(self.maxuoiter)):
1050 # copy the initial unitary transformation and orbital
1051 # dependent energies
1052 W_old_mn = self.W_mn.copy()
1054 # setup the steepest-descent/conjugate gradient
1055 # D_nn: search direction
1056 # K_nn: inverse gradient
1057 # G0 : <K,D> (projected length of D along K)
1058 if Gold != 0.0:
1059 # conjugate gradient
1060 G0 = np.sum(K_mm * K_mm.conj()).real
1061 beta = G0 / Gold
1062 Gold = G0
1063 D_mm = K_mm + beta * D_old_mm
1065 G0 = np.sum(K_mm * D_mm.conj()).real
1066 else:
1067 # steepest-descent
1068 G0 = np.sum(K_mm * K_mm.conj()).real
1069 Gold = G0
1070 D_mm = K_mm
1072 updated = False
1073 minimum = False
1074 failed = True
1075 noise = False
1076 E0 = ESI
1078 # try to estimate optimal step-length from change in length
1079 # of the gradient (force-only)
1080 if (epsstep != 0.0 and norm > self.uomaxres):
1081 step = epsstep
1082 while (True):
1083 U_mm = matrix_exponential(D_mm, step)
1084 self.W_mn = np.dot(U_mm, W_old_mn)
1085 self.update_optimal_states()
1086 self.update_potentials()
1087 E1 = self.esic
1088 K0_mm = K_mm.copy()
1089 V_mm, K_mm, norm = self.calculate_sic_matrixelements()
1091 # projected length of the gradient at the new position
1092 G1 = np.sum(((K_mm - K0_mm) / step) * D_mm.conj()).real
1093 #
1094 if (abs(E1 - E0) < prec and E1 >= E0):
1095 #
1096 eps_works = True
1097 Eeps = E1
1098 noise = True
1099 elif (E1 < E0):
1100 #
1101 # trial step reduced energy
1102 eps_works = True
1103 Eeps = E1
1104 else:
1105 #
1106 # epsilon step did not work
1107 eps_works = False
1108 optstep = 0.0
1109 break
1111 # compute the optimal step size
1112 # optstep = step*G0/(G0-G1)
1113 # print -G0/G1
1114 optstep = -G0 / G1
1116 if (eps_works):
1117 break
1119 # print eps_works, optstep, G0/((G0-G1)/step)
1121 # decide on the method for stepping
1122 if (optstep > 0.0):
1124 # convex region -> force only estimate for minimum
1125 U_mm = matrix_exponential(D_mm, optstep)
1126 self.W_mn = np.dot(U_mm, W_old_mn)
1127 self.update_optimal_states()
1128 self.update_potentials()
1129 E1 = self.esic
1130 if (abs(E1 - E0) < prec and E1 >= E0):
1131 V_mm, K_mm, norm = self.calculate_sic_matrixelements()
1132 ESI = E1
1133 optstep = optstep
1134 lsiter = 0
1135 maxlsiter = -1
1136 updated = True
1137 minimum = True
1138 failed = False
1139 lsmethod = 'CV-N'
1140 noise = True
1141 elif (E1 < E0):
1142 V_mm, K_mm, norm = self.calculate_sic_matrixelements()
1143 ESI = E1
1144 optstep = optstep
1145 lsiter = 0
1146 maxlsiter = -1
1147 updated = True
1148 minimum = True
1149 failed = False
1150 lsmethod = 'CV-S'
1151 else:
1152 # self.K_unn[q] = K_nn
1153 ESI = E0
1154 step = optstep
1155 optstep = 0.0
1156 lsiter = 0
1157 maxlsiter = self.maxlsiter
1158 updated = False
1159 minimum = False
1160 failed = True
1161 lsmethod = 'CV-F-CC'
1162 else:
1163 # self.K_unn[q] = K_nn
1164 ESI = E0
1165 step = optstep
1166 optstep = 0.0
1167 lsiter = 0
1168 maxlsiter = self.maxlsiter
1169 updated = False
1170 minimum = False
1171 failed = True
1172 lsmethod = 'CC'
1173 else:
1174 maxlsiter = 0
1175 lsiter = -1
1176 optstep = epsstep
1177 updated = False
1178 minimum = True
1179 failed = False
1180 lsmethod = 'CC-EPS'
1182 if (optstep == 0.0):
1183 #
1184 # we are in the concave region or force-only estimate failed,
1185 # just follow the (conjugate) gradient
1186 step = dltstep * abs(step)
1187 U_mm = matrix_exponential(D_mm, step)
1188 self.W_mn = np.dot(U_mm, W_old_mn)
1189 self.update_optimal_states()
1190 self.update_potentials()
1191 E1 = self.esic
1192 #
1193 #
1194 if (abs(E1 - E0) < prec and E1 >= E0):
1195 ESI = E1
1196 optstep = 0.0
1197 updated = False
1198 minimum = True
1199 failed = True
1200 lsmethod = lsmethod + '-DLT-N'
1201 noise = True
1202 maxlsiter = -1
1203 elif (E1 < E0):
1204 ESI = E1
1205 optstep = step
1206 updated = True
1207 minimum = False
1208 failed = False
1209 lsmethod = lsmethod + '-DLT'
1210 maxlsiter = self.maxlsiter
1211 elif (eps_works):
1212 ESI = Eeps
1213 E1 = Eeps
1214 step = epsstep
1215 updated = False
1216 minimum = False
1217 failed = False
1218 lsmethod = lsmethod + '-EPS'
1219 maxlsiter = self.maxlsiter
1220 else:
1221 optstep = 0.0
1222 updated = False
1223 minimum = False
1224 failed = True
1225 lsmethod = lsmethod + '-EPS-F'
1226 maxlsiter = -1
1227 #
1228 G = (E1 - E0) / step
1229 step0 = 0.0
1230 step1 = step
1231 step2 = 2 * step
1232 #
1233 for lsiter in range(maxlsiter):
1234 #
1235 # energy at the new position
1236 U_mm = matrix_exponential(D_mm, step2)
1237 self.W_mn = np.dot(U_mm, W_old_mn)
1238 self.update_optimal_states()
1239 self.update_potentials()
1240 E2 = self.esic
1241 G = (E2 - E1) / (step2 - step1)
1242 #
1243 #
1244 if (G > 0.0):
1245 if self.lsinterp:
1246 a = (E0 / ((step2 - step0) * (step1 - step0)) +
1247 E2 / ((step2 - step1) * (step2 - step0)) -
1248 E1 / ((step2 - step1) * (step1 - step0)))
1249 b = (E2 - E0) / (step2 - step0) - a * (step2 +
1250 step0)
1251 if (a < step**2):
1252 optstep = 0.5 * (step0 + step2)
1253 else:
1254 optstep = -0.5 * b / a
1255 updated = False
1256 minimum = True
1257 break
1258 else:
1259 optstep = step1
1260 updated = False
1261 minimum = True
1262 break
1264 elif (G < 0.0):
1265 optstep = step2
1266 step0 = step1
1267 step1 = step2
1268 step2 = step2 + step
1269 E0 = E1
1270 E1 = E2
1271 ESI = E2
1272 updated = True
1273 minimum = False
1275 if (cgiter != 0):
1276 lsmethod = lsmethod + ' CG'
1278 if (cgiter >= self.maxcgiter or not minimum):
1279 Gold = 0.0
1280 cgiter = 0
1281 else:
1282 cgiter = cgiter + 1
1283 D_old_mm[:, :] = D_mm[:, :]
1285 # update the energy and matrixelements of V and Kappa
1286 # and accumulate total residual of unitary optimization
1287 if (not updated):
1288 if optstep > 0.0:
1289 U_mm = matrix_exponential(D_mm, optstep)
1290 self.W_mn = np.dot(U_mm, W_old_mn)
1291 self.update_optimal_states()
1292 self.update_potentials()
1293 ESI = self.esic
1294 V_mm, K_mm, norm = self.calculate_sic_matrixelements()
1295 else:
1296 self.W_mn = W_old_mn
1297 self.update_optimal_states()
1298 self.update_potentials(restore=True)
1299 ESI = self.esic
1300 V_mm, K_mm, norm = self.calculate_sic_matrixelements()
1302 if (lsiter == maxlsiter - 1):
1303 V_mm, K_mm, norm = self.calculate_sic_matrixelements()
1305 E0 = ESI
1307 # orthonormalize the energy optimal orbitals
1308 self.W_mn = ortho(self.W_mn)
1309 K = max(norm, 1.0e-16)
1311 if self.finegd.comm.rank == 0:
1312 if self.logging == 1:
1313 print(" UO-%d: %2d %5.1f %20.5f " %
1314 (self.spin, iter, np.log10(K), ESI * Hartree))
1315 elif self.logging == 2:
1316 print(" UO-%d: %2d %5.1f %20.5f " %
1317 (self.spin, iter, np.log10(K), ESI * Hartree) +
1318 lsmethod)
1320 if ((K < self.uomaxres or
1321 K < self.wfs.eigensolver.error * self.uorelres) or noise or
1322 failed) and not self.maxuoiter < 0:
1323 break