Propagated (fragment) Pipek–Mezey Wannier functions in real-time time-dependent density functional theory

Due to its high accuracy relative to computational demands, density functional theory (DFT)^1,2 is the most popular electronic structure method of quantum chemistry.³ Modern quantum chemistry software implements Born–von Karman periodic boundary conditions (PBC) on unit cells or supercells of a given system to simulate bulk properties in DFT calculations. This method leads to periodic Bloch functions representing crystalline orbitals, as opposed to canonical orbitals (COs) in non-periodic calculations.⁴ In the following, we use the abbreviation CO for both crystalline and canonical orbitals depending on whether the system is periodic or not. The COs are not restricted in their spatial extent in any way and, as a result, often delocalized over the entire system, making them hard to interpret from chemical intuition and in terms of regional phenomena. The established approach to mitigate these issues in post-processing is to perform a localization procedure.^5,6 Localization refers to optimization of a unitary transformation with respect to a certain localization measure, which is then applied to the COs to obtain a set of localized orbitals (LOs) without changing any observables. Under PBC, this results in a transformation into maximally localized Wannier functions (MLWFs).^7–10 The most important of such spread-based functionals include the Foster–Boys (FB) criterion¹¹ and its Berry phase analogue for periodic systems,^12,13 leading to Foster–Boys Wannier functions (FBWFs). This approach is highly successful in making the orbitals more spatially confined (according to the FB criterion) than the COs but has the downside of commonly mixing $σ$- and $π$-type orbitals. Such mixing results in asymmetric hybrid bonds as shown in Fig. 1, occasionally referred to as “banana” bonds. There, on the example of $CO2$⁠, we see that the delocalized canonical highest occupied molecular orbital (HOMO) spans all three atoms, which is not aligned with chemical intuition expecting two separate double bonds. In contrast, the corresponding FBWF⁷ is clearly localized between two atomic centers, but the expected $π$-symmetry is lost due to admixture of the $σ$-bond.

This mixing is often inconsequential to the problem at hand, e.g., when expectation values are summed over many eigenstates. In cases where mixing is undesirable, a different localization criterion is needed. The established Pipek–Mezey (PM) criterion,¹⁴ which is based on a partial charge functional, preserves the $σ$/$π$ symmetry of the LOs. However, there exists ambiguity in the choice of the partial charge estimate: The original formulation used Mulliken charges.¹⁵ Since those suffer from having no proper basis set limit, Lehtola and Jónsson generalized the criterion to other partial charge estimates, including well-defined real-space methods.¹⁶ While this gives us a great deal of flexibility to choose whatever partial charge estimate is most suited to the problem at hand, they did find that all estimates tested yielded similar results. Thus, they settled on a flexible, generalized real-space method.¹⁰ The robustness of this PM localization scheme has recently been demonstrated by Clement et al.,¹⁷ which made use of an atomic charge estimate by projection into a minimal basis set and also introduced the idea of the Canonicalize Phase then Randomize (CPR) protocol for generating the initial guess of the localization in periodic systems. In this article, we will use the acronyms FBWF and PMWF for MLWFs obtained using the FB/Berghold and PM spread functionals, respectively. By design, subsystem density embedding methods also keep the COs local to the respective subsystem, which can be exploited to investigate local properties of larger and periodic systems.^18,19

The theoretical foundations of DFT are generally limited to the electronic ground state. Thus, modifications are required to capture excited state dynamics, which give rise to some of the most interesting material properties, such as photochemistry, an increasingly important field of chemistry, e.g., in the pursuit of green energy. Time-dependent density functional theory (TDDFT)^20–22 can be employed in these cases for detailed insights into excited state dynamics at a low computational cost. The two major approaches to TDDFT are a perturbative treatment,^20,23 which operates in the frequency domain, and the real-time approach, which records the response to a given perturbation in the time domain.²⁴ Real-time TDDFT has gained popularity over the past years, mainly due to its favorable scaling at higher density of states, and the ability of predicting natural linewidths,^25–28 with our group successfully applying it to theoretical spectroscopy²⁹ in the fields of electronic circular dichroism,^30,31 (resonance) Raman,^32,33 and (resonance) Raman optical activity.³² There are also variational DFT-based approaches to excited state dynamics, such as the $Δ$ self-consistent field method.^34–41 

