Coverage for gpaw/directmin/etdm_fdpw.py: 92%
602 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-14 00:18 +0000
1"""
2A class for finding optimal orbitals corresponding to a stationary point of
3the energy functional using direct optimization and exponential transformation
4in FD and PW modes.
6It can be used with Kohn-Sham and hybrid (exact exchange) functionals, and
7can include Perdew-Zunger self-interaction correction (PZ-SIC) in the
8calculations.
10Ground state as well as variational excited state calculations can be
11performed. Ground state calculations involve minimization of the energy in a
12direction tangent to the orbitals without the exponential transformation
13(direct minimization). For excited state calculations, the energy is optimized
14by converging on a saddle point, which involves an inner loop using the
15exponential transformation (direct optimization). PZ-SIC requires an
16additional inner loop to minimize the energy with respect to unitary
17transformation of the occupied orbitals (inner loop localization).
19GPAW Implementation of direct optimization with preconditioned quasi-Newton
20algorithms and maximum overlap method (DO-MOM) for excited state calculations
21FD and PW modes:
23 J. Chem. Theory Comput. 17, 5034–5049 (2021) :doi:10.1021/acs.jctc.1c00157
24 arXiv:2102.06542 [physics.comp-ph]
25"""
27import time
29from ase.parallel import parprint
30from ase.units import Hartree
31from ase.utils import basestring
32import numpy as np
34from gpaw.directmin import search_direction, line_search_algorithm
35from gpaw.directmin.fdpw.etdm_inner_loop import ETDMInnerLoop
36from gpaw.directmin.fdpw.pz_localization import PZLocalization
37from gpaw.directmin.functional.fdpw import get_functional
38from gpaw.directmin.locfunc.localize_orbitals import localize_orbitals
39from gpaw.directmin.tools import get_n_occ, sort_orbitals_according_to_occ
40from gpaw.utilities import unpack_hermitian
41from gpaw.xc import xc_string_to_dict
42from gpaw.xc.hybrid import HybridXC
45class FDPWETDM:
47 def __init__(self,
48 excited_state=False,
49 searchdir_algo=None,
50 linesearch_algo='max-step',
51 use_prec=True,
52 functional='ks',
53 need_init_orbs=None,
54 need_localization=True,
55 localizationtype=None,
56 localizationseed=None,
57 localization_tol=None,
58 maxiter_pz_localization=50,
59 maxiter_inner_loop=100,
60 max_step_inner_loop=0.2,
61 grad_tol_pz_localization=5.0e-4,
62 grad_tol_inner_loop=5.0e-4,
63 restartevery_iloop_notconverged=3,
64 restart_canonical=True,
65 momevery=10,
66 printinnerloop=False,
67 blocksize=1,
68 converge_unocc=False,
69 maxiter_unocc=333
70 ):
71 """Class for direct orbital optimization in FD and PW modes.
73 Parameters
74 ----------
75 excited_state: bool
76 If False (default), perform a minimization in the tangent space of
77 orbitals (ground state calculation), and set need_init_orbs to
78 True. Otherwise, perform outer loop minimization in the tangent
79 space and inner loop optimization with exponential transformation
80 (excited state calculation), and set need_init_orbs to False.
81 searchdir_algo: str, dict or instance
82 Search direction algorithm for the outer loop minimization. Can be
83 one of the algorithms available in sd_etdm.py:
84 'sd': Steepest descent
85 'fr-cg': Fletcher-Reeves conjugate gradient
86 'l-bfgs': Limited-memory BFGS (default)
87 'l-bfgs-p': Limited-memory BFGS with preconditioner presented
88 in :doi:`10.1016/j.cpc.2021.108047`
89 'l-sr1p': Limited-memory SR1 algorithm presented in
90 :doi:`10.1021/acs.jctc.0c00597`
91 The default memory for 'l-bfgs'/'l-bfgs-p' and 'l-sr1p' is 3 and
92 20, respectively, and can be changed by supplying a dictionary:
93 {'name': name, 'memory': memory}, where name should be 'l-bfgs',
94 'l-bfgs-p' or 'l-sr1p' and memory should be an int.
95 linesearch_algo: str, dict or instance
96 Line search algorithm for the outer loop minimization. Can be one
97 of the algorithms available in ls_etdm.py:
98 'max-step': The quasi-Newton step is scaled if it exceeds a
99 maximum step length (default). The default maximum step
100 length is 0.20, and can be changed by supplying a
101 dictionary: {'name': 'max-step', 'max_step': max_step},
102 where max_step should be a float.
103 'swc-awc': Line search with Wolfe conditions
104 use_prec: bool
105 If True (default) use a preconditioner. For the outer loop
106 minimization, the preconditioner is the inverse kinetic energy
107 operator. For the inner loop optimization, the preconditioner is
108 the inverse of a diagonal approximation of the Hessian (see
109 :doi:`10.1021/j100322a012`) apart for 'l-bfgs-p', which uses the
110 composite preconditioner presented in
111 :doi:`10.1016/j.cpc.2021.108047`.
112 functional: str, dict or instance
113 Type of functional. Can be one of:
114 'ks': The functional as specified in the GPAW calculator is
115 used (default)
116 'pz-sic': Apply the Perdew-Zunger self-interaction correction
117 on top of the functional as specified in the GPAW
118 calculator. Dy default full SIC is applied. A scaling
119 factor for SIC can be given by supplying a dictionary:
120 functional={'name': 'pz-sic', 'scaling_factor': (a, a)},
121 where a is the scaling factor (float).
122 need_init_orbs: bool
123 If True (default when excited_state is False), obtain initial
124 orbitals from eigendecomposition of the Hamiltonian matrix. If
125 False (default when excited_state is True), use orbitals stored in
126 wfs object as initial guess.
127 need_localization: bool
128 If True (default), localize initial guess orbitals. Requires a
129 specification of localizationtype, otherwise it is set to False.
130 Recommended for calculations with PZ-SIC.
131 localizationtype: str
132 Method for localizing the initial guess orbitals. Can be one of:
133 'pz': Unitary optimization among occupied orbitals (subspace
134 optimization) with PZ-SIC
135 'ks':
136 'er': Edmiston-Ruedenberg localization
137 'pm': Pipek-Mezey localization (recommended for PZ-SIC)
138 'fb' Foster-Boys localization
139 Default is None, meaning that no localization is performed.
140 localizationseed: int
141 Seed for Edmiston-Ruedenberg, Pipek-Mezey or Foster-Boys
142 localization. Default is None (no seed is used).
143 localization_tol: float
144 Tolerance for convergence of the localization. If not specified,
145 the following default values will be used:
146 'pz': 5.0e-4
147 'ks': 5.0e-4
148 'er': 5.0e-5
149 'pm': 1.0e-10
150 'fb': 1.0e-10
151 maxiter_pz_localization: int
152 Maximum number of iterations for PZ-SIC inner loop localization.
153 maxiter_inner_loop: int
154 Maximum number of iterations of inner loop optimization for
155 excited state calculations. If the maximum number of inner loop
156 iterations is exceeded, the optimization moves on with the outer
157 loop step.
158 max_step_inner_loop: float
159 Maximum step length of inner loop optimization for excited state
160 calculations. The inner loop optimization uses the 'l-sr1p' search
161 direction. Default is 0.20.
162 grad_tol_pz_localization: float
163 Tolerance on the norm of the gradient for convergence of the
164 PZ-SIC inner loop localization.
165 grad_tol_inner_loop: float
166 Tolerance on the norm of the gradient for convergence of the
167 inner loop optimization for excited state calculations.
168 restartevery_iloop_notconverged: int
169 Number of iterations of the outer loop after which the calculation
170 is restarted if the inner loop optimization for excited states is
171 not converged.
172 restart_canonical: bool
173 If True (default) restart the calculations using orbitals from the
174 eigedecomposition of the Hamiltonian matrix if the inner loop
175 optimization does not converge or MOM detects variational
176 collapse. Otherwise, the optimal orbitals are used.
177 momevery: int
178 MOM is applied every 'momevery' iterations of the inner loop
179 optimization for excited states.
180 printinnerloop: bool
181 If True, print the iterations of the inner loop optimization for
182 excited states to standard output. Default is False.
183 blocksize: int
184 Blocksize for base eigensolver class.
185 converge_unocc: bool
186 If True, converge also the unoccupied orbitals after convergence
187 of the occupied orbitals. Default is False.
188 maxiter_unocc: int
189 Maximum number of iterations for convergence of the unoccupied
190 orbitals.
191 """
193 self.error = np.inf
194 self.blocksize = blocksize
196 self.name = 'etdm-fdpw'
197 self.sda = searchdir_algo
198 self.lsa = linesearch_algo
199 self.use_prec = use_prec
200 self.func_settings = functional
201 self.need_init_orbs = need_init_orbs
202 self.localizationtype = localizationtype
203 self.localizationseed = localizationseed
204 self.need_localization = need_localization
205 self.localization_tol = localization_tol
206 self.maxiter_pz_localization = maxiter_pz_localization
207 self.maxiter_inner_loop = maxiter_inner_loop
208 self.max_step_inner_loop = max_step_inner_loop
209 self.grad_tol_pz_localization = grad_tol_pz_localization
210 self.grad_tol_inner_loop = grad_tol_inner_loop
211 self.printinnerloop = printinnerloop
212 self.converge_unocc = converge_unocc
213 self.maxiter_unocc = maxiter_unocc
214 self.restartevery_iloop_notconverged = restartevery_iloop_notconverged
215 self.restart_canonical = restart_canonical
216 self.momevery = momevery
217 self.excited_state = excited_state
218 self.check_inputs_and_init_search_algo()
220 self.eg_count_iloop = 0
221 self.total_eg_count_iloop = 0
222 self.eg_count_outer_iloop = 0
223 self.total_eg_count_outer_iloop = 0
224 self.initial_random = True
226 # for mom
227 self.initial_occupation_numbers = None
229 self.eg_count = 0
230 self.etotal = 0.0
231 self.globaliters = 0
232 self.need_init_odd = True
233 self.initialized = False
235 self.gpaw_new = False
237 def check_inputs_and_init_search_algo(self):
238 defaults = self.set_defaults()
239 if self.need_init_orbs is None:
240 self.need_init_orbs = defaults['need_init_orbs']
241 if self.localizationtype is None:
242 self.need_localization = False
244 if self.sda is None:
245 self.sda = 'LBFGS'
246 self.searchdir_algo = search_direction(self.sda)
248 self.line_search = line_search_algorithm(
249 self.lsa, self.evaluate_phi_and_der_phi,
250 self.searchdir_algo)
252 if isinstance(self.func_settings, basestring):
253 self.func_settings = xc_string_to_dict(self.func_settings)
255 def set_defaults(self):
256 if self.excited_state:
257 return {'need_init_orbs': False}
258 else:
259 return {'need_init_orbs': True}
261 def __repr__(self):
263 sda_name = self.searchdir_algo.name
264 lsa_name = self.line_search.name
265 if isinstance(self.func_settings, basestring):
266 func_name = self.func_settings
267 else:
268 func_name = self.func_settings['name']
270 sds = {'sd': 'Steepest Descent',
271 'fr-cg': 'Fletcher-Reeves conj. grad. method',
272 'l-bfgs': 'L-BFGS algorithm',
273 'l-bfgs-p': 'L-BFGS algorithm with preconditioning',
274 'l-sr1p': 'Limited-memory SR1P algorithm'}
276 lss = {'max-step': 'step size equals one',
277 'swc-awc': 'Inexact line search based on cubic interpolation,\n'
278 ' strong and approximate Wolfe '
279 'conditions'}
281 repr_string = 'Direct minimisation using exponential ' \
282 'transformation.\n'
283 repr_string += ' ' \
284 'Search ' \
285 'direction: {}\n'.format(sds[sda_name])
286 repr_string += ' ' \
287 'Line ' \
288 'search: {}\n'.format(lss[lsa_name])
289 repr_string += ' ' \
290 'Preconditioning: {}\n'.format(self.use_prec)
291 repr_string += ' ' \
292 'Orbital-density self-interaction ' \
293 'corrections: {}\n'.format(func_name)
294 repr_string += ' ' \
295 'WARNING: do not use it for metals as ' \
296 'occupation numbers are\n' \
297 ' ' \
298 'not found variationally\n'
300 return repr_string
302 def reset(self, need_init_odd=True):
303 self.initialized = False
304 self.need_init_odd = need_init_odd
305 self.searchdir_algo.reset()
307 def todict(self):
308 """
309 Convert to dictionary, needs for saving and loading gpw
310 :return:
311 """
312 return {'name': 'etdm-fdpw',
313 'searchdir_algo': self.searchdir_algo.todict(),
314 'linesearch_algo': self.line_search.todict(),
315 'use_prec': self.use_prec,
316 'functional': self.func_settings,
317 'need_init_orbs': self.need_init_orbs,
318 'localizationtype': self.localizationtype,
319 'localizationseed': self.localizationseed,
320 'need_localization': self.need_localization,
321 'localization_tol': self.localization_tol,
322 'maxiter_pz_localization': self.maxiter_pz_localization,
323 'maxiter_inner_loop': self.maxiter_inner_loop,
324 'max_step_inner_loop': self.max_step_inner_loop,
325 'momevery': self.momevery,
326 'grad_tol_pz_localization': self.grad_tol_pz_localization,
327 'grad_tol_inner_loop': self.grad_tol_inner_loop,
328 'restartevery_iloop_notconverged':
329 self.restartevery_iloop_notconverged,
330 'restart_canonical': self.restart_canonical,
331 'printinnerloop': self.printinnerloop,
332 'blocksize': self.blocksize,
333 'converge_unocc': self.converge_unocc,
334 'maxiter_unocc': self.maxiter_unocc,
335 'excited_state': self.excited_state
336 }
338 def initialize_dm_helper(self, wfs, ham, dens, log):
339 self.initialize_eigensolver(wfs, ham)
340 self.initialize_orbitals(wfs, ham)
342 if not wfs.read_from_file_init_wfs_dm or self.excited_state:
343 wfs.calculate_occupation_numbers(dens.fixed)
345 self.initial_sort_orbitals(wfs)
347 # localize orbitals?
348 self.localize(wfs, dens, ham, log)
350 # MOM
351 self.initialize_mom_reference_orbitals(wfs, dens)
353 # initialize search direction, line search and inner loops
354 self.initialize_dm(wfs, dens, ham)
356 def initialize_eigensolver(self, wfs, ham):
357 """
358 Initialize base eigensolver class
360 :param wfs:
361 :param ham:
362 :return:
363 """
364 if isinstance(ham.xc, HybridXC):
365 self.blocksize = wfs.bd.mynbands
367 if self.blocksize is None:
368 if wfs.mode == 'pw':
369 S = wfs.pd.comm.size
370 # Use a multiple of S for maximum efficiency
371 self.blocksize = int(np.ceil(10 / S)) * S
372 else:
373 self.blocksize = 10
375 from gpaw.eigensolvers.eigensolver import Eigensolver
376 self.eigensolver = Eigensolver(keep_htpsit=False,
377 blocksize=self.blocksize)
378 self.eigensolver.initialize(wfs)
380 def initialize_dm(
381 self, wfs, dens, ham, converge_unocc=False):
383 """
384 initialize search direction algorithm,
385 line search method, SIC corrections
387 :param wfs:
388 :param dens:
389 :param ham:
390 :param converge_unocc:
391 :return:
392 """
394 self.searchdir_algo.reset()
396 self.dtype = wfs.dtype
397 self.n_kps = wfs.kd.nibzkpts
398 # dimensionality, number of state to be converged
399 self.dimensions = {}
400 for kpt in wfs.kpt_u:
401 bd = self.eigensolver.bd
402 nocc = get_n_occ(kpt)[0]
403 if converge_unocc:
404 assert nocc < bd.nbands, \
405 'Please add empty bands in order to converge the' \
406 ' unoccupied orbitals'
407 dim = bd.nbands - nocc
408 elif self.excited_state:
409 dim = bd.nbands
410 else:
411 dim = nocc
413 k = self.n_kps * kpt.s + kpt.q
414 self.dimensions[k] = dim
416 if self.use_prec:
417 self.prec = wfs.make_preconditioner(1)
418 else:
419 self.prec = None
421 self.iters = 0
422 self.alpha = 1.0 # step length
423 self.phi_2i = [None, None] # energy at last two iterations
424 self.der_phi_2i = [None, None] # energy gradient w.r.t. alpha
425 self.grad_knG = None
427 # odd corrections
428 if self.need_init_odd:
429 self.odd = get_functional(self.func_settings, wfs, dens, ham)
430 self.e_sic = 0.0
432 if 'SIC' in self.odd.name:
433 self.iloop = PZLocalization(
434 self.odd, wfs, self.maxiter_pz_localization,
435 g_tol=self.grad_tol_pz_localization)
436 else:
437 self.iloop = None
439 if self.excited_state:
440 self.outer_iloop = ETDMInnerLoop(
441 self.odd, wfs, 'all', self.maxiter_inner_loop,
442 self.max_step_inner_loop, g_tol=self.grad_tol_inner_loop,
443 useprec=True, momevery=self.momevery)
444 # if you have inner-outer loop then need to have
445 # U matrix of the same dimensionality in both loops
446 if 'SIC' in self.odd.name:
447 for kpt in wfs.kpt_u:
448 k = self.n_kps * kpt.s + kpt.q
449 self.iloop.U_k[k] = self.outer_iloop.U_k[k].copy()
450 else:
451 self.outer_iloop = None
453 self.total_eg_count_iloop = 0
454 self.total_eg_count_outer_iloop = 0
456 self.initialized = True
457 wfs.read_from_file_init_wfs_dm = False
459 def initial_sort_orbitals(self, wfs):
460 occ_name = getattr(wfs.occupations, "name", None)
461 if occ_name == 'mom' and self.globaliters == 0:
462 update_mom = True
463 self.initial_occupation_numbers = \
464 wfs.occupations.numbers.copy()
465 else:
466 update_mom = False
467 sort_orbitals_according_to_occ(wfs, update_mom=update_mom)
469 def iterate(self, ham, wfs, dens, log, converge_unocc=False):
470 """
471 One iteration of outer loop direct minimization
473 :param ham:
474 :param wfs:
475 :param dens:
476 :param log:
477 :return:
478 """
480 n_kps = self.n_kps
481 psi_copy = {}
482 alpha = self.alpha
483 phi_2i = self.phi_2i
484 der_phi_2i = self.der_phi_2i
486 wfs.timer.start('Direct Minimisation step')
488 if self.iters == 0:
489 # calculate gradients
490 if not converge_unocc:
491 phi_2i[0], grad_knG = \
492 self.get_energy_and_tangent_gradients(ham, wfs, dens)
493 else:
494 phi_2i[0], grad_knG = \
495 self.get_energy_and_tangent_gradients_unocc(ham, wfs)
496 else:
497 grad_knG = self.grad_knG
499 wfs.timer.start('Get Search Direction')
500 for kpt in wfs.kpt_u:
501 k = n_kps * kpt.s + kpt.q
502 if not converge_unocc:
503 psi_copy[k] = kpt.psit_nG.copy()
504 else:
505 n_occ = get_n_occ(kpt)[0]
506 dim = self.dimensions[k]
507 psi_copy[k] = kpt.psit_nG[n_occ:n_occ + dim].copy()
509 p_knG = self.searchdir_algo.update_data(
510 wfs, psi_copy, grad_knG, precond=self.prec,
511 dimensions=self.dimensions)
512 self.project_search_direction(wfs, p_knG)
513 wfs.timer.stop('Get Search Direction')
515 # recalculate derivative with new search direction
516 # as we used preconditiner before
517 # here we project search direction on prec. gradients,
518 # but should be just grad. But, it seems also works fine
519 der_phi_2i[0] = 0.0
520 for kpt in wfs.kpt_u:
521 k = n_kps * kpt.s + kpt.q
522 for i, g in enumerate(grad_knG[k]):
523 if kpt.f_n[i] > 1.0e-10:
524 der_phi_2i[0] += self.dot(
525 wfs, g, p_knG[k][i], kpt, addpaw=False).item().real
526 der_phi_2i[0] = wfs.kd.comm.sum_scalar(der_phi_2i[0])
528 alpha, phi_alpha, der_phi_alpha, grad_knG = \
529 self.line_search.step_length_update(
530 psi_copy, p_knG, wfs, ham, dens, converge_unocc,
531 phi_0=phi_2i[0], der_phi_0=der_phi_2i[0],
532 phi_old=phi_2i[1], der_phi_old=der_phi_2i[1],
533 alpha_max=3.0, alpha_old=alpha, kpdescr=wfs.kd)
534 self.alpha = alpha
535 self.grad_knG = grad_knG
537 # and 'shift' phi, der_phi for the next iteration
538 phi_2i[1], der_phi_2i[1] = phi_2i[0], der_phi_2i[0]
539 phi_2i[0], der_phi_2i[0] = phi_alpha, der_phi_alpha,
541 self.iters += 1
542 if not converge_unocc:
543 self.globaliters += 1
544 wfs.timer.stop('Direct Minimisation step')
545 return phi_2i[0], self.error
547 def update_ks_energy(self, ham, wfs, dens, updateproj=True):
548 """Update Kohn-Sham energy.
550 It assumes the temperature is zero K.
551 """
553 if updateproj:
554 # calc projectors
555 with wfs.timer('projections'):
556 for kpt in wfs.kpt_u:
557 wfs.pt.integrate(kpt.psit_nG, kpt.P_ani, kpt.q)
559 dens.update(wfs)
560 ham.update(dens, wfs, False)
562 return ham.get_energy(0.0, wfs, False)
564 def evaluate_phi_and_der_phi(self, psit_k, search_dir, alpha, wfs, ham,
565 dens, converge_unocc, phi=None, grad_k=None):
566 """
567 phi = E(x_k + alpha_k*p_k)
568 der_phi = grad_alpha E(x_k + alpha_k*p_k) cdot p_k
569 :return: phi, der_phi # floats
570 """
572 if phi is None or grad_k is None:
573 alpha1 = np.array([alpha])
574 wfs.world.broadcast(alpha1, 0)
575 alpha = alpha1[0]
577 x_knG = \
578 {k: psit_k[k] +
579 alpha * search_dir[k] for k in psit_k.keys()}
580 if not converge_unocc:
581 phi, grad_k = self.get_energy_and_tangent_gradients(
582 ham, wfs, dens, psit_knG=x_knG)
583 else:
584 phi, grad_k = self.get_energy_and_tangent_gradients_unocc(
585 ham, wfs, x_knG)
587 der_phi = 0.0
588 for kpt in wfs.kpt_u:
589 k = self.n_kps * kpt.s + kpt.q
590 for i, g in enumerate(grad_k[k]):
591 if not converge_unocc and kpt.f_n[i] > 1.0e-10:
592 der_phi += self.dot(
593 wfs, g, search_dir[k][i], kpt,
594 addpaw=False).item().real
595 else:
596 der_phi += self.dot(
597 wfs, g, search_dir[k][i], kpt,
598 addpaw=False).item().real
599 der_phi = wfs.kd.comm.sum_scalar(der_phi)
601 return phi, der_phi, grad_k
603 def get_energy_and_tangent_gradients(
604 self, ham, wfs, dens, psit_knG=None, updateproj=True):
606 """
607 calculate energy for a given wfs, gradient dE/dpsi
608 and then project gradient on tangent space to psi
610 :param ham:
611 :param wfs:
612 :param dens:
613 :param psit_knG:
614 :return:
615 """
617 n_kps = self.n_kps
618 if psit_knG is not None:
619 for kpt in wfs.kpt_u:
620 k = n_kps * kpt.s + kpt.q
621 kpt.psit_nG[:] = psit_knG[k].copy()
622 wfs.orthonormalize(kpt)
623 elif not wfs.orthonormalized:
624 wfs.orthonormalize()
625 if not self.excited_state:
626 energy = self.update_ks_energy(
627 ham, wfs, dens, updateproj=updateproj)
628 grad = self.get_gradients_2(ham, wfs)
630 if 'SIC' in self.odd.name:
631 e_sic = 0.0
632 if self.iters > 0:
633 self.run_inner_loop(ham, wfs, dens, grad_knG=grad)
634 else:
635 for kpt in wfs.kpt_u:
636 e_sic += self.odd.get_energy_and_gradients_kpt(
637 wfs, kpt, grad, self.iloop.U_k, add_grad=True)
638 self.e_sic = wfs.kd.comm.sum_scalar(e_sic)
639 ham.get_energy(0.0, wfs, kin_en_using_band=False,
640 e_sic=self.e_sic)
641 energy += self.e_sic
642 else:
643 grad = {}
644 n_kps = self.n_kps
645 for kpt in wfs.kpt_u:
646 grad[n_kps * kpt.s + kpt.q] = np.zeros_like(kpt.psit_nG[:])
647 self.run_inner_loop(ham, wfs, dens, grad_knG=grad)
648 energy = self.etotal
650 self.project_gradient(wfs, grad)
651 self.error = self.error_eigv(wfs, grad)
652 self.eg_count += 1
653 return energy, grad
655 def get_gradients_2(self, ham, wfs, scalewithocc=True):
657 """
658 calculate gradient dE/dpsi
659 :return: H |psi_i>
660 """
662 grad_knG = {}
663 n_kps = self.n_kps
665 for kpt in wfs.kpt_u:
666 grad_knG[n_kps * kpt.s + kpt.q] = \
667 self.get_gradients_from_one_k_point_2(
668 ham, wfs, kpt, scalewithocc)
670 return grad_knG
672 def get_gradients_from_one_k_point_2(
673 self, ham, wfs, kpt, scalewithocc=True):
674 """
675 calculate gradient dE/dpsi for one k-point
676 :return: H |psi_i>
677 """
678 nbands = wfs.bd.mynbands
679 Hpsi_nG = wfs.empty(nbands, q=kpt.q)
680 wfs.apply_pseudo_hamiltonian(kpt, ham, kpt.psit_nG, Hpsi_nG)
682 c_axi = {}
683 if self.gpaw_new:
684 dH_asii = ham.potential.dH_asii
685 for a, P_xi in kpt.P_ani.items():
686 dH_ii = dH_asii[a][kpt.s]
687 c_xi = np.dot(P_xi, dH_ii)
688 c_axi[a] = c_xi
689 else:
690 for a, P_xi in kpt.P_ani.items():
691 dH_ii = unpack_hermitian(ham.dH_asp[a][kpt.s])
692 c_xi = np.dot(P_xi, dH_ii)
693 c_axi[a] = c_xi
695 # not sure about this:
696 ham.xc.add_correction(
697 kpt, kpt.psit_nG, Hpsi_nG, kpt.P_ani, c_axi, n_x=None,
698 calculate_change=False)
699 # add projectors to the H|psi_i>
700 wfs.pt.add(Hpsi_nG, c_axi, kpt.q)
701 # scale with occupation numbers
702 if scalewithocc:
703 for i, f in enumerate(kpt.f_n):
704 Hpsi_nG[i] *= f
705 return Hpsi_nG
707 def project_gradient(self, wfs, p_knG):
708 """
709 project gradient dE/dpsi on tangent space at psi
710 See Eq.(22) and minimization algorithm p. 12 in
711 arXiv:2102.06542v1 [physics.comp-ph]
712 :return: H |psi_i>
713 """
715 n_kps = self.n_kps
716 for kpt in wfs.kpt_u:
717 kpoint = n_kps * kpt.s + kpt.q
718 self.project_gradient_for_one_k_point(wfs, p_knG[kpoint], kpt)
720 def project_gradient_for_one_k_point(self, wfs, p_nG, kpt):
721 """
722 project gradient dE/dpsi on tangent space at psi
723 for one k-point.
724 See Eq.(22) and minimization algorithm p. 12 in
725 arXiv:2102.06542v1 [physics.comp-ph]
726 :return: H |psi_i>
727 """
729 k = self.n_kps * kpt.s + kpt.q
730 n_occ = self.dimensions[k]
731 psc = wfs.integrate(p_nG[:n_occ], kpt.psit_nG[:n_occ], True)
732 psc = 0.5 * (psc.conj() + psc.T)
733 s_psit_nG = self.apply_S(wfs, kpt.psit_nG, kpt, kpt.P_ani)
734 p_nG[:n_occ] -= np.tensordot(psc, s_psit_nG[:n_occ], axes=1)
736 def project_search_direction(self, wfs, p_knG):
738 """
739 Project search direction on tangent space at psi
740 it is slighlt different from project grad
741 as it doesn't apply overlap matrix because of S^{-1}
743 :param wfs:
744 :param p_knG:
745 :return:
746 """
748 for kpt in wfs.kpt_u:
749 k = self.n_kps * kpt.s + kpt.q
750 n_occ = self.dimensions[k]
751 psc = self.dot(
752 wfs, p_knG[k][:n_occ], kpt.psit_nG[:n_occ], kpt, addpaw=False)
753 psc = 0.5 * (psc.conj() + psc.T)
754 p_knG[k][:n_occ] -= np.tensordot(psc, kpt.psit_nG[:n_occ], axes=1)
756 def apply_S(self, wfs, psit_nG, kpt, proj_psi=None):
757 """
758 apply overlap matrix
760 :param wfs:
761 :param psit_nG:
762 :param kpt:
763 :param proj_psi:
764 :return:
765 """
767 if proj_psi is None:
768 proj_psi = wfs.pt.dict(shape=wfs.bd.mynbands)
769 wfs.pt.integrate(psit_nG, proj_psi, kpt.q)
771 s_axi = {}
772 for a, P_xi in proj_psi.items():
773 dO_ii = wfs.setups[a].dO_ii
774 s_xi = np.dot(P_xi, dO_ii)
775 s_axi[a] = s_xi
777 new_psi_nG = psit_nG.copy()
778 wfs.pt.add(new_psi_nG, s_axi, kpt.q)
780 return new_psi_nG
782 def dot(self, wfs, psi_1, psi_2, kpt, addpaw=True):
783 """
784 dor product between two arrays psi_1 and psi_2
786 :param wfs:
787 :param psi_1:
788 :param psi_2:
789 :param kpt:
790 :param addpaw:
791 :return:
792 """
794 dot_prod = wfs.integrate(psi_1, psi_2, global_integral=True)
795 if not addpaw:
796 if len(psi_1.shape) == 4 or len(psi_1.shape) == 2:
797 sum_dot = dot_prod
798 else:
799 sum_dot = np.asarray([[dot_prod]])
801 return sum_dot
803 def dS(a, P_ni):
804 """
805 apply PAW
806 :param a:
807 :param P_ni:
808 :return:
809 """
810 return np.dot(P_ni, wfs.setups[a].dO_ii)
812 if len(psi_1.shape) == 3 or len(psi_1.shape) == 1:
813 ndim = 1
814 else:
815 ndim = psi_1.shape[0]
817 P1_ai = wfs.pt.dict(shape=ndim)
818 P2_ai = wfs.pt.dict(shape=ndim)
819 wfs.pt.integrate(psi_1, P1_ai, kpt.q)
820 wfs.pt.integrate(psi_2, P2_ai, kpt.q)
821 if ndim == 1:
822 if self.dtype == complex:
823 paw_dot_prod = np.array([[0.0 + 0.0j]])
824 else:
825 paw_dot_prod = np.array([[0.0]])
827 for a in P1_ai.keys():
828 paw_dot_prod += np.dot(dS(a, P2_ai[a]), P1_ai[a].T.conj())
829 else:
830 paw_dot_prod = np.zeros_like(dot_prod)
831 for a in P1_ai.keys():
832 paw_dot_prod += np.dot(dS(a, P2_ai[a]), P1_ai[a].T.conj()).T
833 paw_dot_prod = np.ascontiguousarray(paw_dot_prod)
834 wfs.gd.comm.sum(paw_dot_prod)
835 if len(psi_1.shape) == 4 or len(psi_1.shape) == 2:
836 sum_dot = dot_prod + paw_dot_prod
837 else:
838 sum_dot = [[dot_prod]] + paw_dot_prod
840 return sum_dot
842 def error_eigv(self, wfs, grad_knG):
843 """
844 calcualte norm of the gradient vector
845 (residual)
847 :param wfs:
848 :param grad_knG:
849 :return:
850 """
852 n_kps = wfs.kd.nibzkpts
853 norm = [0.0]
854 for kpt in wfs.kpt_u:
855 k = n_kps * kpt.s + kpt.q
856 for i, f in enumerate(kpt.f_n):
857 if f > 1.0e-10:
858 a = self.dot(
859 wfs, grad_knG[k][i] / f, grad_knG[k][i] / f, kpt,
860 addpaw=False).item() * f
861 a = a.real
862 norm.append(a)
863 error = sum(norm)
864 error = wfs.kd.comm.sum_scalar(error)
866 return error.real
868 def get_canonical_representation(self, ham, wfs, rewrite_psi=True):
869 """
870 choose orbitals which diagonalize the hamiltonain matrix
872 <psi_i| H |psi_j>
874 For SIC, one diagonalizes L_{ij} = <psi_i| H + V_{j} |psi_j>
875 for occupied subspace and
876 <psi_i| H |psi_j> for unoccupied.
878 :param ham:
879 :param wfs:
880 :param rewrite_psi:
881 :return:
882 """
883 self.choose_optimal_orbitals(wfs)
885 scalewithocc = not self.excited_state
887 grad_knG = self.get_gradients_2(ham, wfs, scalewithocc=scalewithocc)
888 if 'SIC' in self.odd.name:
889 for kpt in wfs.kpt_u:
890 self.odd.get_energy_and_gradients_kpt(
891 wfs, kpt, grad_knG, self.iloop.U_k,
892 add_grad=True, scalewithocc=scalewithocc)
894 for kpt in wfs.kpt_u:
895 # Separate diagonalization for occupied
896 # and unoccupied subspaces
897 bd = self.eigensolver.bd
898 k = self.n_kps * kpt.s + kpt.q
899 n_occ = get_n_occ(kpt)[0]
900 dim = bd.nbands - n_occ
901 if scalewithocc:
902 scale = kpt.f_n[:n_occ]
903 else:
904 scale = 1.0
906 if scalewithocc:
907 grad_knG[k][n_occ:n_occ + dim] = \
908 self.get_gradients_unocc_kpt(ham, wfs, kpt)
909 lamb = wfs.integrate(kpt.psit_nG[:], grad_knG[k][:], True)
910 lamb1 = (lamb[:n_occ, :n_occ] +
911 lamb[:n_occ, :n_occ].T.conj()) / 2.0
912 lumo = (lamb[n_occ:, n_occ:] +
913 lamb[n_occ:, n_occ:].T.conj()) / 2.0
915 # Diagonal elements Lagrangian matrix
916 lo_nn = np.diagonal(lamb1).real / scale
917 lu_nn = np.diagonal(lumo).real / 1.0
919 # Diagonalize occupied subspace
920 if n_occ != 0:
921 evals, lamb1 = np.linalg.eigh(lamb1)
922 wfs.gd.comm.broadcast(evals, 0)
923 wfs.gd.comm.broadcast(lamb1, 0)
924 lamb1 = lamb1.T
925 kpt.eps_n[:n_occ] = evals[:n_occ] / scale
927 # Diagonalize unoccupied subspace
928 evals_lumo, lumo = np.linalg.eigh(lumo)
929 wfs.gd.comm.broadcast(evals_lumo, 0)
930 wfs.gd.comm.broadcast(lumo, 0)
931 lumo = lumo.T
932 kpt.eps_n[n_occ:n_occ + dim] = evals_lumo.real
934 if rewrite_psi: # Only for SIC
935 kpt.psit_nG[:n_occ] = np.tensordot(
936 lamb1.conj(), kpt.psit_nG[:n_occ], axes=1)
937 kpt.psit_nG[n_occ:n_occ + dim] = np.tensordot(
938 lumo.conj(), kpt.psit_nG[n_occ:n_occ + dim], axes=1)
940 wfs.pt.integrate(kpt.psit_nG, kpt.P_ani, kpt.q)
942 if 'SIC' in self.odd.name:
943 self.odd.lagr_diag_s[k] = np.append(lo_nn, lu_nn)
945 del grad_knG
947 def get_gradients_unocc_kpt(self, ham, wfs, kpt):
948 """
949 calculate gradient vector for unoccupied orbitals
951 :param ham:
952 :param wfs:
953 :param kpt:
954 :return:
955 """
956 orbital_dependent = ham.xc.orbital_dependent
958 n_occ = get_n_occ(kpt)[0]
959 bd = self.eigensolver.bd
960 if not orbital_dependent:
961 dim = bd.nbands - n_occ
962 n0 = n_occ
963 else:
964 dim = bd.nbands
965 n0 = 0
967 # calculate gradients:
968 psi = kpt.psit_nG[n0:n0 + dim].copy()
969 P1_ai = wfs.pt.dict(shape=dim)
970 wfs.pt.integrate(psi, P1_ai, kpt.q)
971 Hpsi_nG = wfs.empty(dim, q=kpt.q)
972 wfs.apply_pseudo_hamiltonian(kpt, ham, psi, Hpsi_nG)
973 c_axi = {}
974 for a, P_xi in P1_ai.items():
975 if self.gpaw_new:
976 dH_ii = ham.potential.dH_asii[a][kpt.s]
977 else:
978 dH_ii = unpack_hermitian(ham.dH_asp[a][kpt.s])
979 c_xi = np.dot(P_xi, dH_ii)
980 c_axi[a] = c_xi
981 # not sure about this:
982 ham.xc.add_correction(kpt, psi, Hpsi_nG,
983 P1_ai, c_axi, n_x=None,
984 calculate_change=False)
985 # add projectors to the H|psi_i>
986 wfs.pt.add(Hpsi_nG, c_axi, kpt.q)
988 if orbital_dependent:
989 Hpsi_nG = Hpsi_nG[n_occ:n_occ + dim]
991 return Hpsi_nG
993 def get_energy_and_tangent_gradients_unocc(self, ham, wfs, psit_knG=None):
994 """
995 calculate energy and trangent gradients of
996 unooccupied orbitals
998 :param ham:
999 :param wfs:
1000 :param psit_knG:
1001 :return:
1002 """
1003 wfs.timer.start('Gradient unoccupied orbitals')
1004 n_kps = self.n_kps
1005 if psit_knG is not None:
1006 for kpt in wfs.kpt_u:
1007 k = n_kps * kpt.s + kpt.q
1008 # find lumo
1009 n_occ = get_n_occ(kpt)[0]
1010 dim = self.dimensions[k]
1011 kpt.psit_nG[n_occ:n_occ + dim] = psit_knG[k].copy()
1012 wfs.orthonormalize(kpt)
1013 elif not wfs.orthonormalized:
1014 wfs.orthonormalize()
1016 grad = {}
1017 energy_t = 0.0
1018 error_t = 0.0
1020 for kpt in wfs.kpt_u:
1021 n_occ = get_n_occ(kpt)[0]
1022 k = n_kps * kpt.s + kpt.q
1023 dim = self.dimensions[k]
1024 Hpsi_nG = self.get_gradients_unocc_kpt(ham, wfs, kpt)
1025 grad[k] = Hpsi_nG.copy()
1027 # calculate energy
1028 psi = kpt.psit_nG[:n_occ + dim].copy()
1029 wfs.pt.integrate(kpt.psit_nG, kpt.P_ani, kpt.q)
1030 lamb = wfs.integrate(psi, Hpsi_nG, global_integral=True)
1031 s_axi = {}
1032 for a, P_xi in kpt.P_ani.items():
1033 dO_ii = wfs.setups[a].dO_ii
1034 s_xi = np.dot(P_xi, dO_ii)
1035 s_axi[a] = s_xi
1036 wfs.pt.add(psi, s_axi, kpt.q)
1038 grad[k] -= np.tensordot(lamb.T, psi, axes=1)
1040 minstate = np.argmin(np.diagonal(lamb, offset=-n_occ).real)
1041 energy = np.diagonal(lamb, offset=-n_occ)[minstate].real
1042 norm = []
1043 for i in [minstate]:
1044 norm.append(self.dot(
1045 wfs, grad[k][i], grad[k][i], kpt, addpaw=False).item())
1046 error = sum(norm).real * Hartree ** 2 / len(norm)
1047 error_t += error
1048 energy_t += energy
1050 error_t = wfs.kd.comm.sum_scalar(error_t)
1051 energy_t = wfs.kd.comm.sum_scalar(energy_t)
1052 self.error = error_t
1054 wfs.timer.stop('Gradient unoccupied orbitals')
1056 return energy_t, grad
1058 def run_unocc(self, ham, wfs, dens, max_err, log):
1060 """
1061 Converge unoccupied orbitals
1063 :param ham:
1064 :param wfs:
1065 :param dens:
1066 :param max_err:
1067 :param log:
1068 :return:
1069 """
1071 self.need_init_odd = False
1072 self.initialize_dm(wfs, dens, ham, converge_unocc=True)
1074 while self.iters < self.maxiter_unocc:
1075 en, er = self.iterate(ham, wfs, dens, log, converge_unocc=True)
1076 log_f(self.iters, en, er, log)
1077 # it is quite difficult to converge unoccupied orbitals
1078 # with the same accuracy as occupied orbitals
1079 if er < max(max_err, 5.0e-4):
1080 log('\nUnoccupied orbitals converged after'
1081 ' {:d} iterations'.format(self.iters))
1082 break
1083 if self.iters >= self.maxiter_unocc:
1084 log('\nUnoccupied orbitals did not converge after'
1085 ' {:d} iterations'.format(self.iters))
1087 def run_inner_loop(self, ham, wfs, dens, grad_knG, niter=0):
1089 """
1090 calculate optimal orbitals among occupied subspace
1091 which minimizes SIC.
1092 """
1094 if self.iloop is None and self.outer_iloop is None:
1095 return niter, False
1097 wfs.timer.start('Inner loop')
1099 if self.printinnerloop:
1100 log = parprint
1101 else:
1102 log = None
1104 if self.iloop is not None:
1105 if self.excited_state and self.iters == 0:
1106 eks = self.update_ks_energy(ham, wfs, dens)
1107 else:
1108 etotal = ham.get_energy(
1109 0.0, wfs, kin_en_using_band=False, e_sic=self.e_sic)
1110 eks = etotal - self.e_sic
1111 if self.initial_random and self.iters == 1:
1112 small_random = True
1113 else:
1114 small_random = False
1115 self.e_sic, counter = self.iloop.run(
1116 eks, wfs, dens, log, niter, small_random=small_random,
1117 seed=self.localizationseed)
1118 self.eg_count_iloop = self.iloop.eg_count
1119 self.total_eg_count_iloop += self.iloop.eg_count
1121 if self.outer_iloop is None:
1122 if grad_knG is not None:
1123 for kpt in wfs.kpt_u:
1124 k = self.n_kps * kpt.s + kpt.q
1125 n_occ = get_n_occ(kpt)[0]
1126 grad_knG[k][:n_occ] += np.tensordot(
1127 self.iloop.U_k[k].conj(),
1128 self.iloop.odd_pot.grad[k], axes=1)
1129 wfs.timer.stop('Inner loop')
1131 ham.get_energy(0.0, wfs, kin_en_using_band=False,
1132 e_sic=self.e_sic)
1133 return counter, True
1135 for kpt in wfs.kpt_u:
1136 k = self.iloop.n_kps * kpt.s + kpt.q
1137 U = self.iloop.U_k[k]
1138 n_occ = U.shape[0]
1139 kpt.psit_nG[:n_occ] = np.tensordot(
1140 U.T, kpt.psit_nG[:n_occ], axes=1)
1141 # calc projectors
1142 wfs.pt.integrate(kpt.psit_nG, kpt.P_ani, kpt.q)
1144 self.etotal, counter = self.outer_iloop.run(
1145 wfs, dens, log, niter, ham=ham)
1146 self.eg_count_outer_iloop = self.outer_iloop.eg_count
1147 self.total_eg_count_outer_iloop += self.outer_iloop.eg_count
1148 self.e_sic = self.outer_iloop.esic
1149 for kpt in wfs.kpt_u:
1150 k = self.n_kps * kpt.s + kpt.q
1151 grad_knG[k] += np.tensordot(self.outer_iloop.U_k[k].conj(),
1152 self.outer_iloop.odd_pot.grad[k],
1153 axes=1)
1154 if self.iloop is not None:
1155 U = self.iloop.U_k[k]
1156 n_occ = U.shape[0]
1157 kpt.psit_nG[:n_occ] = \
1158 np.tensordot(U.conj(),
1159 kpt.psit_nG[:n_occ], axes=1)
1160 # calc projectors
1161 wfs.pt.integrate(kpt.psit_nG, kpt.P_ani, kpt.q)
1162 grad_knG[k][:n_occ] = \
1163 np.tensordot(U.conj(),
1164 grad_knG[k][:n_occ], axes=1)
1165 self.iloop.U_k[k] = \
1166 self.iloop.U_k[k] @ self.outer_iloop.U_k[k]
1167 self.outer_iloop.U_k[k] = np.eye(n_occ, dtype=self.dtype)
1169 wfs.timer.stop('Inner loop')
1171 ham.get_energy(0.0, wfs, kin_en_using_band=False,
1172 e_sic=self.e_sic)
1174 return counter, True
1176 def initialize_orbitals(self, wfs, ham):
1177 if self.need_init_orbs and not wfs.read_from_file_init_wfs_dm:
1178 if self.gpaw_new:
1179 def Ht(psit_nG, out, spin):
1180 return wfs.hamiltonian.apply(
1181 wfs.potential.vt_sR,
1182 wfs.potential.dedtaut_sR,
1183 wfs.ibzwfs,
1184 wfs.density.D_asii,
1185 psit_nG,
1186 out,
1187 spin)
1189 for w in wfs.ibzwfs:
1190 w.subspace_diagonalize(Ht, wfs.potential.dH)
1191 else:
1192 for kpt in wfs.kpt_u:
1193 wfs.orthonormalize(kpt)
1194 self.eigensolver.subspace_diagonalize(
1195 ham, wfs, kpt, True)
1196 wfs.gd.comm.broadcast(kpt.eps_n, 0)
1197 self.need_init_orbs = False
1198 if wfs.read_from_file_init_wfs_dm:
1199 self.initial_random = False
1201 def localize(self, wfs, dens, ham, log):
1202 if not self.need_localization:
1203 return
1204 localize_orbitals(wfs, dens, ham, log, self.localizationtype,
1205 tol=self.localization_tol,
1206 seed=self.localizationseed,
1207 func_settings=self.func_settings)
1208 self.need_localization = False
1210 def choose_optimal_orbitals(self, wfs):
1211 """
1212 choose optimal orbitals and store them in wfs.kpt_u.
1213 Optimal orbitals are those which minimize the energy
1214 functional and might not coincide with canonical orbitals
1216 :param wfs:
1217 :return:
1218 """
1219 for kpt in wfs.kpt_u:
1220 k = self.n_kps * kpt.s + kpt.q
1221 if self.iloop is not None:
1222 dim = self.iloop.U_k[k].shape[0]
1223 kpt.psit_nG[:dim] = \
1224 np.tensordot(
1225 self.iloop.U_k[k].T, kpt.psit_nG[:dim],
1226 axes=1)
1227 self.iloop.U_k[k] = np.eye(self.iloop.U_k[k].shape[0])
1228 self.iloop.Unew_k[k] = np.eye(
1229 self.iloop.Unew_k[k].shape[0])
1230 if self.outer_iloop is not None:
1231 dim = self.outer_iloop.U_k[k].shape[0]
1232 kpt.psit_nG[:dim] = \
1233 np.tensordot(
1234 self.outer_iloop.U_k[k].T,
1235 kpt.psit_nG[:dim], axes=1)
1236 self.outer_iloop.U_k[k] = np.eye(
1237 self.outer_iloop.U_k[k].shape[0])
1238 self.outer_iloop.Unew_k[k] = np.eye(
1239 self.outer_iloop.Unew_k[k].shape[0])
1240 if self.iloop is not None or \
1241 self.outer_iloop is not None:
1242 wfs.pt.integrate(kpt.psit_nG, kpt.P_ani, kpt.q)
1244 def check_assertions(self, wfs, dens):
1246 assert dens.mixer.driver.basemixerclass.name == 'no-mixing', \
1247 'Please, use: mixer={\'backend\': \'no-mixing\'}'
1248 assert wfs.bd.comm.size == 1, \
1249 'Band parallelization is not supported'
1250 if wfs.occupations.name != 'mom':
1251 errormsg = \
1252 'Please, use occupations={\'name\': \'fixed-uniform\'}'
1253 assert wfs.occupations.name == 'fixed-uniform', errormsg
1255 def check_restart(self, wfs):
1256 occ_name = getattr(wfs.occupations, 'name', None)
1257 if occ_name != 'mom':
1258 return
1260 sic_calc = 'SIC' in self.func_settings['name']
1261 if self.outer_iloop is not None:
1262 # mom restart ?
1263 astmnt = self.outer_iloop.restart
1264 # iloop not converged?
1265 bstmnt = \
1266 (self.iters + 1) % self.restartevery_iloop_notconverged == 0 \
1267 and not self.outer_iloop.converged
1268 if astmnt or bstmnt:
1269 self.choose_optimal_orbitals(wfs)
1270 if not sic_calc and self.restart_canonical:
1271 # Will restart using canonical orbitals
1272 self.need_init_orbs = True
1273 wfs.read_from_file_init_wfs_dm = False
1274 self.iters = 0
1275 self.initialized = False
1276 self.need_init_odd = True
1278 return self.outer_iloop.restart
1280 def initialize_mom_reference_orbitals(self, wfs, dens):
1281 # Reinitialize the MOM reference orbitals
1282 # after orthogonalization/localization
1283 occ_name = getattr(wfs.occupations, 'name', None)
1284 if occ_name == 'mom' and self.globaliters == 0:
1285 for kpt in wfs.kpt_u:
1286 wfs.pt.integrate(kpt.psit_nG, kpt.P_ani, kpt.q)
1287 wfs.orthonormalize()
1288 wfs.occupations.initialize_reference_orbitals()
1289 wfs.calculate_occupation_numbers(dens.fixed)
1292def log_f(niter, e_total, eig_error, log):
1293 """
1294 log function for convergence of unoccupied states.
1296 :param niter:
1297 :param e_total:
1298 :param eig_error:
1299 :param log:
1300 :return:
1301 """
1303 T = time.localtime()
1304 if niter == 1:
1305 header = ' ' \
1306 ' wfs \n' \
1307 ' time Energy:' \
1308 ' error(ev^2):'
1309 log(header)
1311 log('iter: %3d %02d:%02d:%02d ' %
1312 (niter,
1313 T[3], T[4], T[5]
1314 ), end='')
1316 log('%11.6f %11.1e' %
1317 (Hartree * e_total, eig_error), end='')
1319 log(flush=True)