Application of localization procedures in real-time TDDFT, however, requires an adapted formulation with appropriate optimization algorithms. In this article, we present a PMWF algorithm for propagated orbitals and its implementation in the massively parallel quantum chemistry package CP2K,⁴² based on the Cardoso–Souloumiac algorithm⁴³ we previously employed for propagated FBWFs.^19,44,45 We discuss its effectiveness in $σ$/$π$ separation and computational scaling on the example of stacked slipped carbon dioxide oligomers, polymeric carbon chains, boron and nitrogen (B/N) co-doped graphene, and nitrogen-vacancy doped (NV) diamond. All systems under investigation are subject to PBC, considering the $Γ$ point of the Brillouin zone.

This theoretical discussion adheres to the following conventions:

• Bold lowercase letters signify vectors, while bold capital letters signify matrices.

• Latin subscripts run over COs/LOs, and Greek subscripts run over atomic orbitals (AOs).

• Time- and location-dependence of CO/LO coefficients and wavefunctions is implied but might be made explicit for clarity.

• The Born–Oppenheimer approximation is employed, and atomic positions are fixed in all examples.⁴⁶ 

• Hartree atomic units⁴⁷ are used throughout.

• Repeated indices imply summation following the Einstein summation convention unless explicit sums are given.⁴⁸ 

The general PMWF scheme with an CO-based partial charge estimate scales approximately $nA×nCO2$ since $nA$ matrices of $nCO$ rows and columns need to be diagonalized. As $nCO$ grows linearly with $nA$⁠, the overall scaling of the PMWF localization scheme is approximately $O(n3)$⁠. This is in contrast to spatial schemes, where there are only one (FB) or two (Berghold) $nCO×nCO$ matrices per Cartesian direction, resulting in approximate $O(n2)$ scaling. To counteract this, there have been efforts to compress the $nA$ factor in large systems for computational gain.^57,58

Introduced in Ref. 58, the sequential variant of F-PMWF (sF-PMWF) uses a double loop to forgo any knowledge of the full set of $nMO$ states. Here, we will present an algorithmic overview of this procedure.

The objective of this algorithm is to find $nF$ states that are localized on a selected molecular fragment. To do so, we firstly identify a “core” region of $nF$ states with the largest values for the locality measure $Li$⁠. Any remaining states become part of the “rest” space and are ordered by their locality. The “rest” space is then divided into $nb$ blocks, of which the sizes are given by an arbitrary parameter $nr$ that can range from 1 to the total size of the rest space $nMO−nF$⁠. If the size of the rest space is not divisible by $nr$⁠, the size of the last block will be the modulo of that division. After this initialization step, the outer loop begins, by selecting one of the blocks in the rest space, selected by taking the modulo of the iteration counter and $nb$⁠, ensuring that the rest space is rotated into the actual localization procedure one by one. The inner loop performs a TD-PMWF optimization as outlined before, with only the states enumerated by the core space and the current block being used in the loop over index combinations. After the inner loop terminates, the core and rest spaces are recomputed, the next block of the rest space is rotated into the inner loop, and the next iteration is performed. The outer loop terminates when the spread functional has converged, which is checked after each update of the core and rest spaces. A pseudocode rendition of this TD-sF-PMWF algorithm is given in Algorithm 2. This algorithm opens up a new parallelization opportunity in that one can rotate a different block into the inner loop on each processor and only choose the one with the biggest increase in localization to actually apply to the common set of matrices. It remains to be seen if this strategy performs better or worse than parallelizing the inner loop as implemented before.¹⁹ 

Since the time-dependent MOs are complex-valued, we color-mapped the following heatmaps accordingly: The argument of a complex number $z∈C$ $arg(z)∈[0,2π)$⁠, with $arg(1)=0$⁠), normalized to $[0,1]$⁠, is mapped to the hue, while its magnitude $|z|$⁠, normalized to $[0,1]$⁠, is mapped to the saturation in an HSV color-mapping⁶⁹ where the value parameter V is kept constant at $V=1$⁠. The resulting colormap is shown in Fig. 2. A three-dimensional approach to complex-valued MOs was introduced in Ref. 70.

The following examples show four distinct test cases: first, we demonstrate the method on stacked $CO2$ molecules, which demonstrate the method’s ability to localize states from intermolecular bonds to intramolecular $π$-bonds. We also test on linear acetylenic carbon (LAC), to assess the method’s ability to localize states from a totally delocalized $π$ system to a well-defined $C=C$ double bond. Third, we examine its ability to localize a heteroatomic bond in a totally delocalized $π$ system on the example of B/N-co-doped graphene. Finally, we examine its abilities in a 3D-periodic system on the example of NV diamond. The orbitals shown are the result of localization after $1ℏEh−1$ of propagation time, corresponding to ten steps.

The stacks of $CO2$ molecules shown here are extrapolated from the global minimum geometry of the $CO2$ dimer.⁷¹ Therefore, an attractive interaction is expected, which is symmetrically allowed for either $σ$-type interactions of the in-plane $π$-orbitals oriented head-on or $π$-type interactions of the $π$-orbitals perpendicular to the plane. Here, we investigate two orbitals: the HOMO, which is formed by the $σ$-type interactions of the in-plane $π$-orbitals, and the highest occupied $π$-orbital perpendicular to the molecular plane. Figure 3 shows the evolution of the former as a PMWF over the first 100 time steps of 0.1 ℏ/$Eh$ (⁠$≈$2.42 a.s.) each, while Fig. 4 shows a top-down cross section through the bond axes of the stack, to present a sideways view of the expected $π$-bonds and a sideways view to present a cross section of the perpendicular $π$-bond at the end of the propagation. Without any localization method applied, the stacked $CO2$ molecules display a canonical HOMO of predominantly three-center character as in Fig. 1. In addition, the orbital is significantly delocalized over all molecules of the system, which renders it useless to discuss any local effects. After applying a TD-FBWF localization, the Wannier function is localized entirely on a single molecule. However, like in the isolated molecule, the expectedly $π$-symmetric Wannier functions are significantly mixed with the corresponding $C=O$ $σ$-bonds, leading to a loss of the nodal plane through the bond. The result of this loss in symmetry is a banana bond in both orientations investigated here. When applying a TD-PMWF optimization instead, we successfully recover the expected nodal plane for a much improved $π$-character of the Wannier function.

LAC is an artificial allotrope of carbon consisting of linear chains of repeated $(–C≡C–)n$ units.⁷² While it is of considerable interest to nanotechnology, it should be noted that it is likely to cross-link exothermically when two side chains approach each other, making it difficult to employ in condensed phase applications. Nevertheless, with it being the simplest carbon polymer possible, it can serve here as a toy system for observing a $π$-system delocalized over one full dimension of the simulation box. PMWFs of LAC have previously been reported in Ref. 10. We have chosen a chain of eight carbon atoms as the repeating unit. The expected HOMO in that case is any of the eight degenerate $π$-bonds since they are symmetrically equivalent. Thus, Fig. 5 shows a sideways cross section through the $π$-bond of the arbitrarily chosen state. The COs that we obtained directly after the SCF procedure, indeed, show eight orbitals of equal energy distributed uniformly to all nuclei. After maximal localization, this picture changes drastically. While the FBWF is now confined to a single $C≡C$ bond, it has been mixed with the bond’s $σ$-orbital, resulting in an almost total loss of $π$-character. Only a slight skewedness of the orbital with regard to the rotational symmetry axis remains from the $π$-type bond. In contrast, the PMWF retains a perfect nodal plane and shows a predominantly two-center bond, with a slight antibonding influence on the neighboring $C– C$ pairs.

In the research effort to further mature graphene-based microelectronics and nanotechnology, doping techniques are essential to fine-tune its properties at the molecular level. Here, we investigate a 2D-periodic graphene sheet doped with a pair of boron and nitrogen, which has recently been shown to behave as a peroxidase mimic⁷³ and allows for finely tuned redox behavior when adsorbed on metal surfaces,⁷⁴ depending on the doping configuration. The graphene sheet under investigation in this article has been doped with a directly connected $B– N$ pair, with the focus being the expected $π$-orbital between them. As in the LAC example, Fig. 6 shows a sideways view of the system, to provide a cross section view of the expected $π$-bond. We identified the corresponding Wannier function by the position of its Wannier center. After the SCF optimization, the CO displays delocalization over the whole simulation box, with no discernible $π$-symmetry. The corresponding FBWF becomes strongly confined to the $B–N$ system but is still missing any noticeable $π$-symmetry. Finally, the PMWF procedure is able to both confine the orbital to the bonding pair and reproduce the expected $π$-symmetry.

In contrast to graphene, where the electrons are totally delocalized over the 2D sheet, the diamond allotrope of carbon only possesses localized $σ$-bonds. This property is immensely useful to frustrate electronic systems, as the resulting well-localized electron spin interacts with a wide range of electromagnetic radiation.⁷⁵ Among other applications, the NV center has been used as a qubit.^76,77 Here, we investigate the contribution of the lone pair on the nitrogen center to the vacancy. Therefore, we identified the corresponding Wannier function by its Wannier center and present in Fig. 7 a cross section perpendicular to the longest cell axis in the z-direction, through the NV cavity. In the CO picture, the orbital is already comparatively well-localized within the vacant site and the neighboring carbon and nitrogen atoms, with stronger localization toward the nitrogen site, as expected from their relative electronegativities. Both localization procedures fully agree with each other in that this particular contribution to the vacancy center is fully localized to the nitrogen site, forming an orbital shape strongly resembling an $sp3$ orbital from hybridization theory. This example demonstrates how FBWF and PMWF procedures can both arrive at the same solution in naturally well-localized systems, a finding also noted in Ref. 10.

We used the slipped parallel $CO2$ stacks to carry out a scaling test between our TD-FBWF, TD-PMWF, and TD-F-PMWF implementations. With our computational setup as described, we were able to localize up to 32 stacked $CO2$ molecules before the slowest method at that size (TD-PMWF) failed to finish the localization procedure within one hour. With our pseudopotential-based methodology, each molecule contributes eight states, as well as three partial charge matrices to the localization. Figure 8 shows the average walltime required to perform a single real-time TDDFT step including the localization procedure for each combination of algorithms and system sizes. As expected, we observe an initial parity of computational demand between FBWF and PMWF (with a slight edge for PMWF). At more than eight molecules, the PMWF algorithm shows its weakness of the ever-growing amount of partial charge matrices to consider. This makes for much worse scaling behavior than the FBWF implementation, where the amount of matrices to diagonalize remains constant. In contrast, the fragment-based PMWF method (in which we selected one of the $CO2$ molecules as the fragment) suffers some overhead at low molecule counts but is able to retain a very close scaling behavior to FBWF, improving considerably over basic PMWF. Thus, the applicability of PMWF to larger systems is vastly broadened to systems accessible via FBWF using suitable fragments in F-PMWF.

To investigate the impact of $σ$/$π$ mixing at the application level, we performed a real-time TDDFT simulation with subsystem analysis of $CO2$⁠. The parameters of the simulation were identical to the previous calculations except for the number of propagation steps, which was set to 10 000 for a total propagation time of 1000 ℏ/$Eh$⁠. The subsystem analysis followed the process outlined in Ref. 19. In the TD-PMWF case, we replaced the TD-FBWF procedure with the TD-PMWF subroutine introduced in this work. The spectrum (Fig. 9) shows the signals of the $σ$- and $π$-bonds between the p-orbitals of one $C=O$ bond localized with TD-PMWF and the mixed $σ$/$π$ bond localized with TD-FBWF. We observe that the effect of $σ$/$π$ mixing is significant: The principal difference between the well-separated $σ$ and $π$ orbitals is that the $σ$ orbital absorbs in all three peaks of the depicted range, while the $π$ orbital only absorbs at 86 nm. In contrast, the FBWF orbital displays mixed properties: it does not absorb at 95 nm, just like the well-separated $π$ orbital. At 86 nm, it absorbs with an intermediate intensity between the $σ$ and $π$ orbitals. Then, at 83 nm, it absorbs like the $σ$ orbital but at only half of the intensity. In cases like these, having the choice between different localization procedures can be a significant advantage.

In this article, we have adapted the PMWF localization procedure to real-time TDDFT and presented our implementation into the massively parallel quantum chemistry package CP2K. We have shown its efficacy in localizing time-dependent electronic states in cases where traditional FBWFs fall short, e.g., by $σ$/$π$-mixing, and produce qualitatively equivalent results where they do not. The method performed well in this regard as well as computational efficiency for all tested systems, including stacked $CO2$ oligomers, 1D carbon polymers, B/N-co-doped graphene, and NV diamond. Inclusion of the sequential fragment variant proved vital to make the TD-PMWF approach feasible on large, periodic systems. We especially see potential in their application to the resolution of real-time TDDFT-derived theoretical spectroscopy by individual PMWFs on user-selected atomic fragments, which should be explored on more complex systems in further research efforts.

This work was supported by the University of Zurich and the Swiss National Science Foundation (Grant No. 200021-197207). We thank the Swiss National Supercomputing Centre for computing resources (Project No. 1036).

The authors have no conflicts to disclose.

Lukas Schreder: Data curation (lead); Formal analysis (lead); Investigation (lead); Validation (lead); Visualization (lead); Writing – original draft (lead). Sandra Luber: Funding acquisition (lead); Project administration (lead); Resources (lead); Supervision (lead); Writing – review & editing (lead).

The data that support the findings of this study are available from the corresponding author upon reasonable request.