Octopus, a computational framework for exploring light-driven phenomena and quantum dynamics in extended and finite systems
Nicolas Tancogne-Dejean, Micael J. T. Oliveira, Xavier Andrade, Heiko Appel, Carlos H. Borca, Guillaume Le Breton, Florian Buchholz, Alberto Castro, Stefano Corni, Alfredo A. Correa, Umberto De Giovannini, Alain Delgado, Florian G. Eich, Johannes Flick, Gabriel Gil, Adrián Gomez, Nicole Helbig, Hannes Hübener, René Jestädt, Joaquim Jornet-Somoza, Ask H. Larsen, Irina V. Lebedeva, Martin Lüders, Miguel A. L. Marques, Sebastian T. Ohlmann, Silvio Pipolo, Markus Rampp, Carlo A. Rozzi, David A. Strubbe, Shunsuke A. Sato, Christian Schäfer, Iris Theophilou, Alicia Welden, Angel Rubio
OOctopus, a computational framework for exploring light-driven phenomena andquantum dynamics in extended and finite systems
Nicolas Tancogne-Dejean, a) Micael J. T. Oliveira, b) Xavier Andrade, Heiko Appel, Carlos H. Borca, Guillaume Le Breton, Florian Buchholz, Alberto Castro,
4, 5
StefanoCorni,
6, 7
Alfredo A. Correa, Umberto De Giovannini, Alain Delgado, Florian G. Eich, Johannes Flick,
9, 10
Gabriel Gil,
6, 11
Adri´an Gomez, Nicole Helbig, Hannes H¨ubener, Ren´e Jest¨adt, Joaquim Jornet-Somoza, Ask H. Larsen, Irina V. Lebedeva, MartinL¨uders, Miguel A. L. Marques, Sebastian T. Ohlmann, Silvio Pipolo, MarkusRampp, Carlo A. Rozzi, David A. Strubbe, Shunsuke A. Sato,
1, 18
ChristianSch¨afer, Iris Theophilou, Alicia Welden, and Angel Rubio
1, 13, 10, c)1)
Max Planck Institute for the Structure and Dynamics of Matter,Luruper Chaussee 149, D-22761 Hamburg, Germany Quantum Simulations Group, Lawrence Livermore National Laboratory,Livermore, California 94551, USA D´epartement de Physique, ´Ecole Normale Sup´erieure de Lyon, 46 All´ee d’Italie,Lyon Cedex 07, France Institute for Biocomputation and Physics of Complex Systems,University of Zaragoza, Calle Mariano Esquillor, 50018 Zaragoza,Spain ARAID Foundation, Avda. de Ranillas 1-D, 50018 Zaragoza,Spain Dipartimento di Scienze Chimiche, Universit`a degli studi di Padova,via F. Marzolo 1, 35131 Padova, Italy CNR – Istituto Nanoscienze, via Campi 213a, 41125 Modena,Italy Xanadu, 777 Bay Street, Toronto, Ontario, M5G 2C8,Canada John A. Paulson School of Engineering and Applied Sciences, Harvard University,Cambridge, MA 02138, USA
Center for Computational Quantum Physics, Flatiron Institute, 162 5th Avenue,New York, NY 10010.
Instituto de Cibern´etica, Matem´atica y F´ısica, Calle E 309, 10400 La Habana, a r X i v : . [ phy s i c s . c o m p - ph ] D ec uba Nanomat/Qmat/CESAM and ETSF, Universit´e de Li`ege, B-4000 Sart-Tilman,Belgium
Nano-Bio Spectroscopy Group and ETSF, Universidad del Pa´ıs Vasco,20018 San Sebasti´an, Spain
Institut f¨ur Physik, Martin-Luther-Universit¨at Halle-Wittenberg,06120 Halle (Saale), Germany
Max Planck Computing and Data Facility, Gießenbachstraße 2, 85741 Garching,Germany
Universit de Lille, CNRS, Centrale Lille, ENSCL, Universit d ArtoisUMR 8181 UCCS Unit de Catalyse et Chimie du Solide, F-59000, Lille,France
Department of Physics, School of Natural Sciences, University of California,Merced, CA, 95343, USA
Center for Computational Sciences, University of Tsukuba, Tsukuba 305-8577,Japan ab initio computer simulation tool that enables a reliable and accurate simulationof light-induced changes in the physical and chemical properties of complex systemsis of utmost importance. The first principles real-space-based Octopus project wasborn with that idea in mind, providing an unique framework allowing to describenon-equilibrium phenomena in molecular complexes, low dimensional materials, andextended systems by accounting for electronic, ionic, and photon quantum mechani-cal effects within a generalized time-dependent density functional theory framework.The present article aims to present the new features that have been implementedover the last few years, including technical developments related to performance andmassive parallelism. We also describe the major theoretical developments to ad-dress ultrafast light-driven processes, like the new theoretical framework of quantumelectrodynamics density-functional formalism (QEDFT) for the description of novellight-matter hybrid states. Those advances, and other being released soon as part ofthe Octopus package, will enable the scientific community to simulate and character-ize spatial and time-resolved spectroscopies, ultrafast phenomena in molecules andmaterials, and new emergent states of matter (QED-materials). a) Electronic mail: [email protected] b) Electronic mail: [email protected] c) Electronic mail: [email protected] . INTRODUCTION It is a general challenge in the electronic structure community to develop accurate andefficient methods of modeling materials of ever increasing complexity in order to predicttheir properties. In this respect, time-dependent density functional theory (TDDFT) andrelated methods have become a natural choice for modeling materials and complex systemsin and out of their equilibrium.During the last years, novel directions of research emerged in the field of chemistry,physics, and material science that required the development of novel simulation tools neededto face the new challenges posed by those experimental advances and emerging phenomena.Among these new fields of research, one could cite some examples: the strong couplingbetween light and matter (including materials embedded in cavities), new states of matter(hidden phases and topological solids and molecules), and the strong field dynamics inperiodic systems. Whereas the strong-field dynamics of atoms and molecules is now wellunderstood, the strong-field dynamics in solids is an active field of research. Real-timeTDDFT represents a natural tool to study highly non-linear phenomena in solids andlow-dimensional materials and the development of efficient numerical methods to performreal-time TDDFT in such periodic or semi-periodic systems is crucial to explore this newphenomena (including the quantum nature of light and phonons). Indeed, it allows todescribe highly nonlinear processes without having to resort to the perturbation theory ontop of equilibrium DFT calculations.Recent years have seen tremendous experimental progress in the field of strong light-matter interactions, where the strong coupling of light to chemical systems, quantummaterials, or nanoplasmonic systems, among others, has been demonstrated. In this regime,light and matter meet on the same footing and the electron-photon interaction has to beexplicitly considered.
A theoretically novel approach accelerating this field is quantum-electrodynamical density-functional theory (QEDFT), which complements TDDFTwith the photonic degrees of freedom and provides reliable and predictive simulations in thisemerging field of research.In this paper we explore the recent advances in the Octopus code project.
We focusour attention on the recently added features, and particularly on the ones that have notbeen described in the previous papers.
These new features include the implementation4f new levels of self-consistent and microscopic couplings of light and matter, the treatmentof solvent effects through the polarizable continuum model, the implementation of variousmethods to treat van der Waals interactions, new methods to calculate magnons, conductiv-ities, and photoelectron spectroscopy from real-time TDDFT, and the calculation of orbitalmagneto-optical responses. Advances in numerical algorithms and methods, such as newpropagators and the use of iterative eigensolvers in the context of reduced density matrixfunctional theory, are also discussed. Finally, recent improvements in the treatment of peri-odic systems, as well as more technical code improvements, are also presented. For a moredetailed explanation about how to use these newly introduced features in practice, we referthe readers to the Octopus webpage, as new tutorials and examples are regularly beingadded to it.This paper is organized as follow. First, we present in Sec. II to Sec. XII new implemen-tations of physical theories and algorithms that allow us to deal with new non-equilibriumphenomena in materials and nanostructures. This is followed by a second set of sections,from Sec. XIII to Sec. XVI, dealing with technical developments that are fundamental in im-proving the code performance and the algorithm’s stability. Finally, we draw our conclusionsin Sec. XVII. Unless otherwise stated, atomic units are used throughout the paper. II. COUPLED MAXWELL-KOHN-SHAM EQUATIONS
In most cases when light-matter interactions are considered, a decoupling of light andmatter is performed at the outset. Either the electromagnetic fields are prescribed and thenproperties of the matter subsystem are determined, as frequently found in, e.g., quantumchemistry or solid-state physics, or the properties of matter are prescribed and then theproperties of the photon subsystem are determined, as done in, e.g., quantum optics orphotonics.
The microscopic interaction of light and matter in Octopus has followed this decouplingstrategy and has been treated so far only in forward-coupling direction. In this approxima-tion an external classical laser pulse or kick is prescribed and the response of the system iscomputed by evolving the Kohn-Sham orbitals. The back-action of the matter subsystemon the laser pulse, and the subsequent effect of this modified pulse on matter, and so on, isignored. This forward-coupling approximation is highly accurate when the generated total5urrent in the system is comparably small, such as in atoms or smaller molecules. This hasbeen exploited in Octopus and many different types of spectroscopies have been computedsuccessfully in the past using this approach.In contrast, in classical electromagnetic modeling the opposite view is taken. Here thematerial properties are prescribed and the resulting electromagnetic fields are computed. Inpractice, the material properties are routinely approximated by a local continuum model orthe dielectric function of the system, such as a Debye, Lorentz, or Drude model and thenMaxwell’s equations are solved for linear dielectric media arranged in appropriate geome-tries. It is clear that both perspectives, a focus on matter dynamics alone or a focus on elec-tromagnetic field dynamics alone, break down when the total induced currents becomelarge and when electromagnetic near-field effects on the scale of the material system arenot negligible anymore. Prime examples for such cases are nanoplasmonic systems, surfaceplasmon-polaritons, or tip-enhanced spectroscopies. In these cases the back-action of thematerial response on the system itself has to be taken into account leading to screeningand retardation effects. The proper theoretical framework which encompasses all these ef-fects is quantum electrodynamics. Starting with a generalized Pauli-Fierz field theory forthe combined system of electrons, nuclei, and photons, we have recently derived differentlevels of self-consistent and microscopic couplings of light and matter. In the classical limitthis results in a coupled set of Ehrenfest-Maxwell-Pauli-Kohn-Sham equations. To im-plement these equations, we have added a Maxwell solver to the Octopus code which wecouple self-consistently to the dynamics of the electrons and nuclei. In the following, webriefly summarize the basic ingredients for this implementation and show an example ofself-consistent light-matter interactions for a nanoplasmonic system. Further details of theimplementation and nano-optical applications can be found in Ref. 21.Since over the years Octopus has been optimized heavily to solve time-dependentSchr¨odinger and Kohn-Sham equations, we have exploited the fact that Maxwell’s equa-tions can be formulated in Schr¨odinger form to benefit from the efficient time-evolutionin the code. This reformulation is based on the Riemann-Silberstein vector, which is acombination of the electric E ( r , t ) and magnetic field B ( r , t ) F ± ( r , t ) = (cid:114) (cid:15) E ( r , t ) ± i (cid:114) µ B ( r , t ) . (1)6he sign of the imaginary part of the Riemann-Silberstein vector corresponds to differenthelicities. The reformulation of Maxwell’s equations in Schr¨odinger form is purely algebraicand starts out with the microscopic Maxwell’s equations ∇ · E = ρ(cid:15) , ∇ · B = 0 , (2) ∇ × B = 1 c ∂ E ∂t + µ J , ∇ × E = − ∂ B ∂t , (3)where E and B are the classical electric and magnetic fields, ρ and J are the charge andcurrent densities, (cid:15) and µ are the vacuum permittivity and permeability, and c = ( (cid:15) µ ) − / is the speed of light. Using the Riemann-Silberstein vector, the electric and magnetic Gausslaws may now be combined in real and imaginary part ∇ · F = 1 √ (cid:15) ρ (4)and, likewise, the Faraday and Ampere law can be combined into one evolution equation forthe Riemann-Silberstein vector i (cid:126) ∂ F ∂t = c (cid:18) S · (cid:126) i ∇ (cid:19) F − i (cid:126) √ (cid:15) J . (5)Here S = ( S x , S y , S z ) denotes a vector of spin one matrices S x = − i i , S y = i − i , S z = − i i , (6)which are analogous to the Pauli matrices and show the spin-one character of the pho-ton. Having cast Maxwell’s equations as an inhomogeneous Schr¨odinger equation, it is nowstraightforward to use the time-evolution algorithms in Octopus to time-evolve the Riemann-Silberstein vector. The only difference to the matter propagation is that we are now dealingwith the “Maxwell Hamiltonian” H EM = c (cid:18) S · (cid:126) i ∇ (cid:19) , (7)which acts on six orbitals of the Riemann-Silberstein vector corresponding to the threecomponents of the electric and magnetic field vectors. As in the matter case in Octopus, thediscretization of the gradient in the Maxwell Hamiltonian is performed with finite-difference7tencils and the domain parallelization of Octopus can be used seamlessly for the Maxwellcase as well. While also finite-difference discretizations are used for the Maxwell solver inOctopus, the difference to finite-difference time-domain (FDTD) codes based on the Yeealgorithm is that we employ not two shifted grids for the electric and magnetic fields, butrather a single grid for the Riemann-Silberstein vector. This simplifies the coupling to matterand allows us to use higher-order finite-difference discretizations for the gradient. Sincethe spatial discretization is connected to the temporal discretization through the Courantcondition, this in turn allows to take larger time-steps and, from our experience, a unifiedgrid also improves the stability compared to FDTD.Instead of using the constitutive relations, we couple Maxwell’s equations directly tothe microscopic current density of the matter subsystem, consisting of the usual paramag-netic current term, the diamagnetic current term, and the magnetization current term. Forthe coupling of the electromagnetic fields to the matter subsystem, we are relying on thePower-Zienau-Woolley transformation, which leads to a multipole expansion. We haveimplemented the first two orders of this expansion: the dipole approximation in lowest or-der and electrical-quadrupole and magnetic-dipole coupling in the next order. In addition,we are currently working on implementing the full minimal coupling with a full positiondependence of the vector potential.The time evolution of the Kohn-Sham orbitals and the Maxwell fields is performed side-by-side and the two subsystems are coupled self-consistently in each time-step, as illustratedin Fig. 1. To propagate different subsystems with different Hamiltonians and different sets oforbitals simultaneously, we have implemented a multi-system framework in Octopus (moredetails about this can be found in Subsection XVI C).Similar to the matter propagation, also in the Maxwell case outgoing waves that reachthe box boundary of the simulation box have to be absorbed to avoid artificial reflectionsand backscattering at the boundaries. Our first attempt was to use also the mask functionsthat are used for the matter propagation in Octopus. However, in the electromagnetic casethe mask absorption of outgoing waves turned out to be not efficient enough so that weimplemented a perfectly matched layer (PML) for the Maxwell propagation.When considering incoming electromagnetic fields with optical wavelengths, the couplingto atomistic or nano-scale systems leads to a multi-scale problem. The optical wavelengthof the radiation is in this case much larger than the de Broglie wavelength of the matter.8 nly forward coupling Self-consistent forward and backward coupling Figure 1. The figure on the left-hand side illustrates the standard forward-coupling approximation:the electromagnetic fields (in blue) propagate freely and only influence the propagation of thematter (in red). The back reaction of the matter currents on the electromagnetic fields is neglected.In the figure on the right-hand side we illustrate a fully self-consistent predictor-corrector scheme fora coupled Maxwell-Pauli-Kohn-Sham time-stepping. As before, the electromagnetic fields influencethe propagation of the matter (forward coupling). However, in addition, here the currents fromthe matter propagation also influence the propagation of the electromagnetic fields (backwardcoupling). A given time-step for the matter wavefunctions and the electromagnetic fields is repeateduntil self-consistency is found (self-consistent forward-backward coupling).
Likewise, the electromagnetic waves are traveling with the speed of light, which requiressub-attosecond time-steps. We have therefore implemented different multi-scale couplingsin space and time. For example, the Maxwell simulation box can be on the same scale as thematter box. In this case the electromagnetic waves are represented as incoming analyticaltime-dependent boundary conditions and are propagated numerically inside the simulationbox. Alternatively, the Maxwell simulation box can be much larger than the matter boxto fully encompass laser pulses with optical wavelengths. In this case prolongations andinterpolations have to be used similar to multi-grid methods. Since the electronic and nuclearmotion is much slower than the time-evolution of the electromagnetic waves, we also haveimplemented a multi-scale approach for the real-time propagation. The Riemann-Silbersteinvector is propagated with frozen electronic current from the last point of interaction for manyintermediate time-steps before a coupling to the matter subsystem takes place. The numberof intermediate steps is a convergence parameter and depends on the physical situation athand.Since we now include the description of classical electromagnetic fields explicitly in our9eal-time simulations, we have directly access to the outgoing electromagnetic radiation.This allows to define electromagnetic detectors at the box boundaries which accumulatethe outgoing electromagnetic waves. We have implemented such electromagnetic detectorsin Octopus and this allows to run simulations in close analogy with experiments and todirectly observe the outgoing radiation. For example, it is then no longer needed to Fouriertransform the time signal of the matter dipole to get optical spectra, but we rather haveaccess to the spectrum directly on the Maxwell grid.As an example of a coupled Ehrenfest-Maxwell-Pauli-Kohn-Sham propagation with Oc-topus we have selected a nano-optical application. We consider in this example two almostspherical sodium nanoparticles with 297 sodium atoms each, which are arranged in a dimerconfiguration. This system is excited with an incoming laser pulse which excites either theinternal dipole or quadrupole plasmon motion of the dimer. In Fig. 2 we show the resultingelectromagnetic field enhancements for different levels of light-matter coupling. In panels a)and b) we show the temporal profile of the incoming laser pulse and the resulting currentdensity at the center point between the two nanoparticles. The field enhancement can beseen in panels c)-e). Including a self-consistent back reaction in the light-matter coupling,the field enhancement is reduced at the center point between the two nanoparticles, as canbe seen in panel c), while far away from the dimer, as shown in panels d) and e), the fieldenhancement is larger than in the forward coupled case. Furthermore, frequency shifts canalso be observed, which are more pronounced in the near-field. We have found that thefield enhancement is also sensitive to the coupling terms of the multipole expansion whichare included and that the quantitative difference of switching from LDA to PBE for theexchange-correlation functional is in this case smaller than including the back reaction inthe light-matter coupling, cf. Ref. 21.To conclude, with our new efficient implementation for coupled Ehrenfest-Maxwell-Pauli-Kohn-Sham equations in Octopus we have now a very versatile tool at hand which allows tocompute fully self-consistent forward-backward light-matter coupling in real-time and real-space for a vast set of applications and in close analogy to experiments. As the coupling toMaxwell’s equations of motion represents the classical limit of the light-matter interaction,this development leads to the classical limit of QEDFT.10 igure 2. Electric field values and current density in the z -direction for a dimer of sodiumnanoparticles. The centers of the two nanoparticles are located along the z -axis and thedistance between the two effective spheres of the nanoparticles is 0.5 nm. The electric fieldvalues are calculated at two different points: the center point r cp between the two nanoparti-cles, located at the origin, and an off-center point r ocpx located along the x -axis at a distanceof 1.957 nm from r cp . The first panel a ) illustrates the incident cosinusoidal laser pulsewith frequency ω = 3 .
05 eV (0 .
112 a.u.), λ = 406 . .
84 a.u.), and amplitude of E z = 5 . × V/m (10 − a.u.), which drives the system. The second panel b ) displaysthe current density at r cp . In panels c ) − e ), we show the electric field enhancement at r cp ,the electric field enhancement at r ocpx , and the average of the electric field over the detectorsurface close to the box boundary, respectively. The curve in bright gray in panel c ) has beenadded to simplify the comparison and is identical to the laser pulse in panel a ). The period T = 1 .
36 fs corresponding to the laser frequency ω is indicated with grey vertical lines. II. STRONG ELECTRON-PHOTON INTERACTIONS IN REAL SPACE:QUANTUM-ELECTRODYNAMICAL DENSITY-FUNCTIONAL THEORY(QEDFT)
The nascent field of strong light-matter interaction expanded over the past decades fromsmall atomic structures to chemistry and material science. This development necessitatespredictive first-principles methods capable to describe light and matter on the same footing.We have introduced for the first time a time-dependent density functional theory for quantumelectrodynamics (QEDFT) to treat ab-initio weak and strong light-matter interactions,and its applications to chemistry and materials, that provides a unique framework to explore,predict, and control new states of matter out of equilibrium. This generalization of thetime-dependent density-functional method allows us, for the first time, to explore the effectsof dressing electronic states with photons, while retaining the electronic properties of realmaterials.A general non-relativistic Hamiltonian for light-matter systems treating N interactingelectrons coupled to N p photonic modes in the case of the so-called length-gauge, and afteremploying the long-wavelength (dipole) approximation, reads as follows:ˆ H = N (cid:88) k =1 (cid:2) − ∂ r k + v ( r k , t ) (cid:3) + (cid:88) k (cid:54) = l w ( r k , r l ) + 12 N p (cid:88) α =1 − ∂ ∂q α + (cid:32) ω α q α − λ α · N (cid:88) k =1 r k (cid:33) + 2 j ( α )ext ( t ) ω α q α . (8)Here the first two terms on the right-hand side correspond to the usual electronic many-bodyHamiltonian, while the last term describes the photon mode, which is characterized for eachphoton mode α by its elongation q α , frequency ω α , and electron-photon coupling strengthvector λ α that includes the polarization of the photon mode and introduces the couplingto the total dipole (cid:80) Nk =1 r k of the electronic system. The external variable for the photonsystem is the time-dependent current j ( α )ext ( t ).QEDFT is structurally similar to time-dependent density-functional theory in thatit is based on a one-to-one correspondence between internal and external variables. If nowphotons are also considered, the set of internal variables has to be expanded. In the frame ofEq. 8, the internal variables become the density n ( r , t ) and the mode-resolved contributionsto the electric displacement field q α ( t ). By introducing and exploiting the bijective mappingof these internal and the external variables ( v ext ( r , t ) and j ext ( t )), the auxiliary Kohn-Sham12ystem can be set up which is characterized by the electronic Kohn-Sham equations aswell as Maxwell’s equations, leading to no exchange-correlation contribution in the photonsubsystem ( j ( α )xc ( t ) = 0). This reformulation subsumes the ’quantumness’ of the light-matterinteraction solely into the local exchange-correlation potential that now features a componentdue to the electron-photon interactions, in addition to the part due to electron-electroninteraction part, i.e., v xcσ ( r , t ) = v eexcσ ( r , t ) + v epxcσ ( r , t ). Extending the coupled Maxwell-Kohn-Sham equations to consider quantum photons is thus condensed into the calculationof the local exchange-correlation potential.In practice, QEDFT requires the construction of an additional exchange-correlation po-tential that describes the electron-photon interaction. First attempts to construct this po-tential are based on many-body perturbation theory within the exact exchange approxima-tion by utilizing the optimized-effective potential (OEP) method. A. The optimized effective potential (OEP)
The OEP photon energy has been introduced in Ref. 29 and depends on occupied andunoccupied orbitals of the system. Alternatively, the OEP photon energy can also be formu-lated using occupied orbitals and orbitals shifts only. The full energy expression is givenby E xc = E ( ee ) xc + N p (cid:88) α =1 E ( α ) x , (9)where E ( ee ) xc describes the electronic exchange-correlation energy and E ( α ) x the exchange en-ergy due to the interaction of the electrons with the photon mode α . Avoiding unoccupiedorbitals is computationally much more favorable for larger systems and the photonic inducedexchange energy can be correspondingly expressed as sum over the N σ occupied orbitals ofspin channel σE ( α ) x = (cid:88) σ = ↑ , ↓ N σ (cid:88) i =1 (cid:114) ω α (cid:104) Φ (1) iσ,α | ˆ d α | ϕ iσ (cid:105) + 14 (cid:104) Φ (2) iσ,α | ˆ d α | ϕ iσ (cid:105) + c.c. , (10)where ω α describes the α th mode of the electromagnetic field and ˆ d α = λ α · r describesthe dipole operator and the electron-photon coupling strength. We can now reformulatethe problem in terms of two electron-photon orbital shifts. The Kohn-Sham orbitals ϕ iσ (1) iσ,α and Φ (2) iσ,α that can be calculatedusing the Sternheimer equations. The first electron-photon orbital shift can be obtainedexplicitly by the solution of a linear equation, i.e., a Sternheimer equation (cid:104) ˆ h sσ − ( (cid:15) iσ − ω α ) (cid:105) Φ (1) iσ,α ( r ) = − (cid:114) ω α d α ϕ iσ ( r ) + (cid:114) ω α N σ (cid:88) k =1 d ( α ) kiσ ϕ kσ ( r ) , (11)with the matrix element d ( α ) ijσ = (cid:104) ϕ iσ | ˆ d α | ϕ jσ (cid:105) . In contrast, the second electron-photon orbitalshift Φ (2) iσ,α ( r ) can be defined explicitly as followsΦ (2) iσ,α ( r ) = ˆ d α ϕ iσ ( r ) − N σ (cid:88) k =1 d ( α ) kiσ ϕ kσ ( r ) . (12)From Eq. (10), we can now deduce the potential using v xcσ ( r ) = δE xc δn σ ( r ) . (13)Doing these reformulations, we find for the final OEP equation including electron-electroneffects as well as electron-photon effects N σ (cid:88) i =1 ψ ∗ iσ ( r ) ϕ iσ ( r ) − Λ iσ ( r ) + c.c. = 0 . (14)where the homogeneity Λ iσ ( r ) is given byΛ iσ ( r ) = 12 N p (cid:88) α =1 (cid:34) Φ (1) ∗ iσ,α ( r )Φ (1) iσ,α ( r ) − (cid:104) Φ (1) iσ,α | Φ (1) iσ,α (cid:105) ϕ ∗ iσ ( r ) ϕ iσ ( r ) (cid:35) In Eq. 14 we defined a third orbital shift, the exchange-correlation orbital shift, that willbe used to obtain the corresponding exchange correlation potential and is defined along thelines of the orbital shift usually used in OEP calculations. We can obtain ψ ∗ iσ ( r ) using aSternheimer equation (cid:16) ˆ h sσ − (cid:15) iσ (cid:17) ψ ∗ iσ ( r ) = M ∗ iσ ( r ) − (cid:104) M iσ | ϕ iσ (cid:105) ϕ ∗ iσ ( r ) , (15)where M ∗ iσ ( r ) now consist of the electron-photon orbital shifts and the Kohn-Sham orbitals,as described in Ref. 13. Accordingly we define this quantity as M ∗ iσ ( r ) = − ( v xσ ( r ) − u xiσ ( r )) ϕ ∗ iσ ( r ) (16)+ N p (cid:88) α =1 (cid:34) ˆ d α (cid:32)(cid:114) ω α (1) ∗ iσ,α ( r ) + 12 ˆ d α ϕ ∗ iσ ( r ) (cid:33) − N σ (cid:88) k =1 d ( α ) ikσ (cid:32)(cid:114) ω α (1) ∗ kσ,α ( r ) + ˆ d α ϕ ∗ kσ ( r ) (cid:33)(cid:35) .
10 0 10y (˚A) − − x ( ˚A ) (a) n ( r ) −
10 0 10y (˚A) − × − ∆ n ( r )(b) OEPKLI 0 50 100 150 − − − (c) | S ( r ) | constBB Figure 3. Ground-state density of two sodium dimers affected by the vacuum-field of a cavity:In (a) we show the electron density, in (b) the difference of the electron density inside and outsidethe cavity for the OEP and the KLI approximations (see main text for details). We find for theexchange energy E ( α ) x = 6 .
52 meV and the number of photons in the correlated light-matter ground-state n pt = (cid:104) a † α a α (cid:105) = 2 . × − for the OEP case, and E ( α ) x = 6 .
67 meV and n pt = 2 . × − forKLI. Panel (c) shows the convergence of | S ( r ) | as defined by Eq. 17 using a constant with c = 20,and using the Barzilai-Borwein (BB) scheme. The simulation is set up as described in Refs. 12 and13, with (cid:126) ω α = 2 .
19 eV, λ α = 2 .
95 eV / nm − , λ α = λ α e y , and the real-space grid is sampled as31 . × . × .
75 ˚A , with a grid spacing of 0 .
265 ˚A. and include here effects of the electron-electron interaction in the quantity u xiσ ( r ). Forinstance in exchange-only calculations this quantity is defined as u xiσ ( r ) = ϕ ∗ iσ ( r ) δE ( ee ) x [ { ϕ jτ } ] δϕ iσ ( r ) ,where E ( ee ) x is the usual Fock exchange energy.Eq. 15 has to be solved self-consistently with Eq. (11). By this reformulation, we havereplaced the problem of calculating the OEP equation using all unoccupied states by aproblem of solving N p +1 Sternheimer equations that each only invoke occupied orbitals. Inthis way, the formulation of the problem becomes similar to the one of Ref. 31 for electronsonly, and which can be easily extended.For the practical implementation, we reformulate the OEP equation in the following form,as it is commonly done to construct the electronic OEP : S σ ( r ) = N σ (cid:88) i =1 ψ ∗ iσ ( r ) ϕ iσ ( r ) − Λ iσ ( r ) + c.c. (17)15nd update the potential with v ( new ) xσ ( r ) = v ( old ) xσ ( r ) + c ( r ) S σ ( r ) . (18)The quantity S σ ( r ) becomes a measure for convergence, since it is vanishing in the caseof convergence (compare Eq. (17) and Eq. (14)). For the function c ( r ), we have differentpossibilities, like using a constant or using the inverse of the electron density, as used inRef. 31. Other methods are also possible, such as the Barzilai-Borwein method. We foundstable algorithms when using a constant and when using the Barzilai-Borwein method.While we will show the computational feasibility of this approach in the following, oftena simplified solution is beneficial as starting point for the self-consistency procedure. Sucha simplified approximation can be deduced by reformulating Eq. (15) into the followingequivalent form: v xσ ( r ) = 12 n σ ( r ) N σ (cid:88) i =1 (cid:18) (cid:104) ϕ iσ | v xσ | ϕ iσ (cid:105) | ϕ iσ ( r ) | + (cid:20) { M ∗ iσ ( r ) + v xσ ( r ) ϕ ∗ iσ ( r ) }− (cid:104){ M iσ + v xσ ϕ iσ }| ϕ iσ (cid:105) ϕ ∗ iσ ( r ) − (ˆ h sσ ( r ) − (cid:15) iσ ) ψ ∗ iσ ( r ) (cid:21) ϕ iσ ( r ) (cid:19) + c.c. (19)and subsequently assuming (ˆ h sσ ( r ) − (cid:15) iσ ) ψ ∗ iσ ( r ) = 0 to start the iterative process. In situa-tions where Λ iσ ( r ) = 0, such as pure electronic exact exchange, this approximation is exactfor a single electron and is referred to as the Krieger-Li-Iafrate (KLI) approximation. By multiplying Eq. (19) by | ϕ jσ ( r ) | and integrating over the spatial coordinates, we arriveat a linear equation which in turn can be solved for the approximate v xσ . This approxima-tion scheme has proven to often deliver sufficiently accurate results for electronic structurecalculation with significantly lower computational effort and reliable stability. As this ap-proximation can be seen as a diagonal approximation to the response function, it unavoidablyfails in accurately describing polarization features. In the context of light-matter correlatedground-states, this leads to a slight unbalancing when approximating components includingphotonic excitations (introduced by Φ (1) ) and self-polarization interaction (introduced byΦ (2) ). This results in a violation of translational invariance and introduces an ar-tificial dependence on the permanent dipole. When performing KLI calculations includinglight-matter interaction, we thus suggest moving the set of coordinates into the electroniccenter-of-charge instead of the center-of-mass frame. To reduce the effect of this dependenceduring a self-consistent calculation, the optional input parameter KLIpt coc has been intro-16uced in the code. When activated, this option defines the dipole operator ˆ d α with respectto the electronic center-of-charge and can improve the stability of the algorithm.Finally, in Fig. 3 we show the capabilities of the new implementation, where we calculatetwo sodium dimers in the weak coupling regime under light-matter coupling. Figure 3 (a)shows the electron density and Fig. 3 (b) shows the comparison of OEP and KLI results.As shown in Ref. 13, in the weak-coupling regime the KLI is close to the OEP result. InFig. 3 (c) we show the convergence behavior when using a constant in Eq. (18) and whenusing the Barzilai-Borwein method. Extensions of QEDFT to the regime of vibrational strong coupling, the linear-responseregime, as well as multitrajectory methods that capture quantum fluctuations are cur-rently work in progress and will further strengthen the capabilities of the Octopus code forthe real-space description of strong light-matter coupled systems. To describe effects of theultra-strong coupling regime, one can use an alternative method that is presented in thenext section. IV. DRESSED REDUCED DENSITY MATRIX FUNCTIONAL THEORYFOR ULTRA-STRONGLY COUPLED LIGHT-MATTER SYSTEMS
The accurate description of the (ultra-)strong coupling regime of light-matter systemsis a formidable task. In many cases the known functionals for QEDFT (see Sec. III)are inaccurate and for complex electronic systems typical few-level approximations becomeunreliable.
In this section, we present the Octopus implementation of an alternative real-space ab-initio method for coupled light-matter systems. Dressed Reduced Density MatrixFunctional Theory (dressed RDMFT) extends standard electronic RDMFT to coupledlight-matter systems similarly to how QEDFT extends DFT. First tests on simple modelsystems suggest that dressed RDMFT remains accurate from the weak to the ultra-strongcoupling regime. A proper introduction of the theory, examples, and convergence studiescan be found in Ref. 44.This approach allows for the description of an interacting N -electron system coupled toone photonic mode within the dipole approximation. The respective Hamiltonian is givenin Eq. (8) of Sec. III. Note that we set j ext = 0 throughout this section and that the groundstate Ψ of Eq. (8) depends on 4 N + 1 coordinates, i.e. , Ψ = Ψ( r , σ , ... r N , σ N , p ), where17 σ i } denote the spin degrees of freedom and p is the elongation of the photon mode.Within dressed RDMFT, the original Hamiltonian (8) of N electrons in d dimensions andone mode is replaced with an extended auxiliary Hamiltonian of N dressed fermions in d + 1dimensions with coordinates z = ( r , q ) ∈ R d +1 . This auxiliary Hamiltonian readsˆ H (cid:48) = N (cid:88) k =1 (cid:2) − ∆ (cid:48) k + v (cid:48) ( z k ) (cid:3) + (cid:88) k (cid:54) = l w (cid:48) ( z k , z l ) (20)and gives access to the same physics (see Ref. 44, Sec. 4). For a d = 1 matter subsystem,the operators introduced in Eq. (20) read: the dressed Laplacian ∆ (cid:48) = ∂ ∂x + ∂ ∂q , the dressedlocal potential v (cid:48) ( z ) = v ( x ) + (cid:104) ω q − ω √ N λqx + ( λx ) (cid:105) , (21)and the dressed interaction kernel w (cid:48) ( z , z (cid:48) ) = w ( x, x (cid:48) ) + (cid:104) − ω √ N λqx (cid:48) − ω √ N λq (cid:48) x + λ xx (cid:48) (cid:105) . (22)The ground state Ψ (cid:48) ( z , σ , ..., z N , σ N ) = Ψ( x , σ , ..., x N , σ n ) ⊗ χ ( p , ..., p N ) of ˆ H (cid:48) is a prod-uct of the original physical ground state Ψ and the ground state of χ , which in turn is theproduct of N − γ (cid:48) ( z , z (cid:48) ) = N (cid:88) σ ,...,σ N (cid:90) d N − z Ψ (cid:48)∗ ( z (cid:48) σ , z σ , ..., z N σ N )Ψ (cid:48) ( z σ , z σ , ..., z N σ N ) . (23)To apply RDMFT theory on the auxiliary system, we have to replace the total energyfunctional of electronic RDMFT, given in Ref. 45, with the newly introduced quantities ofthe dressed system, i.e., the auxiliary Hamiltonian of Eq. (20) (approximately) evaluated bythe dressed 1RDM γ (cid:48) of Eq. (23). By that, common approximations for the two-body energyexpression in terms of the 1RDM can be directly transferred from electronic theory to thedressed system. The minimization is performed like in the electronic case and is basedon the RDMFT implementation of Octopus, though the convergence of the dressed systemrequires a more complicated protocol that can be found in the Supplement of Ref. 44. Thecurrent implementation in Octopus approximates the conditions under which the dressed18 igure 4. Differences of dressed HF (dHF) and dressed RDMFT (dRDMFT) from the exactground state energies (in hartree, left) and from the exact photon number (right) as a functionof the coupling g/ω for the stretched H molecule in the dressed orbital description. DressedRDMFT improves considerably upon dressed HF for the energy, which is especially due to thebetter description of the electronic correlation. The photon number, an example of a photonicobservable, is captured similarly with both methods. However, the auxiliary wave function exhibits also another exchange symmetrywith respect to the auxiliary coordinates, which is currently neglected. For the practicalvalidity of this approach, the reader is referred to Ref. 44 (Sec. 5 and the Supplement).As an example, we consider the one-dimensional H molecule in a soft-Coulomb potentialwith a slightly stretched bond-length of b = 2 . that is modeled by the local potential v H ( x ) = − (cid:112) ( x − b/ + 1 − (cid:112) ( x + b/ + 1 + 1 √ b + 1 , (24)and the soft Coulomb interaction w ( x, x (cid:48) ) = 1 (cid:112) | x − x (cid:48) | + 1 . (25)In Fig. 4, we show the total energy and the total photon number of the dressed RDMFT,the dressed HF, and the exact many-body calculations for different coupling strengths.19e see that for small couplings, both observables are captured well by dressed RDMFTand dressed HF. With increasing coupling strength, both approximations fail to capture thestrongly increasing photon number. For the total energy instead, dressed RDMFT remainsvery close to the exact result, whereas the deviations to dressed HF increase with increas-ing coupling strength. This shows the potential of dressed RDMFT to describe correlatedelectron systems that are strongly coupled to a cavity mode. In the future, we plan toinvestigate better approximations to the polaritonic N-representability conditions that alsoaccount for the symmetry of the many-body wave function with respect to the exchange ofphoton coordinates. V. TOWARDS DYNAMICS OF STRONGLY CORRELATED SYSTEMS:THE TDDFT+ U FUNCTIONAL
It was soon recognized that the standard local and semilocal functionals of DFT tend toover-delocalize the electrons, as usual approximations are based on the homogeneous electrongas. This leads to several failures of DFT for materials in which the localization of electronsplays a critical role in dictating the system’s properties. This is, for instance, the case fortransition metals oxides. The DFT+ U method was originally proposed to compensate forsome of the failures of the LDA for such materials. In essence, the DFT+ U method aimsat a better description of the local electron-electron interaction, which is achieved by addingthe mean-field Hubbard model on a chosen localized subspace to the DFT total energy.The double counting of electron interaction in this localized subspace is then removed. TheDFT+ U total energy functional reads E DFT+U [ n, { n I,σmm (cid:48) } ] = E DFT [ n ] + E ee [ { n I,σmm (cid:48) } ] − E dc [ { n I,σmm (cid:48) } ] , (26)where E ee is the usual electron-electron interaction energy, and E dc accounts for the doublecounting of the electron-electron interaction already present in E DFT . This double-countingterm is not known in the general case and this is a general problem to all +U methods.Several approximated forms have been proposed along the years.
In the Octopus code,we implemented the most commonly used double-counting terms: the fully localized limit(FLL) and the around-mean field (AMF) double-counting terms. They respectively read20s E FLL dc [ { n I,σmm (cid:48) } ] = U N ( N − − J (cid:88) σ N σ ( N σ −
1) (27)and E AMF dc [ { n I,σmm (cid:48) } ] = U N ↑ N ↓ − ( U − J ) 2 l l + 1 (cid:88) σ N σ , (28)where N σ = (cid:80) m n σmm and N = N ↑ + N ↓ . The E ee and E dc energies depend on the densitymatrix of a localized orbital basis set { φ I,m } , which are localized orbitals attached to theatom I . In the following we refer to the elements of the density matrix of the localizedbasis as occupation matrices and we denote them { n I,σmm (cid:48) } . In the rotational-invariant formof DFT+U proposed by Dudarev et al. , and for the FLL double-counting term, we obtainthe E U energy to be added to the DFT total energy, which only depends on an effectiveHubbard U parameter U eff = U − J : E U [ { n I,σmm (cid:48) } ] = E ee [ { n I,σmm (cid:48) } ] − E dc [ { n I,σmm (cid:48) } ] = (cid:88) I,n,l U eff I,n,l (cid:88) m,σ (cid:16) n I,n,l,σmm − (cid:88) m (cid:48) n I,n,l,σmm (cid:48) n I,n,l,σm (cid:48) m (cid:17) , (29)where I is an atom index, σ is the spin index, and n , l , and m refer to the principal,azimuthal, and angular quantum numbers, respectively. In the case of a periodic system,the occupation matrices n I,n,l,σmm (cid:48) are given by n I,n,l,σmm (cid:48) = (cid:88) n BZ (cid:88) k w k f σn k (cid:104) ψ σn, k | φ I,n,l,m (cid:105) (cid:104) φ I,n,l,m (cid:48) | ψ σn, k (cid:105) , (30)where w k is the k -point weight and f σn k is the occupation of the Bloch state | ψ σn, k (cid:105) . Here, | φ I,n,l,m (cid:105) are the localized orbitals that form the basis used to describe the electron local-ization. Details of the implementation can be found in Ref. 58. We recently extended ouroriginal implementation to be able to construct a localized subspace from localized states,such as Wannier states, and to treat the intersite interaction. In its usual formulation, the DFT+ U method is an empirical method, in which the effec-tive U is a parameter of the calculation. However, it recently became possible to evaluate U and J fully ab initio and self-consistently, using the ACBN0 functional. We also imple-mented this method in Octopus and extended it to the time-dependent case in order tobe able to investigate strongly correlated materials out-of equilibrium. We showed that theabsorption spectra of transition metal oxides, such as NiO or MnO, are well reproduced byour TDDFT+ U simulations. igure 5. Self-consistent dynamics of Hubbard U for the Ni 3 d orbitals (bottom panel) for pumpintensities as indicated. The top panel represents the time-dependent vector potential and thevertical dashed lines indicates the extrema of the vector potential, i.e. , the minima of the drivingelectric field. Results from Ref. 63. Figure 5 shows the calculated time profile of the effective Hubbard U eff = U − J forthe 3 d orbitals of Ni, for light-driven NiO. The top panel shows the time profile of thedriving vector potential. This shows that strongly-driven correlated materials cannot bedescribed by simply assuming that the effective electronic parameters (here, the effectiveHubbard U ) remain constant out of equilibrium. Moreover, the possibility to tune theseeffective electronic parameters offers opportunities for light-driven phase transitions, suchas light-induced magnetic Weyl semimetals. VI. VAN DER WAALS INTERACTIONS
The van der Waals (vdW) interactions arise from correlations between electrons and arein principle described by the exact energy functional through the exchange and correlationenergy E xc [ n ( r )]. However, the vdW interactions are inherently non-local and long-range22nd, by construction, cannot be described by usual local and semi-local functionals. There-fore, much work has been devoted to finding consistent ways to enhance available functionalsto correctly describe the necessary correlations yielding vdW forces.Within a pure DFT approach, the inclusion of vdW interactions should be done throughthe exchange and correlation functional. To this end, a family of non-local density functionalshas been proposed. The so-called vdW density functionals (vdW-DF) are derived ab-initio from the screened response of the homogeneous electron gas. Another route to include vdWinteractions in DFT calculations is by adding explicit corrections to the energy and forcesbased on atomic parametrizations. The most well-known method of this type is probablythe vdW-D3 scheme of dispersion corrections from Grimme. Although computationallymore favorable than the vdW-DF method, the major disadvantage of this type of approach,particularly from the perspective of time-dependent simulations, is that it fails to correctlydescribe effects that cannot simply be understood from atomic configurations. A trivialexample would be a lone electron traveling through space. Another scheme based on explicitcorrections to the energy and forces was proposed by Tkatchenko and Scheffler (vdW-TS). This scheme has the advantage of retaining much of the low computational cost of the vdW-D3 scheme, while making the atomic parametrization dependent on the electronic density.Recently, these three schemes (vdW-DF, vdW-TS, and vdW-D3) were implemented in theOctopus code to deal with vdW interactions in isolated and periodic systems.
A. vdW-DF
Octopus supports vdW-DF functionals through the libvdwxc library. The vdW-DFfunctionals are expressed as a sum, E vdW-DFxc [ n ] = E LDAc [ n ] + E GGAx [ n ] + E nlc [ n ] , (31)of the LDA correlation energy, a GGA exchange energy, and a fully non-local correlationterm. The latter is the integral over a kernel function φ ( q , q (cid:48) , r ), E nlc = 12 (cid:90) (cid:90) n ( r ) φ ( q ( r ) , q ( r (cid:48) ) , | r − r (cid:48) | ) n ( r (cid:48) ) d r d r (cid:48) , (32)where q ( r ) depends on the local density and its gradient. Explicit evaluation of this 6-dimensional integral is very expensive and scales as volume squared, O ( N ). Roman-P´erez23nd Soler proposed an efficient method which approximates it as a sum of 3-dimensionalintegrals. The method works by expressing the integrand as a convolution using a lim-ited set of helper functions, then applying the Fourier convolution theorem for a scalingof O ( N log N ) and much lower prefactor. This has since become the standard method forevaluating vdW-DF functionals.As implemented in Octopus, the density is redistributed from its normal uniform grid ofarbitrary shape onto a uniform 3D grid forming a cube or parallelepiped, which is suitablefor 3D Fourier transforms. libvdwxc then evaluates the non-local energy and contributionsto the potential, relying on the FFTW library for efficient parallel Fourier transforms.After calculating the energy and potential, the potential is redistributed back to the originalform.Octopus supports the standard functionals vdW-DF1, vdW-DF2, and vdW-DF-cx, as well as other common forms that differ by substituting a different GGA exchangefunctional in Eq. (31). Some common supported variations are vdW-DF-optPBE, vdW-DF-optB88, and vdW-DF-C09. B. vdW-TS
Since the vdW-TS approach depends on the density, the effect of the van der Waalsinteraction can be observed in properties other than the forces. In particular, we expect toobserve an effect in the excited states of systems that interact through vdW forces. We usethe hydrogen fluoride dimer as a simple model system for a proof-of-concept application ofthe modular implementation of the TS-vdW functional correction on TDDFT calculations.The dimer geometry is setup as shown in Fig. 6. The hydrogen fluoride monomers areplaced in anti-parallel fashion, each one with its main symmetry axis oriented along the y -axis. The hydrogen fluoride bond in each monomer is 0.92 ˚A long and the moleculesare separated by 2.8 ˚A along the z -axis. At this distance, the van der Waals interactionbetween the monomers is the strongest, according to the TS-vdW model. We choose thisdimer model because it is a conveniently small system in which the effects of including thedispersion correction can be shown at an affordable computationally cost.To calculate absorption, we excite the system with an infinitesimal electric-field pulse,and then the time-dependent Kohn-Sham equations are propagated for 30.38535 (cid:126) /eV. The24 A b s o r p t i o n [ a . u . ] Energy [eV]
LDALDA+TS-vdW 0 1 2 3 4 5 9.1 9.2 9.3 9.4 9.5
Figure 6. Absorption cross section spectrum of the hydrogen fluoride dimer with and without vander Waals effects. A zoom-in of the absorption feature at around 133 nm illustrates a small redshift when van der Waals effects are considered. Both hydrogen fluoride molecules are placed onthe same plane at a separation of 2.8 ˚A. singlet dipole spectrum is evaluated from the time-dependent dipole moment. The strengthof the perturbation is set to 0.01 ˚A − and it is polarized in the z -axis. The time evolution iscarried out using the Enforced Time-Reversal Symmetry (ETRS) propagator, with (default)time steps of 0.03352 (cid:126) /eV.The results on Fig. 6 show a small van der Waals-induced bathochromic-like (red) shift inthe optical spectrum of the hydrogen fluoride dimer calculated with the LDA. This exampleopens the door for a new series of applications in supra-molecular chemistry, structuralbiology, polymer science, etc., that incorporate van der Waals effects on real-time electrondynamics. C. vdW-D3
Octopus also supports the DFT-D3 van der Waals correction. This correction dependsonly on the atomic positions and does not depend on the atomic density. It is implemented25n Octopus by linking to the DFT-D3 library provided by the authors. As we have validatedOctopus results with the reference data provided with the library, the results are guaranteedto be consistent with other codes that implement this correction.The only modification we did to the library was to move the very large set of coefficientsfrom a source file to a stand alone text file that is parsed only when necessary. This speedsup compilation and reduces the size of the binaries. The modified library is distributed withOctopus, so users need not compile it separately.
VII. POLARIZABLE CONTINUUM MODEL
The Polarizable Continuum Model (PCM) comprises a family of implicit-solvent ap-proaches to tackle quantum-mechanical calculations of molecules in solution. PCM as-sumes that (i) the solvent is a continuum and infinite dielectric medium characterized bya frequency-dependent dielectric function; and that (ii) a void cavity of appropriate shapeand size encapsulates the solute molecule, separating it from the solvent by a sharp inter-face. The numerical implementation of PCM relies on the Apparent Surface Charge (ASC)approach and the Boundary Element Method (BEM). In this framework, the solvent polar-ization response induced by the molecule’s charge density is modeled by a reaction potentialdefined by a set of point charges q = { q , q , ..., q T } that spread over a tessellated cavitysurface consisting of T finite surface elements or tesserae . The reaction potential is theobject through which the complex dielectric environment is accounted for. The implemen-tation of the PCM in Octopus rests on the Integral Equation Formalism (IEF). The keyIEF-PCM equation for computing the polarization charges is q = QV , (33)where q and V are T × Q is the T × T PCM response matrix, which depends on the geometry of the cavity and the dielectricfunction of the solvent. A realistic description of the molecular cavity is key to capture accurate electronic andoptical properties of molecules in solution. In principle, the cavity should (i) exclude thesolvent, (ii) comprise most of the solute electronic density, and (iii) conform with the molec-ular shape. Octopus uses the GEPOL algorithm, which builds up the van der Waals26avity from the union of interlocking spheres with element-specific radii centered at eachatom position (by default no spheres are built around hydrogen atoms). Within GEPOL,the BEM tessellation of the solute-solvent interface is done starting from a 60- or 240-facecircumscribed polyhedron per sphere, selecting only exposed tesserae and properly reshap-ing those that are partially exposed. The BEM tessellation is a surface grid constructedindependently from the real-space three-dimensional grid used in Octopus to represent boththe Kohn-Sham electronic Hamiltonian and molecular orbitals. The mismatch between thetwo discrete representations might cause numerical problems arising from the Coulomb sin-gularities whenever tesserae and grid points are close to each other. The implementation ofthe PCM in Octopus regularizes such singularities by using normalized spherical Gaussianfunctions to smooth the discretized polarization charges q . Having briefly described the implementation, we now look at specific and relevant casesfor chemical reactions, namely solvation energies and ground-state stabilization. The mostevident effect of the presence of a solvent is to change the total energy of a system com-pared to its value in vacuum. In the framework of DFT, the Kohn-Sham Hamiltonian of thesolvated molecule contains the solvent reaction potential, which is a functional of the elec-tronic and nuclear densities of the solute molecule. The latter implies that the Kohn-Shamand IEF-PCM equations become coupled and the polarization charges q e = QV Hartree [ ρ e ],induced by the molecule’s electronic density, have to be computed at each Kohn-Sham iter-ation until convergence is reached. The resulting electronic density can be used to computethe electrostatic contribution to the solvation free energy ∆ G el = G [ ρ e solv ] − E [ ρ e vac ], where G and E denote the free and total energy functionals of the molecule in solvent and in vacuum,respectively. The current PCM implementation in Octopus has been thoroughly tested ona benchmark set of organic molecules and compared with analogous calculations performedwith the quantum chemistry package Gamess . As an example of the solvent stabilizationeffect in the ground state, we studied the nitrobenzene molecule in water (static dielectricconstant set to (cid:15) s = 80).In Fig. 7, we show the convergence of the solvation free energy as a function of the Kohn-Sham iteration. The solvation free energy is stabilized already after 10 Kohn-Sham iterationsout of the 19 required to optimize the molecular orbitals of the solvated molecule. We havealso verified that the numerical error inherent to the discretization of the solute cavity surfaceis very small. For example, the total polarization charge Q PCM = (cid:80) Ti =1 q ei at each Kohn-27 [Å] Figure 7. (Left panel) Nitrobenzene molecule in water ( (cid:15) s = 80) surrounded by the ground statePCM charges distributed on the solvation cavity. The actual distribution of the charges in 3Dcan be seen in the supplementary animation SM1. (Right panel) Electrostatic contribution tothe solvation free energy of the computed at each Kohn-Sham iteration. Kohn-Sham equationswere solved in real-space as implemented in Octopus using the GGA-PBE approximation to theexchange-correlation energy. The simulation box was built using spheres of radius 5 ˚A and anuniform spacing of 0.19 ˚A between grid points. The radii of the spheres used to build the cavitysurface of the solute molecule are 2.4 ˚A for carbon, 1.8 ˚A for oxygen and 1.9 ˚A for nitrogen. Sham iteration only deviates by only 0 .
03% from the actual number of valence electronsin the nitrobenzene molecule (the relation between q and solute charge is determined byGauss’s theorem ). However, the magnitude of this error will depend in general on the sizeand the geometry of the molecule.Now we move forward and describe the extension of the previous PCM method to the timedomain, using real-time PCM and non-equilibrium solvation. Ground state PCM is unableto capture the complex dynamical interactions between the solvent and the solute when thelatter is in an excited state. Fortunately, the PCM admits a generalization to account for sol-vation dynamics. The time-dependent PCM (TD-PCM) implementation in Octopus comesin three flavors of increasing complexity, all coupled with real-time TDDFT calculations ofthe solute molecules. The first one, called equilibrium TD-PCM, assumes that the solventis fast enough to instantaneously equilibrate the solute charge density fluctuations and topolarize accordingly. This approach is physically sound for weakly polar solvents, having28imilar values for the static and dynamic dielectric constants. The second is a nonequilibriumapproach called inertial TD-PCM and amounts to partitioning the solvent response in a fast(dynamic) and a slow (inertial) part. Faster degrees of freedom respond instantaneouslyto the changes of the applied potential (either of molecular or external in origin), whereasslower degrees of freedom remain “frozen” and in equilibrium with the initial value of thefield. This approximation works well when the solvent relaxation times are large enoughwith respect to the electronic excitations in the solute molecule (e.g., within the picosecondscale). The third TD-PCM approach, called equation of motion (EOM) TD-PCM, considersthe full history-dependent evolution of the solvent polarization through a set of equations ofmotion for the polarization charges. Nonequilibrium polarization effects of this sort origi-nate from the frequency-dependent dielectric response of the solvent, encoding the fact thatit takes a different non-negligible time to adjust to fast or slow electrostatic perturbations.Solvation dynamics affect strongly and non-trivially the absorption spectrum of molecules,especially for fast and polar solvents, by inducing solvatochromic shifts of the peaks in theUV-Vis absorption spectrum and modifying their relative amplitudes. The details aboutthe implementation of all of these schemes and a detailed discussion of their effects can befound in Refs. 81 and 89.Here we show the differences among the TD-PCM flavors described above for the case ofthe nitrobenzene molecule in aqueous solution, excited with an electrical dipole perturbationin the x direction. We take the dynamic dielectric constant of water, entering in the inertialand EOM TD-PCMs, as (cid:15) d = 1 . invacuo , all of the TD-PCM methods produce shifts of the features (in this case, toward lowerexcitation energies). The shifts are not rigid overall, but depend on the excited state andon the specific solvation scheme (equilibrium, inertial, EOM). Within Debye’s model, theequilibrium and inertial TD-PCMs are limiting cases for the dynamics when the relaxationtime is zero or infinite, respectively. The EOM TD-PCM, with a finite relaxation time,interpolates between these limits. Water has a large relaxation time of 3.37 ps and, therefore,the absorption spectrum is almost coincident for the EOM and the inertial TD-PCMs, ascan be seen in Fig. 8. Whenever the solvent is as slow as water, there is not much gainin selecting EOM over inertial TD-PCM and the latter is the method of choice in terms ofcomputational performance. Instead, for faster solvents, the EOM results depart from those29f the inertial algorithm, approaching to those of the equilibrium TD-PCM. In Fig. 8, wehave considered an effective solvent having the same static and dynamic dielectric constantsas water, but with a 1 fs relaxation. The EOM TD-PCM for such a model solvent producesalmost the same excitation energy shifts of equilibrium TD-PCM and peak intensities inbetween the inertial and the equilibrium TD-PCMs, as expected.A real-time representation of the molecular dipole coupled with the PCM surface chargescan be appreciated in the supplementary movie SM2, which also highlights how the differentTD-PCM approaches affect the time-evolution of the dipole.When studying the evolution of a solvated molecule under the effect of an external time-dependent electromagnetic field, a separate treatment is required to take into account thatthere are two solvent polarization contributions interacting with the solute molecule, namely s t r eng t h f un c t i on [ / e V ] excitation energy [eV] in vacuoEquilibrium TD-PCMInertial TD-PCMEOM TD-PCMEOM TD-PCM w/ faster relax. C FF / C FF m a x excitation energy [eV] Figure 8. Absorption spectrum of nitrobenzene in vacuo and in water ( (cid:15) s = 80, (cid:15) d = 1 . max = 1 .
27 and 2.55, respectively. Real-time TDDFT simulations of 20 fswith a time step of 1.7 × − fs were performed to obtain the absorption spectra. An electric dipoleperturbation within the linear regime and oriented along x was applied as an initial perturbation.The rest of the computational details are shared with the ground-state DFT calculations leadingto Fig. 7. The reaction field comes from the polarization inducedby solute molecule itself, while the cavity field arises as a polarization induced by the externalelectric field. Both reaction and cavity field effects can be accounted for in static and real-time quantum-mechanical calculations of molecules within PCM and TD-PCM simulationsin Octopus. All of the TD-PCM calculations shown here were performed using both reaction- andcavity-field effects, although the redshift of the peaks in each TD-PCM scheme is mainly afeature of the reaction-field, as ascertained in test calculations (not shown) and as expected. Cavity field effects impact directly on the peak intensities, by making the absorption morefavorable depending on how large the effective local field acting on the molecule is allowedto be by solvent dielectric properties and the geometry of the cavity (normally, reassem-bling the molecular shape). Still, cavity field effects can have a non-trivial influence inthe absorption spectrum shape when considering non-equilibrium solvation dynamics (EOMTD-PCM), by changing the relative peak intensities. This effect can be seen in the inset ofthe Fig. 8, where the normalized cavity field factor (CFF) – the ratio between absorptioncross-section with and without cavity field effects – is plotted against the excitation energyfor the EOM TD-PCM simulations with water and the aforementioned faster water-like sol-vent. The large difference in absorption peak intensities between the equilibrium and therest of the TD-PCMs in Fig. 8 is also related to the modification of the probing electro-magnetic field inside the dielectric. The photo-absorption cross-section is, by definition, theratio between the absorbed and the incoming power of light. In our case, both increasewith the dielectric constant: (i) the larger the dielectric constant of a medium the largerthe CFF; but also (ii) a larger dielectric constant implies a smaller phase velocity of light,therefore increasing the power of the traveling electromagnetic wave. The impact of the lat-ter effect ( ∝ refractive index = (cid:112) ( | (cid:15) ( ω ) | + (cid:60) (cid:15) ( ω )) /
2) is stronger than the former ( ∝ CFF,e.g., for a spherical cavity ∝ (cid:12)(cid:12)(cid:12) (cid:15) ( ω )2 (cid:15) ( ω )+1 (cid:12)(cid:12)(cid:12) ) for a large enough dielectric constant such as theone corresponding to the equilibrium TD-PCM for water ( (cid:15) s = 80).In conclusion, Octopus is now capable of including an implicit dielectric continuum modelas an environment for quantum mechanical calculations of excited states in the time-domain.The different versions of the TD-PCM scheme implemented allow us to select the mostsuitable to capture the relevant physics accounting for a full range of different relaxationand response times. 31 III. MAGNONS FROM REAL-TIME TDDFT
In the last couple of years the first studies investigating magnetization dynamics fromfirst principles in real time have emerged.
Here we are interested in transverse magneticexcitations, specifically magnons, which are long wavelength collective excitations with atypical energy of tens to hundreds of milli electronvolts. We recently developed an alternativeto the linear-response TDDFT formulation to compute the spin susceptibilities ofmagnetic systems, based on real-time TDDFT. Our approach follows the work of Bertsch et al. for optical excitations. In the original work of Bertsch et al. the system is perturbedby a sudden change in the vector potential, which induces a charge- and current-densityresponse. To investigate magnons, we employ a “transverse magnetic kick”, which inducesmagnetic fluctuations in the system. To be more precise, our perturbation corresponds toan infinitely short application of a Zeeman term, δ ˆ H q ( t ) = B δ ( t ) (cid:90) d r (cid:2) e − i q · r ˆ σ + ( r ) + e i q · r ˆ σ − ( r ) (cid:3) , (34)where q is the momentum of the spin wave we are exciting, B is the strength of the pertur-bation, and σ ± = σ x ± iσ y are Pauli matrices. In this expression, the z axis is taken to bealong the direction of the magnetization of the system before excitation. In case the systemhas a preferred magnetization direction due to the presence of spin-orbit coupling, usuallyreferred as an “easy axis”, we perform a kick in the transverse direction with respect to this“easy axis” of the system.The subsequent time evolution of the spin magnetization m ( r , t ), governed by the time-dependent Schr¨odinger equation, is then computed and analyzed in Fourier space to spinsusceptibilities. If we perturb our system from its ground state and we assume linear re-sponse, we directly have that m + ( q ; ω ) = χ + − ( q ; ω ) B , (35)where χ + − ( q ; ω ) is the spin susceptibility we want to extract. A similar expression is ob-tained for χ − + ( q ; ω ).In order to access finite momenta which are a fraction of the Brillouin zone, one typicallyhas to employ large supercells to perform the dynamics, which can be computationally veryexpensive when a few meV energy resolution is needed. Fortunately, there is a way to cir-cumvent the construction of supercells, the so-called generalized Bloch theorem (GBT). The32BT has been introduced by Sandradskii for the calculation of ground-state spin wavesand requires the implementation of specific boundary conditions. Therefore, we also imple-mented the GBT for the real-time calculation of magnons, taking advantage that Octopusis a real-space finite differences code, for which we can easily specify any type of bound-ary conditions. However, it is important to note that this applies only in the absence ofspin-orbit coupling. The boundary condition, described below, depends on the momentum q of the perturbation and acts differently depending on whether a state was originally “up”or “down” with respect to the unperturbed magnetization. This is determined by the signof (cid:104) Φ n k | S z | Φ n k (cid:105) just before the perturbation, for each spinor state | Φ n k (cid:105) . If we label thesestates α and β , the boundary condition readsΦ α, k n ( r , t ) = e i k · r u ↑ α, k n ( r , t ) e i q · r u ↓ α, k n ( r , t ) , (36)Φ β, k n ( r , t ) = e i k · r e − i q · r u ↑ β, k n ( r , t ) u ↓ β, k n ( r , t ) . We checked that performing the simulation using the GBT or the supercell approachwas leading to the same results, up to numerical precision. We tested our approach againstcubic Ni, Fe, and Co, which have been widely studied using linear response TDDFT, andwe found that our results are in very good agreement with previous works, thus validatingour implementation. Figure 9 shows the results obtained for bulk nickel using the adiabaticlocal-density approximation. We used here a real-space grid spacing of 0.27 bohr, norm-conserving pseudopotentials, and we employed a 16 × × k -point grid shifted 4 times inorder to resolve momenta q which are multiples of 2 π/ (16 a ), where a is the lattice parameterof Ni, which we took to be 3 .
436 ˚A. In order to obtain the response within linear response,we took B = 0 .
02 and we propagated during 435.4 fs, using a time-step of 1.81 as. To reducethe numerical burden, symmetries were employed.Let us now comment on the interest of the proposed approach. So far, we only used thisnew method for investigating weak magnetic kicks, where we recover the results from linearresponse theory. However, our approach does not rely on the assumptions of small perturba-tions and can be used to investigate nonlinear phenomena induced by strong magnetic field,as well as out-of-equilibrium situations where the system is kicked from an excited state,which is a strength of a real-time method.
Moreover, we can directly investigate the33 ΓX ΓX ΓX X m ag n o n f r e q u e n c y [ m e V ] ω ( q ) = Dq ( D ≈ .
89 eV˚A )EXP − − − l og | I m [ χ + − ( q , ω ) ] | Figure 9. Calculated magnon dispersion: Spin susceptibility of bulk Ni along the Γ X direction,displayed in logarithmic scale. Numerical details are given in the main text. The experimentaldata are taken from Ref. 101. coupling with other degrees of freedom, such as phonons or photons, without any new theoryor code development. Using the real-space method, we benefit from the favorable scalingof the time propagation, which is linear in the number of states, whereas sum-over-statesapproaches usually scale quadratically with the number of states. Finally, our approachoffers the great advantage of not requiring the use of an exchange-correlation kernel, as onlythe exchange-correlation potential is needed to perform a time propagation. This is veryinteresting in order to test new functionals and theory levels, for which deriving the expres-sion of the exchange-correlation kernel can become complicated. These different aspects willbe investigated in future works. IX. ORBITAL MAGNETO-OPTICAL RESPONSE OF SOLIDS ANDMOLECULES FROM A STERNHEIMER APPROACH
In the present section we review the implemented routines for magneto-optical phenom-ena, which arise from the loss of symmetry between left and right circularly polarized lightin the presence of a magnetic field.
While magneto-optical spectra can be straight-forwardly computed for molecules, the theory for solids has been developed only re-34ently.
The reason is that external electromagnetic fields break the translational symmetryof periodic systems, which is formally expressed through the unboundness of the position op-erator.
The orbital response to magnetic fields is especially complicated to describe, assuch fields lead to non-perturbative changes in wavefunctions and introduce vector couplingto electron dynamics.
Here we will focus on calculations of changes in the linear opticalresponse in the presence of the magnetic field within the Sternheimer approach.
To treat uniform magnetic fields in periodic systems, we use the approach based onperturbation theory for the one-particle density matrix.
Such an approach allows usto work under purely periodic boundary conditions and to automatically take into accountgauge invariance. A periodic and gauge-invariant counterpart ˜ O is distinguished for anyoperator O = O r r defined for two points r and r in real space O r r = ˜ O r r exp (cid:18) − ic (cid:90) r r A ( r )d r (cid:19) , (37)where A is the vector potential associated with external electric ( E = − c − ∂ t A , c is thespeed of light) and magnetic ( B = ∇ × A ) fields and the integral is taken along the straightline between points r and r .The time-dependent Liouville equation for the gauge-invariant counterpart ˜ ρ of the one-particle density matrix to the first order in E and B takes the form − i∂ t ˜ ρ + [ H , ˜ ρ ] = − (cid:26) E + 1 c V × B , [ r , ˜ ρ ] (cid:27) − [ δ ˜ H, ˜ ρ ] . (38)Here the commutator and anticommutator of operators O (1) and O (2) are denoted as[ O (1) , O (2) ] and {O (1) , O (2) } , respectively, the velocity operator V = − i [ r , ˜ H ] is computedwith account of all non-local contributions to the Hamiltonian, such as from non-localpseudopotentials, and the Hamiltonian is represented as ˜ H = H + δ ˜ H , where the differ-ence δ ˜ H between the gauge-invariant counterpart ˜ H of the Hamiltonian and unperturbedHamiltonian H is related to the local-field effects. Unlike the singular position operator r ,the commutator [ r , ˜ ρ ] of the position operator with the periodic function ˜ ρ is well definedin Eq. (38) and corresponds to the derivative with respect to the wave vector, i∂ k ˜ ρ k , inreciprocal space. Differentiating Eq. (38), one finds the derivatives of the density matrix˜ ρ ( P ) = ∂ ˜ ρ/∂P with respect to external perturbations P ( E , B , etc.).For TDDFT calculations in Octopus, ρ is the Kohn-Sham density matrix. The n -thorder derivative ˜ ρ ( P ) describing the joint response to the perturbations P = P P ...P n is35ivided into four blocks within and between the occupied (V) and unoccupied subspaces (C):˜ ρ ( P )VV = P v ˜ ρ ( P ) P v , ˜ ρ ( P )CC = P c ˜ ρ ( P ) P c , ˜ ρ ( P )VC = P v ˜ ρ ( P ) P c , and ˜ ρ ( P )CV = P c ˜ ρ ( P ) P v , where P v = ρ (0) and P c = 1 − P v are the projectors onto the occupied and unoccupied bands. In accordance withthe density matrix perturbation theory, to get the elements ˜ ρ ( P )CV , Eq. (38) is projectedonto unperturbed Kohn-Sham wavefunctions | ψ (0) v k (cid:105) of occupied bands v to give an equationfor the response function | η ( P ) v k (cid:105) = ˜ ρ ( P )CV (Ω) | ψ (0) v k (cid:105) L v k (Ω) | η ( P ) v k (cid:105) = P c R ( P ) [ ˜ ρ ( n − , ..., ρ (0) , n ( P ) ] | ψ (0) v k (cid:105) . (39)Here the operator on the left-hand side is given by L v k (Ω) = Ω + H − (cid:15) v k , where Ω is thefrequency considered and (cid:15) v k is the energy of the unperturbed state | ψ (0) v k (cid:105) . The operator R comes from the right-hand side of Eq. (38) and is determined by the derivatives of the densitymatrix of the previous orders. If the local-field effects are taken into account, the right-handside R also depends on the derivative of the electron density n ( P ) ( r ) = ρ ( P ) ( r , r ) δ ( r − r ).In this case, Eq. (39) needs to be solved self-consistently. The calculations in practice workwith the periodic parts of the wavefunctions, | u (0) v k (cid:105) . The commutator [ r , ˜ ρ ] corresponding to i∂ k ˜ ρ k in reciprocal space is computed within the k · p theory. Equation (39) is solvedusing the efficient Sternheimer approach, where the function | η ( P ) v k (Ω) (cid:105) that fits intoEq. (39) is found iteratively at each frequency Ω. To avoid divergences at resonances, a smallbut finite imaginary frequency iδ is added to the frequency Ω of the external perturbationso that Ω = Ω + iδ .Once the solution of Eq. 39 is known, the elements ˜ ρ ( P )CV are obtained as˜ ρ ( P )CV (Ω) = (cid:90) BZ d k (2 π ) (cid:88) v | η ( P ) v k (Ω) (cid:105)(cid:104) ψ (0) v k | . (40)The elements of ˜ ρ ( P ) between the occupied and unoccupied subspaces are found using therelation ˜ ρ ( P )VC (Ω) = ( ˜ ρ ( P )CV ( − Ω ∗ )) ∗ and, for that, Eq. 39 is also solved for the frequency − Ω ∗ .To find the elements within the occupied ˜ ρ ( P )VV and unoccupied ˜ ρ ( P )CC subspaces, the idem-potency condition ρ = ρρ is used. In terms of the periodic counterpart ˜ ρ of the densitymatrix and to the first order in the magnetic field, it is written as ˜ ρ = ˜ ρ ˜ ρ + i c B · [ r , ˜ ρ ] × [ r , ˜ ρ ] . (41)The contribution α νµ,γ to the polarizability in the presence of the magnetic field ( α νµ = α νµ + α νµ,γ B γ ) is finally obtained from the current response as α νµ,γ (Ω) = i Ω Tr (cid:2) V ν ˜ ρ ( E µ B γ ) (Ω) (cid:3) , (42)36here indices ν , µ , and γ are used to denote components of the vectors V , E , and B . Thedielectric tensor in the presence of the magnetic field is computed as (cid:15) νµ = δ νµ + 4 πα νµ /w ,where w is the unit cell volume. Note that according to the “2 n + 1” theorem, there isno need to calculate explicitly the second-order derivative ρ ( E µ B γ ) . Instead, α νµ,γ is expressedthrough the first-order derivatives to the perturbations P = E µ , B γ and a supplementaryperturbation corresponding to a vector potential P = A ν at frequency − Ω. To test the developed formalism for solids, it has been applied to bulk silicon and thecorresponding results are shown in Fig. 10. The full details of the calculations can befound in Ref. 113. It is seen in Fig. 10 (b) that even without account of excitonic effects, thecalculated spectra Re / Im (cid:15) xy for the transverse component of the dielectric tensor are alreadyqualitatively similar to the experimental curves at the direct absorption edge. To modelexcitonic effects, we have used the model from Ref. 126. The spectra computed with accountof excitonic effects show a better agreement with the experimental results (Figs. 10(a) and10(b)). Although the magnitudes of the peaks in the magneto-optical spectra are about afactor of two smaller than in the magneto-optical measurements, they can be correctedby reducing the linewidth δ assumed in the calculations.In the limit of a large supercell, this formalism becomes equivalent to the simpler standardformulation for finite systems, which we have also implemented for reference. In this case,the Liouville equation for the density matrix is written in the Coulomb gauge asΩ ρ + [ H , ρ ] = [ d · E + m · B − δH, ρ ] , (43)where d = − r is the electric dipole moment, m = − r × V / c is the orbital magnetic dipolemoment, and δH describes local-field effects. The change in the polarizability α νµ,γ B γ inthe presence of the magnetic field is calculated from the dipole response α νµ,γ (Ω) = Tr (cid:2) d ν ρ ( E µ B γ ) (Ω) (cid:3) . (44)As in the periodic case, based on the “2 n + 1” theorem, explicit calculation ofthe second-order derivative ρ ( E µ B γ ) is avoided. Instead, a supplementary electric field atfrequency − Ω, E (cid:48) ν , is introduced. In the case of real wavefunctions, as can be chosenfor finite systems in real space in the presence of time-reversal symmetry, | η ( E (cid:48) ν ) v ( − Ω) (cid:105) =( | η ( E ν ) v ( − Ω ∗ ) (cid:105) ) ∗ and there is no need to solve additionally the Liouville equation for thissupplementary electric field. 37 igure 10. Calculated components (a) (cid:15) xx and (b) (cid:15) xy of the dielectric tensor of silicon for differentfrequencies of light Ω (in eV) and the magnetic field of 1 T along the z axis (with and withoutaccount of excitonic effects). The linewidth δ = 0 . The experimental data for (cid:15) xx and (cid:15) xy (scaled by 1/2) are indicated by symbols. The efficiency of the present scheme to compute magneto-optics is comparable to standardlinear-response calculations of simple optical polarizability in the absence of the magneticfield.
When local-field effects are included self-consistently, the calculations of magneto-optical spectra for solids take only twice as long as those of polarizability. For finite systems,the computational effort is the same as for the simple optical polarizability.
X. TIME-DEPENDENT ANGULAR RESOLVED PHOTOELECTRONSPECTROSCOPY
The real-space representation of the dynamics of the electronic structure allows for aseamless and straightforward description of dynamics processes outside the material, i.e.,processes where electrons are excited into the vacuum. Beyond simply describing the ioniza-tion process, Octopus has routines implemented that compute photoelectron spectroscopy in38ifferent flavors. Photoelectron spectroscopy is particularly ubiquitous for the characteriza-tion of the electronic structure in solids, because it provides a direct observable of the energyand momentum distribution of electronic states, known as the bandstructure. In contrastto other electronic structure methods that compute quantities linked to the photoelectronspectrum, most commonly the single-quasiparticle spectral function obtained through the GW approximation to the many-body self-energy, the approach described here does not con-sider unit cells of bulk materials, but instead directly computes the energy and momentumresolved ionization probability. Besides the study of the electronic structure of solids,photoelectron spectra of atoms and molecules are of equal interest and the implementationin Octopus is capable to describe accurately all such systems on equal footing.The most detailed quantity available in the experiments is the momentum-resolved photo-electron probability P ( p ), i.e., the probability to detect an electron with a given momentum p . In some cases the experimental setup offers the possibility to measure the spin polar-ization along a given axis. The formalism we are going to outline can easily accommodatea non-collinear spin structure and, therefore, calculate spin-resolved quantities, but for thesake of simplicity, in the following, we are going to restrict ourselves to closed shell systemswith collinear spins, i.e., where all the orbitals are doubly occupied with electrons havingopposite spins. The reader interested in the most general case can find a detailed descriptionin Ref. 131. From P ( p ) one can obtain other derived quantities by simple manipulation. Forinstance, the energy-resolved spectrum P ( E ), used to identify the occupied energy levels,can be obtained with a change of variable using the free electron energy dispersion relation E = p /
2. The angle-resolved photoelectron spectrum (ARPES), normally employed tomeasure the quasiparticle bandstructure, is simply obtained by taking the energy resolvedspectrum as a function of the electron momentum parallel to the surface P ( p (cid:107) , E ).The t-SURFF method was first proposed by Scrinzi for one-electron systems and laterextended to many electrons with TDDFT for periodic and non-periodic systems. It isbased on the assumption that the Kohn-Sham Hamiltonian describing the full experimentalprocess, i.e., including the ionization and detection, can be decomposed into the sum of twoHamiltonians acting into complementary spatial regions, inner and outer , and that can beapproximated in different ways. In the inner region surrounding the system the electrondynamics is governed by the interacting Kohn-Sham Hamiltonian ˆ H KS ( r , t ) while in the outer region, electrons are free from the Coulomb tails of the parent system and behave as39ndependent particles driven by an external field and therefore are described by the VolkovHamiltonian ˆ H V ( r , t ) = 1 / − i ∇ − A ( t ) /c ) . We express the field with a time dependentvector potential potential A ( t ) in the dipole approximation, i.e. by discarding the spatialdependence of the field A ( r , t ) ≈ A ( t ). The advantage of this approach is provided by thefact that the time-dependent Schr¨odinger equation associated with the Volkov Hamiltonian,can be solved analytically and that the solutions, the Volkov waves χ p ( r , t ) = 1(2 π ) e i p · r e i Φ( p ,t ) , (45)with Φ( p , t ) = (cid:90) t d τ (cid:18) p − A ( t ) c (cid:19) , (46)are eigenstates of the momentum operator. We can therefore expand the Kohn-Sham orbitals ϕ i as a superposition of detector states in the form of Volkov waves ϕ i ( r , t ) = (cid:90) d p b i ( p , t ) χ p ( r , t ) (47)and obtain the photoelectron probability in terms of the expansion coefficients by summingup the contribution of all the orbitals: P ( p ) = lim t →∞ /N (cid:80) N/ i =1 | b i ( p , t ) | .Using the continuity equation, t-SURFF allows us to express the coefficients b i , and thus P ( p ), as a time integral of the photo-current flux through the surface S separating thedomains of the two Hamiltonians. More specifically b i ( p , t ) = − (cid:90) t d τ (cid:73) S d s · (cid:104) χ p ( τ ) | ˆ j | ϕ i ( τ ) (cid:105) , (48)with the single-particle current density operator matrix element given by (cid:104) χ p ( τ ) | ˆ j | ϕ i ( τ ) (cid:105) = 12 (cid:26) iϕ i ( r , τ ) ∇ χ ∗ p ( r , τ ) − iχ ∗ p ( r , τ ) ∇ ϕ i ( r , τ ) − A ( t ) c χ ∗ p ( r , τ ) ϕ i ( r , τ ) (cid:27) . (49)The flexibility offered by the definition of the boundary surface S allows us to easilyadapt t-SURFF to periodic and non-periodic systems. For periodic systems we choose S as a plane parallel to the material surface, while for the non-periodic case the most naturalchoice is a sphere as shown in Fig. 11 (a) and (b).For periodic systems one can apply the Bloch theorem to the Volkov Hamiltonian andfurther reduce Eq. (48) to an expression where only the periodic part of the Bloch waves areemployed. (b)(a) S periodic dimensions n o n - p e r i o d i c d i m e n s i o n (d)(c) -2 -1 0 1 2x10(b) p y [ a . u .] -3-2-10123 3210-1-2-3 p x [a.u.]min max -16-12-8-4 min max Figure 11. Calculated photoelectron spectrum with t-SURFF, as implemented in Octopus.
The geometries employed to calculate the flux of the photoelectron current are depicted in (a)for periodic and (b) non-periodic systems. (c) Pump-probe ARPES spectrum of monolayer h-BNdriven by a laser field resonant with the gap at K, adapted from Ref. 131. (d) Photoelectronangular distribution obtained by strong-field ionization of C with an IR field polarized along x ,adapted from Ref. 133. All the electrons rescattering at the same time in one period of the field endup with final momenta forming a ring (dashed) centered at the value of A ( t r ) at the rescatteringtime, t r . The two arrows represent the graphical decomposition of the final momentum of thephotoelectron in the vector potential at the moment of rescattering (horizontal arrow) plus therescattering momentum. S for the evaluation of the flux andin practice this has to be converged by placing it at successively larger distances form thesystem. This method needs an explicit description of the vacuum around the probed system,which for solids requires an explicit construction of the electronic structure of the surfacelayer and, if bulk properties are probed, one needs to converge a slab of material. Comparedto unit cell calculations of solids, this puts the present method at a disadvantage. However,when aiming at simulating specific experiments, it natively includes signals coming from thesurface layers and thus can capture more of the experimental reality of the spectroscopicprocess compared with standard perturbative many-body approaches which are directlyperformed in the bulk unit cell.In practical calculations, in order to avoid spurious reflection of electrons from the bound-aries of the simulation box, one has to employ absorbing boundaries and, depending on thecondition on their transparency, this can add up to the size of the simulation box. Otherthan the increased size of the computational box, this method is not computationally costly,as it requires only the evaluation of the gradient operator on surface points S and straight-forwardly supports existing k-points, states, and mesh parallelizations.Pump-probe configurations can be naturally simulated since there is no restriction onthe functional form of the vector potential A ( t ) and, therefore, we can accommodate anylinear combination of pulses. As an example, in Fig. 11 (c) we show the result of a simu-lated pump-probe time-resolved ARPES measurement of monolayer hexagonal boron nitride(hBN) driven by a field resonant with the gap at the K -point in the Brillouin zone. As it isapparent from the simulation, the excited state population transfer to the conduction bandis well observed in the resulting ARPES spectrum.The method is not limited to simple photoelectron spectroscopy, but can be used to sim-ulate complex experimental techniques, such as the reconstruction of attosecond beating byinterference of two-photon transitions (RABBITT). Furthermore, due to the versatility ofthe real-time description in Octopus, the underlying excitation does not need to be coming42rom an optical pulse (i.e., an external vector potential) and, since the ions can move ac-cording to the TDDFT-Ehrenfest dynamics, one can also simulate features coming fromthe ionic or lattice motion, specifically electron-phonon coupling signature in ARPES.
Pump-probe simulations are not limited to dynamical process of excitations, but can also beemployed to study steady-state modifications of driven electronic states.
This opens thepossibility of studying Floquet physics from an ab-initio perspective and to directly simulatethe effect that a periodic force would have on the dressed electronic structure – an impor-tant aspect to underline given the growing interest in the field of Floquet engineering andFloquet analysis.
Finally, t-SURFF is particularly suited to simulating strong field ionization processessuch as in laser-induced photoelectron diffraction (LIED) experiments, where the ionizationtakes place by direct tunnelling into the continuum and the field is so strong as to driveelectrons in trajectories recolliding with the parent system. In fact, since in Eq. (48) weaccumulate the flux through S over the time of the propagation, with t-SURFF electronscan seamlessly flow back and forth through the surface without producing any artifact inthe final spectrum. As an example, in Fig 11 (d) we present the photoelectron angulardistribution of C ionized by a strong IR pulse capable of inducing rescattering dynamics,which then gets imprinted in the photoelectron spectrum as characteristic rings centered inthe value of the vector potential at the instant of rescattering. In this regime, the result ofsimulations obtained with Octopus are in excellent agreement with the experiments.
XI. ELECTRIC AND THERMAL CONDUCTIVITIES
Electrical conductivity in real materials can be described to a first approximation byOhm’s Law, which is a linear relationship between the applied electric field and the cur-rent density generated in response. The standard method to study electrical and thermalconductivity is to carry out a DFT calculation at equilibrium and apply Kubo-Greenwoodlinear response theory to evaluate the conductivities. Although this method has been usedsuccessfully in the past to calculate transport properties, the application is limited tosystems in the linear response regime. It is known that in many cases, especially in thepresence of strong external fields, Ohm’s law is no longer valid and the system exhibitsnon-linear behavior. In order to capture these complex behaviors in materials, one must43o beyond simple, linear approximations and describe electron interactions in a non-trivialmanner.One promising route to go beyond linear response is to use density functional methodsto directly study thermoelectric transport. This topic has received less attention, but hasrecently been implemented in Octopus and applied to liquid aluminum . For detaileddescriptions and derivations of non-equilibrium thermoelectric phenomena using densitybased methods, we refer the reader to Ref. 147.In Octopus, we are able to calculate the current density and heat current density at eachstep during a time-dependent simulation. A current is induced in this study by applying anelectric field of the form E ( t ) = E δ ( t ), where E describes the magnitude and direction ofthe field. The electric field is induced through a time-dependent vector potential ( E ( t ) = − c − ∂ t A ( t ), c is the speed of light) in the Hamiltonian. This vector potential satisfies theperiodic boundary conditions of an extended system. It should be noted that there is nolimitation on the form of the applied electric field in general. We are able to evaluate themacroscopic current density as J ( t ) = − i Ω (cid:90) d r N (cid:88) j ϕ i ( r , t )[ ˆ H ( t ) , ˆ r ] ϕ i ( r , t ) , (50)where Ω is the volume of the unit cell. The energy current density, ˆ J h ( r , t) is expressed asˆ J h ( r , t ) = ˆ J t ( r , t ) + ˆ J v ( r , t ) + ˆ J u ( r , t ) + ˆ J f ( r , t ) , (51)where ˆ J t ( r ) is the kinetic energy contribution given as[ˆ J t ( r , t )] i = i (cid:16) [ ∂ i ˆ ϕ † ][ ∇ ˆ ϕ ] − [ ∇ ˆ ϕ † ][ ∂ i ˆ ϕ ] − [ ∂ i ∇ ϕ † ] · [ ∇ ϕ ] − [ ∇ ˆ ϕ † ] · [ ∂ i ∇ ˆ ϕ ] (cid:17) (52)and ϕ = ϕ ( r , t ) are understood to be the time-dependent states. The output of the time-dependent simulation can be further analyzed to yield the frequency-dependent conductivity.The electrical (or thermal) conductivity σ can be found by Fourier transforming the corre-sponding current as σ ij ( ω ) = 1[ E ] i (cid:90) ∞ dte − iωt J j ( t ) . (53)In general, it is possible to study the conductivity of an extended system in the presenceof an applied electric field. Here, we illustrate the use of this method on hydrogen at 1400 Kand 400 GPa. At this temperature and pressure, hydrogen is in its liquid metallic phase. A44 C u rr e n t( a . u . ) time (fs)CurrentHeat Current 00.010.020.030.040.050.060.07 0 0.2 0.4 0.6 0.8 1 C o ndu c t i v i t y ( a . u . ) frequency (a.u.)ConductivityHeat Conductivity Figure 12. (Left) Time-dependent current density for one ionic snapshot of a supercell of 128atoms of hydrogen in the liquid phase at 1400 K and 400 GPa with an initial electric field E of0.1 a.u. The blue line represents the current density as given in Eq. (50) and the red line is the heatcurrent evaluated from Eq. (52). Both currents are shown decaying to zero. (Right) The Fouriertransform of the respective currents gives the frequency-dependent conductivities. molecular dynamics simulation was carried out in VASP to produce several ionic snapshotsof the system. For each ionic snapshot, we performed a TDDFT simulation with a timestep of 0.05 a.u. (0.00121 fs) for a total time of 3 fs. The initial electric field E was 0.1 a.u.The calculation was carried out on a 3x3x3 k-point grid to ensure that the current densitydecays to zero.In Fig. 12, we have plotted the time-dependent current density and heat current densityas a function of time after an initial electric field is applied at time zero. Both currentsare shown to decay to zero by the end of the simulation. It should be noted that the heatcurrent corresponds to the kinetic energy contribution given in Eq. (52). In general, thisheat current is related to the Peltier coefficient, since it describes the heat generated fromthe response to an electric current. The remaining contributions in Eq. (51) are not yetimplemented in Octopus. Both J u and J f arise from electron-electron interactions and J v is the potential energy. In future work it would be interesting to evaluate these remainingcontributions to determine if they are large or can reasonably be ignored for some systems.The electrical conductivity and the heat conductivity can be evaluated using Eq. (53). Thesefrequency-dependent conductivities are shown in the right hand panel of Fig. 12. One may45lso calculate the DC conductivity by taking ω = 0.This suite of tools implemented in Octopus will allow for the study of time-dependentthermoelectric phenomena using TDDFT. This approach to study current and conductivityis more general and widely applicable than the standard Kubo-Greenwood approach. WhileKubo-Greenwood is a linear response theory, this method can be applied to study materialswhere non-linear conduction effects can appear. In a recent study, this TDDFT method wasbeen used to illustrate non-linear conductivity effects in liquid aluminum, which wouldnot be possible by applying standard linear response theory. In the future, we plan to extendthese tools by implementing a thermal vector potential that can induce a heat current in anextended system during a time-dependent simulation. XII. LOCAL DOMAIN CONTRIBUTION TO PHYSICAL OBSERVABLES
Usually, we are interested in studying systems consisting of different atoms, molecules,regions, or domains. In such cases, we might want to understand the different contributionsfrom different parts of the system towards a specific observable. Based on the Hohengberg-Kohn theorems for the DFT, and its extension, the Runge-Gross theorem, for TDDFT,we can state that any physical observable of the system is a functional of the electronicdensity, either static (ground state) or time dependent. In other words, the expectationvalue of an operator ˆ O can be expressed as functional of the electronic density O [ n ] = (cid:104) Ψ[ n ] | ˆ O| Ψ[ n ] (cid:105) . (54)Let’s now consider the case where we write the total electronic density as a sum of N densities n ( r ) = N (cid:88) i n i ( r ) , (55)and no further constraints are imposed on n nor on n i . Such partitioning of the density isat the basis of subsystem DFT, which allows to divide a system into several Kohn-Shamsubsystems that interact with each other in a theoretically well justified manner (see Ref. 150and references therein). However, in our case, we are only interested in this kind of partitionfor post-processing purposes where we assign a density n i to a specific part of the systemin order to identify how it contributes to a given observable. We are thus interested in46perators for which the following condition is true O [ n ] = N (cid:88) i O [ n i ] . (56)Such operators are said to be additive. One example of an additive operator that is partic-ularly relevant in the context of TDDFT is the time-dependent dipole d ( t ) = (cid:90) r n ( r , t ) d r = N (cid:88) i (cid:90) r n i ( r , t ) d r , (57)as this is the relevant observable for obtaining the optical absorption cross section of finitesystems.The range of operators that fulfill Eq. (56) can be further expanded if we now considernon-overlapping densities, that is, densities n i that are non-zero only in a given domain V i and that the domains do not overlap. We refer to this type of partitioning as a localdomain partitioning. In such case, any (semi-)local operator that depends (semi-)locally onthe density should fulfill Eq. (56). In this context, an example of a relevant observable isthe exchange-correlation energy within the LDA approximation E LDAxc [ n ] = (cid:90) n ( r ) e LDAxc ( n ( r )) d r = N (cid:88) i (cid:90) V i n i ( r ) e LDAxc ( n i ( r )) d r , (58)where e LDAxc is the exchange-correlation energy per unit particle.Currently, several different options can be used in Octopus to define specific regions ofthe simulation box. These options include simple geometric shapes, such as spheres centeredon particular points of space, atom-dependent domains, such as an union of spheres centredon the atoms, or the definition of Bader volumes. The latter option follows the QuantumTheory of Atoms in Molecules (QTAIM) and associates a density region to a given atomthrough a density gradient path, such that the boundaries of each volume are defined as thesurfaces through which the charge density gradient has a zero flux. Note that some of theseoptions allow for the user to specify overlapping regions. In that case the overlap must betaken into account when analyzing the results and extra care is required when comparingresults from different domains.Let us exemplify the applicability of the local domain partitioning by decomposing theoptical response for a coupled chromophore system. Similar to Frozen Density Embedding47eal-time TDDFT (FDE-rt-TDDFT), we decompose the total optical spectrum as a sumof the local response within each domain α ( ω ) = (cid:88) i α i ( ω ) , (59)where α is the dynamic polarizability. In addition, from each local dynamic polarizabiltytensor we can compute for each fragment the corresponding cross-section. This is justifiedby the fact that the relevant operator to calculate α is the dipole operator and, as shownabove, this is an additive operator that fulfills Eq. (56).As an example, we performed a series of optical spectrum simulations using real-timeTDDFT for a benzene-fulvene dimer with different π -stacking separation, ranging from 4 ˚Ato 8 ˚A. We use the standard PBE exchange-correlation functional and the OptimizedNorm-Conserving Vanderbilt PBE pseudopotential (sg15) set. The real-space grid is de-fined as a parallelepiped box with length 16 ˚A, 17 ˚A and 20 ˚A in the three Cartesian axes.The spacing between points is set to be 0.13 ˚A.For each intermolecular separation we carry out a single ground-state calculation andthree time-propagations. For each time-propagation a dipolar electric perturbation is appliedalong one of the Cartesian axes.
We let the perturbed Kohn-Sham states evolve for apropagation time of T = 24 (cid:126) /eV (15.8 fs), with resolution of 0.26 eV (2 π/T ). The ground-state electron density is fragmented following the Bader atomic decomposition and eachmolecular domain is defined as the sum of these atomic volumes. Then, the correspondingdipole operator is applied over each defined domain. Finally, by Fourier transform of thelocal time-dependent dipole moment, the local polarizability tensor is recovered.Figure 13 shows the spectral decomposition of the benzene-fulvene system as a functionof the intermolecular separation. The full system at large intermolecular separation (8 ˚A)shows two major peaks located at the same frequencies as for the isolated molecules. As theintermolecular distance becomes smaller, the global spectrum changes due to electrostaticeffects induced by the electronic cloud of the neighbouring molecule. We can see that thepeak located at around 5 eV suffers a red-shift, while the stronger peak between 6.5 and7 eV reduces its intensity, giving rise to an oscillator strength transfer to excited states athigher energy. The local domain analysis reveals that the major change is caused by thereduction of the oscillator strength of the fulvene peak located at around 6.5 eV.These results are in agreement with the FDE-rt-TDDFT calculations of Ref. 152, further48 .01.02.03.04.0 4 5 6 7 8 9 A b s o r p t i on [ Å ] ω [eV] FulveneBenzene4Å5Å6Å8Å A b s o r p t i on [ Å ] ω [eV] Fulvene4Å5Å6Å8Å A b s o r p t i on [ Å ] ω [eV] Benzene4Å5Å6Å8Å
Figure 13. Top: Schematic representation of the benzene-fulvene dimer with π -stacking along the z axis and a intermolecular distance of 4 ˚A. From left to right, the electronic density for the globalsystem for an isosurface of 0.006 ˚A -3 and the QTAIM local density for the fulvene and benzenemolecules respectively. Visualized with USCF Chimera software.
Bottom: Photo-absorptioncross section obtained from, from left to right, the full system density, the fulvene local density,and the benzene local density. validating our strategy. In addition, the local domains methodology does not require anyprevious fragmentation, selection, or localization of the basis set, allowing also to treat verylarge systems, such as biological molecules like the major light harvesting complex LHC-II.
More recently, we have shown that this technique also allows for exciton coupling calculationswhen combined with a new formulation of the transition densities from real-time TDDFT.
XIII. NEW PROPAGATORS FOR REAL-TIME TDDFT
At the core of most Octopus calculations lies the need to integrate the time-dependentKohn-Sham equations (TDKS). These are a set of non-linear equations, given the dependenceof the Kohn-Sham Hamiltonian with the electronic density. Upon the discretization of theelectronic Hilbert space, they take the generic form of a first-order ordinary differentialequation (ODE) system ˙ ϕ = f ( ϕ ( t ) , t ) , (60)49 ( t ) = ϕ , (61)where t is the initial time, ϕ is an array containing all the Kohn-Sham orbitals, { ϕ Nm } ( ϕ is its initial value), and f is a vector function ( f = ( f , . . . , f N )) given by the action of theKohn-Sham Hamiltonian ˆ H [ n ( t ) , t ] f i ( ϕ ( t ) , t ) = − i ˆ H [ n ( t ) , t ] ϕ i ( t ) , ( i = 1 , . . . , N ) (62)Note that: (1) If the nuclei are also to be propagated, the state should be supplemented withtheir position and momenta, and the equations with the corresponding nuclear equationsof motion. (2) Likewise, the formalism for solids includes a polarization vector field thatmust also be included in the system definition and propagation. (3) Strictly speaking, theTDKS equations are not ordinary but “delay differential equations” (i.e. equations for whichthe derivative of the unknown function at a certain time depends on the function values atprevious times), due to the dependence of the exact exchange and correlation potentialfunctional on past densities. This is ignored in the adiabatic approximation, which is almostalways assumed.A myriad of numerical propagators for ODEs are available, all of them theoreticallyapplicable to the TDKS equations. Ideally, however, one should choose a propagator thatrespects all the mathematical properties of the equations that describe the problem at hand,for example, the preservation of the norm of the orbitals. But also, in the case of the TDKSequations in the adiabatic approximation, another such property is symplecticity (a fact thatis demonstrated in Ref. 159). This property is usually stated in terms of the conservationof the volume of the flow of the system of ODEs in the phase space, although the precisedefinition is as follows:Any system of ODEs can be viewed as a “flow”, a differentiable map g : R p → R p thattransforms the state y ( t ) at some point in time t into the state at time t , y ( t ) = g ( y ( t )).For systems described with complex variables such as ours, since it can be split into a realand an imaginary part, p is even: g : R N → R N . A map with an even number of variablessuch as this one is defined to be symplectic if and only if ∂g∂x T J ∂g∂x = J , J = I − I , (63)where I is the unit matrix of dimension N , and x ∈ R N (conventionally, the first N variablesof x are called the “coordinates”, and the second half are the “momenta”, but no physical50eaning should be assumed for them in this purely mathematical definition). A numericalpropagator is also a differentiable map that relates the solution y ( t ) to y ( t + ∆ t ). As such,it may respect or not the symplecticity – and other properties – of the original flow.The symplecticity also means that the system has to be “Hamiltonian”: one may find aset of coordinates q ∈ R N , p ∈ R N (in this case, the real and imaginary part of the orbitalscoefficients), and some scalar function H ( q, p ) (the Hamiltonian expectation value), thatpermits us to rewrite the system into the well known form of Hamilton’s equations:˙ q i = ∂H ( p, q ) ∂p i , (64a)˙ p i = − ∂H ( p, q ) ∂q i . (64b)The relevance of the symplectic numerical propagators stems from the fact that they presentsome features, such as a better conservation of the energy at long times, where the valueof the energy ends up oscillating around the true value instead of diverging. For a moredetailed discussion about symplecticity on numerical propagators, check Ref. 160.The Octopus code has several propagator options, many of them already described inRef. 161, such as the Crank-Nicolson, standard Runge-Kutta, exponential midpoint ruleand variations (e.g., the “enforced time-reversal symmetry” scheme), split operator tech-niques, etc. In a recent paper, some of the present authors have studied propagators ofdifferent families that had been scarcely (or not at all) tested for the TDKS equations: mul-tistep, exponential Runge-Kutta, and commutator-free Magnus (CFM) expansions. Afterconsidering both the accuracy, stability, and the performance of the propagators, the CFMtechniques were identified as suitable schemes for TDDFT problems and implemented inOctopus. In this section, we make a brief description.Developed in Ref. 162 for linear non-autonomous systems, the CFM expansion offer analternative to the “standard” Magnus expansion, which requires expensive application ofnested commutators of the Hamiltonian with itself at different times. In essence, a CFMexpansions Γ( t + ∆ t, t ) consists of substituting the propagator ˆ U ( t + ∆ t, t ) by a product ofexponentials Γ( t + ∆ t, t ) = m (cid:89) i =1 exp( ˆ D i ) . (65)The m linear operators ˆ D i are either the Hamiltonian at different times within the interval[ t, t + ∆ t ], or parts of it. 51e have implemented in Octopus an order four ( q =4) version, given for example inEq. 43 of Ref. 162, which we call hereafter CFM4. This propagator requires two exponentials( m = 2), ϕ ( t + ∆ t ) = exp (cid:16) − i ∆ t ( α ˆ H [ t ] + α ˆ H [ t ]) (cid:17) exp (cid:16) − i ∆ t ( α ˆ H [ t ] + α ˆ H [ t ]) (cid:17) ϕ ( t ) , (66)where α = 3 − √ , t = t + (cid:32) − √ (cid:33) ∆ t , (67) α = 3 + 2 √ , t = t + (cid:32)
12 + √ (cid:33) ∆ t . (68)ˆ H [ t ] and ˆ H [ t ] are the Hamiltonians at times t and t , which are in fact unknown becausethey depend on the Kohn-Sham states through the density: we are dealing with a non-linearproblem and the CFM expansions were in fact developed for linear systems. We have variousoptions to extend them for our non-linear problem: for example, one could define ˆ H [ t ] andˆ H [ t ] as interpolated Hamiltonians from ˆ H [ t ] and ˆ H [ t + ∆ t ], in which case we end up withan implicit equation for ϕ ( t + ∆ t ) that we would have to solve at a substantial cost.The alternative that we have implemented, however, is to approximate ˆ H [ t i ] via an extrap-olation from the Hamiltonian at various previous time steps (in practice, it is the Hartree,exchange, and correlation parts that must be extrapolated). The resulting method is thenexplicit, i.e., no linear or non-linear algebraic equations need to be solved. The fourth orderaccuracy is preserved as long as the extrapolation is also done at order four.As an example, Fig. 14 shows the performance of this CFM4 method against the wellknown exponential midpoint rule (EMR). The benchmark system consists of a benzenemolecule. It is subject to an instantaneous perturbation at the beginning of the propagation,and then it evolves freely for some fixed interval of time. The propagations are performedat varying values of ∆ t . The plot displays the cost of the propagation vs. the accuracyachieved in each case. As we can see, the CFM4 outperforms the EMR for all the valuesof the error examined, although the differences become more obvious when higher accuracyis demanded. The key difference between the two methods is the fourth-order accuracy ofCFM4 scheme – at the cost of requiring two exponentials, instead of just one, while theEMR is only second-order accurate. This is reflected in the different slopes of the curves asthe error becomes smaller, i.e., as ∆ t →
0. 52 l og ( C o s t ) log (Error)CFM4EMR Figure 14. Wall-time computational cost of the method (seconds) against the error in the wavefunction (defined as Error( T, ∆ t ) = (cid:113)(cid:80) m || ϕ m ( T ) − ϕ exactm ( T ) || , where the ϕ exactm are the values ofthe Kohn-Sham orbitals computed using the explicit fourth-order Runge-Kutta integrator with avery small time-step) for the EMR and CFM4 propagators. XIV. CONJUGATE GRADIENT IMPLEMENTATION IN RDMFT
The RDMFT implementation in Octopus has been described in detail in a previouspaper. In this Section we briefly review the existing implementation, providing some newinsights, and introduce a recently implemented method to solve the RDMFT equations thatis better suited to the real-space grids used in the code.The optimization of the natural orbitals in RDMFT is subject to an orthonormalizationconstraint for the orbitals. One minimizes the functional E total − M (cid:88) j,k =1 λ jk (cid:18)(cid:90) d r ψ ∗ j ( r ) ψ k ( r ) − δ jk (cid:19) , (69)where E total and M denote the total energy and the number of natural orbitals ψ j included inthe calculation, respectively. λ jk are the Lagrange multipliers which ensure the orthonormal-ity of the natural orbitals at the solution point. We note that, for a RDMFT calculation, thenumber of natural orbitals has to be larger than the number of electrons in the system (coreelectrons which are included in the pseudopotential do not count). The exact number M is53ystem-dependent and should be treated as an additional parameter with respect to whicha convergence study should be carried out, just like it is done for the basis set or the gridparameters. Typically, the optimization is performed using the so-called Piris method, which was the method previously implemented in Octopus. Within this approach one usesthe orthonormality constraint of the natural orbitals, which implies that a certain matrixconstructed from the Lagrange multipliers λ jk is diagonal at the solution point. As an im-mediate consequence of the Piris method, the natural orbitals at the solution point are linearcombinations of the orbitals used as starting point for the minimization. In other words, theinitial orbitals serve as a basis. Consequently, the necessary matrix elements for differentenergy contributions can be calculated for the basis functions (initial orbitals) before theiterative optimization of the natural orbitals and their occupation numbers is started. Inaddition, the optimization of the occupation numbers can be turned off completely, resultingin a Hartree-Fock calculation in the basis of the initial orbitals.For the existing implementation in Octopus of the Piris method, the initial orbitals aretaken to be the solutions obtained with a different level of theory, like independent particlesor DFT. In order to better understand the effect of the choice of basis, we have tested thefollowing choices: (i) independent particles, density functional theory within (ii) the localdensity approximation (LDA) or (iii) the exact exchange (EXX) approximation, as wellas (iv) the Hartree-Fock approximation. In all cases we have to ensure that the numberof unoccupied states in the calculation is sufficient to cover all the natural orbitals whichwill obtain significant occupation in the the following RDMFT calculation. The results forthe convergence of the total energy of a one-dimensional (1D) hydrogen molecule using theM¨uller functional are given in Fig. 15. The calculations were performed on a 1D gridextending from − . . .
03 bohr. The nuclear potentialfor the 1D molecule reads v ( x ) = − (cid:112) ( x − d ) + 1 − (cid:112) ( x + d ) + 1 , (70)with d = 1 .
628 bohr, which corresponds to the equilibrium geometry. The electron-electroninteraction in one dimension is described by the soft-Coulomb interaction w ( x, x (cid:48) ) = 1 (cid:112) ( x − x (cid:48) ) + 1 . (71)Since Octopus performs all calculations on a finite grid, we typically obtain a finitenumber of bound states in the calculations for the basis set and any additional orbitals54 igure 15. Total energy for the RDMFT calculation for one-dimensional H using the Pirismethod with different basis sets and using the conjugate gradient implementation. The insetshows a zoom into the area where convergence is reached. We employ the M¨uller functional forall calculations. extend over the whole grid, i.e., they are unbound and therefore delocalized. However, allnatural orbitals with non-zero occupation, because they decay with the ionization potentialof the system, are localized on the system. Hence, the extended basis states will onlycontribute with very small coefficients, if at all, and their inclusion in the basis set does notlead to a significant improvement of the results. In the past, this problem was addressedby performing an additional step to localize the initial states before starting the RDMFTcalculation. However, further investigations showed that this only improves the results for asmall number of natural orbitals in the calculation. Testing the convergence with respect tothe number of natural orbitals then showed that the additional localization step slows downthe convergence with respect to the number of basis functions. The fastest convergenceand lowest total energies are obtained by using the results from an independent particlescalculation, as those yield the largest number of localized orbitals from all the different basissets that were tested, as shown in Fig. 15. As the additional localization step proved to beunnecessary and even hinders convergence, it has been removed from the new version of thecode.Since the natural way of representing quantities in Octopus is directly on the real-space55rid, and to circumvent the limitations with the quality of the basis sets available for thePiris method, we have decided to implement a conjugate gradient optimization of the naturalorbitals. This implementation follows the procedure for DFT explained in Ref. 167, whichwe adapted for RDMFT. This procedure allows us to take advantage of the full flexibility of areal-space grid and provides a systematic way of improving the results by enlarging the gridand reducing the spacing between grid points. The conjugate gradient algorithm requiresa set of initial orbitals to start the self-consistent calculation, however, at convergence theresults are independent of that starting point. Therefore, while the calculation using thePiris method requires a set of initial states which serve as the basis, the conjugate gradientalgorithm can be used starting from a initial set of random states. In our tests of theconjugate gradient implementation, the quality of the initial states only had an influenceon the number of iterations necessary for the convergence, but not on the final result. Wesuggest to use the orbitals obtained from an independent particles calculation as initial statessince they can be obtained for small numerical costs and simultaneously can serve as a basisset in the Piris implementation.Since we are not solving an eigenvalue equation to obtain the natural orbitals, they arenot automatically orthogonal. As mentioned above, the orthogonality is taken into accountvia a constraint. Compared to Ref. 167, the non-diagonal Lagrange multipliers λ jk lead totwo modifications in the conjugate gradient procedure. First, the steepest-descent direction(cf. Eq. (5.10) of Ref. 167) reads here ζ mi = − ˆ Hψ mi + M (cid:88) k =1 λ ik ψ mk , (72)where λ ik = (cid:104) ψ mi | ˆ Hψ mk (cid:105) . (73)Second, one parameter of the line-minimization (cf. Eq. (5.26) of Ref. 167) is changed to ∂E total ∂ Θ = (cid:104) φ (cid:48) im | ˆ Hψ mi (cid:105) + (cid:104) ψ mi | ˆ Hφ (cid:48) im (cid:105) − (cid:88) k ( λ ki (cid:104) φ (cid:48) im | ψ mi (cid:105) + λ ik (cid:104) ψ mi | φ (cid:48) im (cid:105) ) , (74)where Θ parameterizes the descend in the direction of the gradient and the single-particleHamiltonian acting on a state | φ (cid:105) is defined asˆ H | φ (cid:105) = ∂E total ∂φ ∗ . (75)56or details of the notation, we refer the reader to the paper by Payne et al . The convergence study with respect to the number of natural orbitals for the conjugategradient algorithm is also included in Fig. 15. As one can see, a smaller number of naturalorbitals needs to be included in the calculation than for the basis set implementations withthe Piris method. As the number of natural orbitals equals the number of basis functionsin these calculations, this is mostly due to the fact that the available basis sets are ofrather poor quality. In addition, the converged total energy for the conjugate gradientalgorithm is slightly lower than for all the basis sets using the Piris method (see inset ofFig. 15), which shows that the conjugate gradient algorithm exploits the full flexibility ofthe grid implementation, allowing for contributions to the natural orbitals which are notcovered by any of the basis sets. We have also verified that the converged result of theconjugate gradient method is indeed independent of the choice of initial state for the M¨ullerfunctional employed in our calculations. The M¨uller functional is known to be convexwhen all infinitely many natural orbitals are included.
This property is most likely notshared by all available RDMFT functionals. In addition, the number of natural orbitals inany practical calculation is always finite. Consequently, in practice one needs to test theconvergence for the different starting points, as the appearance of local minima cannot beexcluded.
XV. PERIODIC SYSTEMS AND SYMMETRIES
Electronic structure in periodic systems is usually described using plane waves, but real-space grids have been shown to be a viable alternative when performing DFT and TDDFTcalculations.
Unfortunately, the discretization introduced by the real-space grid oftenbreaks the direct connection between the physical system and the basis set, as symmetriesand translations are, in most of the cases, not compatible with the discretized grid. However,real-space grids offer many advantages, such as natural mixed periodic-boundary conditionsfor semi-periodic systems, and the calculation of the exchange and correlation term of DFTis straightforwardly obtained on the grid.One of the main challenges for treating periodic systems in real space is the generationof the real-space grid, and the corresponding weights for the finite differences. This be-comes relevant when the grid is generated along the primitive axes of the Wigner-Seitz cell57primitive cell) of a solid, where the generating axes are usually non-orthogonal.In Octopus, the grid points are generated along the primitive axes, and the calculation ofthe finite-difference weights follows the implementation described in Ref. 170. Once the gridis generated, and the weights for the finite differences (gradient and Laplacian) are obtained,it is necessary to deal with the discretization of the reciprocal space. Indeed, for periodicsystems, the full crystal is replaced by the primitive cell of the crystal in real-space, thanksto the Bloch theorem, but is then complemented by the Brillouin zone, which also needs tobe sampled. This grid, usually called k -point grid, is common to any type of basis sets, aslong as one decides to reduce the crystal to its primitive cell.In order to reduce the numerical effort associated with the description of periodic systems,we make use of the symmetries to reduce the Brillouin zone to its irreducible Brillouin zone,which can drastically reduce the number of k -points. The space group of the crystal, as wellas its symmetries, are obtained thanks to the spglib library. In the present implementa-tion, we restrict the symmetries to the symmorphic symmetries (inversion, rotations andmirror planes), leaving for later the so-called non-symmorphic symmetries, i.e., symmetriesinvolving a fractional translation, as they are not compatible with arbitrary real-space grids.In order to assess the validity of our implementation, we used the so-called Delta-factortest.
Using the Schlipf-Gygi ONCVPSP 2015 pseudopotential set, we obtained the valueof 1 .
50 meV/atom, which is very close to the value obtained by other codes using the samepseudopotential set.When investigating the electron dynamics driven by a laser field, or any type of symmetry-breaking perturbation (vector potentials, kicks with a finite momentum, strain, etc), someof the original symmetries are lost. To deal with that, Octopus finds the small groupof symmetries corresponding to the original space-group of the solid, retaining only thesymmetries that leave invariant the perturbation direction. This defines the symmetriesthat are used for a time-dependent calculation.One important aspect when using symmetries, is the symmetrization (in real-space) ofthe charge and current densities, as well as other observables, such as kinetic energy densityfor instance. In complement to the reduction of the Brillouin zone to its irreducible Brillouinzone, we implemented a real-space symmetrization. We found that performing the real-spacesymmetrization is very important to achieve good and stable numerical results when takinginto account symmetries. 58or a comprehensive description of periodic systems, the ion dynamics has to be consid-ered in addition to the electron dynamics. For isolated systems, it is often described with theTDDFT-Ehrenfest dynamics method, where the ions obey the Newton equation withforce fields computed by the TDDFT. In contrast, for periodic systems, the ion dynamicshas to be described with two sets of equations: One is the Newton equation for ions on thereduced coordinates in the primitive cell, and the other is the equation of motion for theprimitive cell itself, as the ionic coordinates of periodic systems are described by the combi-nation of the reduced coordinates with the lattice vectors of the primitive cell. The latticedynamics is often treated with the Parrinello-Rahman method, and the key ingredientof the equation of motion is the stress tensor. Therefore, to realize the ab initio simulationsfor electron-ion-lattice dynamics based on the TDDFT, the calculation of the stress-tensorhas been implemented based on Ref. 176, and the implementation of the lattice dynamics isunder way.
XVI. ADDITIONAL TECHNICAL CODE IMPROVEMENTS
In order to make all the new developments and applications described in the above sec-tions possible, the code needs to be accurate, efficient, and reliable. Also, when a codereaches the size and complexity of Octopus (currently more than 200,000 lines of sourcecode), the amount of time required for maintenance and for adding new features becomesconsiderable. Therefore, continuous efforts at optimizing, validating, and improving thesource code quality are needed. These efforts are essential, but often do not get the atten-tion they deserve. In this section we present some noteworthy developments to improve thecode reliability and efficiency, to make the project easier to maintain, and to tackle newcomputer architecture challenges.
A. Web application for analyzing regression tests
In a complex code like Octopus, every new development may have unintended side effects,possibly introducing errors to already existing features. To avoid such a regression, everychange to the code is required to pass a suite of tests before being accepted. This suite isexecuted automatically using Buildbot and is integrated with the continuous integration59ramework of the gitlab platform where Octopus is hosted for several years now. The testsuiteis executed with 30 different toolchains spanning a variety of architectures, compilers, andMPI libraries. It contains currently about 180 tests that execute Octopus roughly 740 timesand that contain about 12000 comparisons to reference values. The overall coverage of thetests is determined using the codecov.io service and lies at about 71% at the moment.Although the testsuite is efficient in avoiding regressions, analyzing why a test failedhas been difficult so far due to the large amount of data generated. Thus, we have imple-mented an interactive analysis and visualization of the testsuite results as a web applicationavailable at https://octopus-code.org/testsuite/ . This application makes it easier tounderstand why a test failed: is it a problem with just one toolchain, e.g. , a particular ar-chitecture, or is it simply a larger numerical variation of the results? Moreover, it facilitatesupdating the tests and improving the testsuite itself.The application has been implemented using the web framework Django coupled to aPostgres database. The testsuite results are automatically uploaded by the Buildbot serviceas soon as they are available. The application allows to analyze single toolchain runs andsingle comparison matches and to compare all toolchain runs for a single commit in git. Itprovides histograms to judge if there are outliers or if it is a broad distribution; moreoverthis allows to determine, e.g. , if there is a difference between MPI and non-MPI toolchains.This application has already helped to identify bugs causing regressions and will continueto be useful to understand failed tests, update them, and to improve the testsuite.
B. Improving ground-state calculations
To improve the reliability of ground-state calculations with Octopus, the default eigen-solver used inside the self-consistent field (SCF) cycle has been improved and differentreal-space preconditioners for the eigensolvers have been evaluated.The default eigensolver, a conjugate-gradients algorithm, has been improved over theprevious implementation by now following closely Ref. 167, which greatly improves thereliability of the solver. The updated implementation differs from Ref. 167 in only onepoint: by default, the current band is not orthogonalized against all bands, instead it is onlyorthogonalized against previously computed bands with lower energies. According to ourtests, in Octopus, this leads to a faster convergence for most cases.60reconditioners for eigensolvers are an integral part of SCF calculations because theygreatly accelerate convergence. They achieve this by applying an approximate inverse ofthe Hamiltonian in each iteration that leaves the solution invariant and brings the systemcloser to the real solution. As Octopus uses a real-space grid (as opposed to plane waveslike in Ref. 167), different preconditioners are needed. A range of preconditioners has beencompared for a wide variety of systems, comprising molecules, semiconductors, metallicsystems, and surfaces to test convergence for very different cases.The default preconditioner in Octopus is a low-pass filter obtained by adding the neigh-bouring values in each dimension to the current value at a grid point, weighted by a certainfactor αψ (cid:48) i,j,k = αψ i,j,k + (1 − α ) ( ψ i − ,j,k + ψ i +1 ,j,k + ψ i,j − ,k + ψ i,j +1 ,k + ψ i,j,k − + ψ i,j,k +1 ) . (76)This preconditioner was first described by Saad et al. with α = 0 .
5. It has proven to bethe most effective preconditioner in our comparison because it is quite cheap to apply and itnevertheless decreases the number of SCF iterations noticeably. It can also be understoodas two weighted Jacobi iterations (up to a prefactor) to solve for the inverse of the kineticterm of the Hamiltonian (i.e., 0 . ψ (cid:48) = ψ ). From the convergence radius of the Jacobiiterations, we can conclude that the allowed values for α are between 0.5 and 1. Fromtheoretical considerations of damping of different spatial wavelengths during Jacobi itera-tions (see Ref. 179, Secs. 4.1.3, 4.2), the ideal value should be 0.75 to most effectively damphigh spatial frequencies. In practice, we find that this depends on the system. Moreover,increasing the number of Jacobi iterations does not make the preconditioner more effective.Also, using a single Jacobi iteration (i.e., dividing by the diagonal of the Laplacian) doesnot speed up the convergence significantly.A multigrid method similar to the one implemented in GPAW also uses Jacobi iter-ations to solve for the kinetic Hamiltonian, but employs different grids to reach a fasterdecrease of the error. Although this method reduces the number of SCF iterations, it iscomputationally more expensive because the Laplacian is applied several times for eachiteration and thus less effective than the filter preconditioner.A preconditioner specifically targeting real-space methods was proposed by Seitsonen et al. It uses the ratio between the difference of the energy and the potential to thekinetic energy (their Eq. (3)) with a preconditioning function (their Eq. (4)) originally from61ef. 167. However, this preconditioner did not speed up convergence significantly when usedin Octopus.Another way to get an approximate inverse of the Hamiltonian is to solve the Poissonequation associated with the kinetic part of the Hamiltonian. This reduces the number ofiterations needed to converge the SCF calculation, but it increases the total calculation time,as each iteration is much more costly.In summary, we find that for our real-space code, the filter preconditioner is most effectivein reducing the total time needed to compute ground states because it reduces the totalnumber of iterations without being computationally too expensive.Nonetheless, sometimes certain states cannot be fully converged using a preconditioner,especially for calculations with many unoccupied states. In this case, restarting withoutpreconditioner is needed to obtain full SCF convergence. The reasons are still under inves-tigation.
C. Novel multi-system framework
In the way Octopus was originally designed, the entire code was structured around theidea that there was only one system (albeit of arbitrary size and complexity) with an associ-ated Hamiltonian. This was typically a combined system of electrons and nuclei, where theformer were described using DFT and the latter were treated classically. All possible sorts ofinteractions of this system with some external source (e.g., external electromagnetic fields,solvents described with the PCM, etc) were then added in an ad hoc fashion. Although thisapproach worked well for many features and applications, its limitations became evidentwith several recent developments, in particular the coupled Maxwell-Kohn-Sham equationsdescribed in Sec. II, where the Kohn-Sham orbitals need to be time-propagated alongsidethe Maxwell fields.This has prompted us to start a major refactoring of the code, where the full physicalsystem is treated as several subsystems interacting with each other. Such subsystems canbe electrons, nuclei, Maxwell fields, etc. While this framework was mainly motivated by theMaxwell-TDDFT coupling, it has been implemented very generically so that all the existingfeatures can be converted to it and that, in the future, many other developments whichrequire the coupling of several subsystems can be based on this.62 . Memory layout
A new memory layout for storing the orbitals was introduced in Ref. 182, where allstates are stored in a number of smaller batches with the innermost index being the stateindex instead of the real-space grid index. Thus, for all states in a batch, exactly the sameoperations can be executed when looping over the grid while accessing memory contiguously.This allows efficient parallelization on GPUs, where all threads in a warp need to executethe same instructions, as well as vectorization on CPUs, where one instruction can operateon several data points. To fully utilize these instructions, the kernel for computing finitedifferences (e.g., the Laplacian) has been specialized to explicitly use SSE, AVX, or AVX512instructions. Now, more parts of the code have been ported to use this new layout to increasetheir performance and, whenever possible, the code will use this layout by default.
E. GPUs
The first GPU implementation has been described in Ref. 182. This was based on OpenCLand was limited to a small selection of numerical algorithms and calculation modes. Theseincluded some of the most commonly used features of the code or algorithms that wereparticularly suited for GPU porting, like the RMM-DIIS eigensolver for ground-state calcu-lations or the enforced time-reversal-symmetry propagator for time-dependent calculations.Since then, the GPU implementation has been expanded in several different ways. Mostnotably, it now supports CUDA through an additional compatibility layer. The number ofsupported features and algorithms has also been increasing steadily, such that most time-dependent calculations can now be run efficiently on GPUs. The implementation has alsobeen expanded to support multiple devices per host. Using the packed storage format de-scribed above, Octopus is able to store the states fully in the GPU memory, provided thememory is large enough, thus reducing memory transfers to a minimum. Recently, this wasimproved further by removing frequent allocations and deallocations of temporary variableson the GPU, now using a custom memory management for those. This proved very effectivein improving the scaling to several GPUs, also to several nodes with GPUs.We show two examples of time-dependent runs in Fig. 16. For the first example (leftpanel), β -cyclodextrine was used as an input system and the simulation was run on the GPU63 Sp ee d - u p IdealSkylake + 2 V100 1 2 4 8GPUs1248 Sp ee d - u p Ideal8 V100 machine
Figure 16. Scaling plots for the GPU implementation of Octopus. The left panel shows speed-upfor a machine with two 20-core Intel Xeon 6148 Gold sockets (Skylake architecture) and two V100GPUs (PCIe) per node. The right panel shows speed-up for a Supermicro server with 8 V100GPUs (NVLink). island of the COBRA supercomputer at the Max Planck Computing and Data Facility. Eachnode in this island consists of two 20-core Intel Xeon 6148 Gold sockets (Skylake architecture)with two Nvidia Volta 100 GPUs (PCIe); the nodes are connected with an Omnipath fabric(100 Gbit/s). When executed on a full node with GPUs, the average time step for theexample is a factor of 4.8 faster than on one full CPU node, i.e., the time to solution isreduced by a factor of 4.8. This also means that the GPU version consumes less energy:although one GPU node draws more power than a CPU node (ca. 950 W vs 450 W), thefaster execution time reduces the overall energy to solution by a factor of 2.3 when runningon a GPU node. When scaling to 2, 4, and 8 nodes, the speed-up is 1.26, 2.14 and 3.37,respectively, and the corresponding parallel efficiency is 63%, 54%, and 42%. Althoughthe scaling is not yet perfect, it has improved considerably and also hints at the need forimproving inter-node communication. For the second example (Fig. 16, right panel), a time-dependent simulation of a part of a chlorophyll complex was executed on a Supermicroserver with two 8-core Intel Xeon 6134 Gold CPUs (Skylake architecture) and eight NvidiaVolta 100 GPUs interconnected with NVLink. Here, the speed-up of the average time stepof the simulation when using 2, 4, and 8 GPUs is 1.95, 3.70, and 5.97 which corresponds toparallel efficiencies of 97%, 93%, and 75%. These efficiencies are much better than for themulti-node example, probably because data is only communicated within one node. For bothexamples, all states were stored in the GPU memory, and the parallelization was achieved64y distributing the states only. Distributing the grid (i.e., parallelization over domain) isnot as effective, because it leads to more frequent communication, and thus has a largeroverhead.Our current and future efforts regarding the GPU implementation are focused on twoaspects. First, we are planning to enhance the existing implementation by improving thecommunication, possibly by overlapping computation and communication, and by optimiz-ing the existing kernels. Second, we want to port more algorithms that are currently onlyimplemented for the CPU version, especially for spin-dependent calculations and ground-state runs. This is a long-term effort and the ultimate goal is to have as many features andalgorithms as possible running efficiently on GPUs.
XVII. CONCLUSIONS
It has been almost 20 years since the development of Octopus started. During this periodthe code has matured and expanded to cover an ever-growing range of methods, theories,applications, and systems. It has also continuously adapted to the available computer archi-tectures and computing paradigms, allowing researchers to tackle increasingly challengingproblems in electronic structure theory.Although many of the code capabilities have been routinely used by many groups aroundthe world, mainly to study electronic excited states properties and dynamics, which stillremain to a large extent the core of Octopus, we believe the main reason for the code’s successlies elsewhere. From the beginning, the code has been designed to take full advantage of theflexibility and versatility offered by the use of real-space grids and provide developers with aframework to easily implement and test new ideas and methods that can be later on adoptedby other codes (as it has already been the case with quite a few features previously developedwithin Octopus). This is demonstrated by the number of new theoretical methodologies andframeworks presented in this paper to deal with non-equilibrium phenomena of complexsystems and to address the combined dynamics of electrons, phonons, and photons that gobeyond what can be found in other electronic structure codes. This includes several novelapproaches to treat the coupling of the electronic systems to the photons, like the coupledMaxwell-Kohn-Sham equations, the OEP approach to the electron-photon coupling, and thedressed RDMFT. Other examples include the description of magnons in real-time and the65rbital magneto-optical response in solids and molecules using the Sternheimer approach.In the case of magnons, supercells need to be employed, and the scalability of real-spacegrid methods make this approach very promising for future applications in more complexcorrelated magnetic systems. We expect that many of these methods and approaches willbe integrated into other electronic-structure codes and will become standard tools in thenear future.The code flexibility is also demonstrated by the variety of systems that it can efficientlytreat, as the applications showcased in this paper include molecules, nanoparticles, modelsystems, solvents, solids, monolayers, etc. Particularly noteworthy is the efficient treatmentof periodic systems, which traditionally have been described using plane-wave basis-sets,and this has been achieved without loss of accuracy, as shown by our results for the Delta-factor test. In the end, it is our hope that combining this flexibility with a growing rangeof methodologies will provide researchers with the necessary tools to study new challengingphenomena, like novel correlated materials, out of equilibrium physics, or coupled electron-boson systems.We have also discussed recent improvements in performance and scalability, with par-ticular emphasis on the GPU support. These developments are crucial in the view of thenew challenges that electronic-structure applications are facing with the upcoming exaflopsupercomputers.Finally, several other implementations are in the pipeline, such as Floquet and cavity QEDmaterials engineering, multitrajectory methods to deal with the nonadiabatic electron-iondynamics, treatement of open quantum dissipative systems, spectroscopies with entangledphotons, etc. All of these should be available to the users in the next years, so we invite youto stay tuned to the Octopus webpage at https://octopus-code.org/ . ACKNOWLEDGMENTS
The authors would like to thank all the people that have contributed to the develop-ment of Octopus over the last two decades. They would also like to thank Lin Lin foruseful and interesting discussions and acknowledge the open discussions about real spacemethods with the group of Prof. Chelikowsky. This work was supported by the EuropeanResearch Council (ERC-2015-AdG694097), the Cluster of Excellence “Advanced Imaging of66atter”(AIM), Grupos Consolidados (IT1249-19) and SFB925. The Flatiron Institute is adivision of the Simons Foundation. XA, AW and AC acknowledge that part of this workwas performed under the auspices of the U.S. Department of Energy at Lawrence LivermoreNational Laboratory under Contract DE-AC52-07A27344. JJS gratefully acknowledges thefunding from the European Union Horizon 2020 research and innovation program underthe Marie Sklodowska-Curie Grant Agreement No. 795246-StrongLights. JF acknowledgesfinancial support from the Deutsche Forschungsgemeinschaft (DFG ForschungsstipendiumFL 997/1-1). DS acknowledges University of California, Merced start-up funding.
REFERENCES M. A. L. Marques, A. Rubio, E. K. Gross, K. Burke, F. Nogueira, and C. A. Ullrich,Time-dependent density functional theory, Vol. 706 (Springer Science & Business Media,2006). M. A. L. Marques, N. T. Maitra, F. M. Nogueira, E. K. Gross, and A. Rubio,Fundamentals of time-dependent density functional theory, Vol. 837 (Springer Science &Business Media, 2012). T. W. Ebbesen, “Hybrid light–matter states in a molecular and material science perspec-tive,” Accounts of Chemical Research , 2403–2412 (2016). M. Ruggenthaler, J. Flick, C. Pellegrini, H. Appel, I. V. Tokatly, and A. Rubio,“Quantum-electrodynamical density-functional theory: Bridging quantum optics andelectronic-structure theory,” Phys. Rev. A , 012508 (2014). J. Flick, M. Ruggenthaler, H. Appel, and A. Rubio, “Kohn-sham approach to quantumelectrodynamical density-functional theory: Exact time-dependent effective potentials inreal space,” Proc. Natl. Acad. Sci. U. S. A. , 15285–15290 (2015). J. Flick, N. Rivera, and P. Narang, “Strong light-matter coupling in quantum chemistryand quantum photonics,” Nanophotonics , 1479–1501 (2018). J. Feist, J. Galego, and F. J. Garcia-Vidal, “Polaritonic chemistry with organicmolecules,” ACS Photonics , 205–216 (2017). R. F. Ribeiro, L. A. Mart´ınez-Mart´ınez, M. Du, J. Campos-Gonzalez-Angulo, andJ. Yuen-Zhou, “Polariton chemistry: controlling molecular dynamics with optical cav-ities,” Chemical science , 6325–6339 (2018).67 A. F. Kockum, A. Miranowicz, S. De Liberato, S. Savasta, and F. Nori, “Ultrastrongcoupling between light and matter,” Nature Reviews Physics , 19–40 (2019). I. V. Tokatly, “Time-dependent density functional theory for many-electron systems in-teracting with cavity photons,” Phys. Rev. Lett. , 233001 (2013). J. Flick, Exact nonadiabatic many-body dynamics: Electron-phonon coupling in photoelectron spectroscopy and light-matter interactions in quantum electrodynamical density-functional theory,Ph.D. thesis, Humboldt-Universit¨at zu Berlin Berlin (2016). M. Ruggenthaler, N. Tancogne-Dejean, J. Flick, H. Appel, and A. Rubio, “From aquantum-electrodynamical light–matter description to novel spectroscopies,” Nature Re-views Chemistry , 0118 (2018). J. Flick, C. Sch¨afer, M. Ruggenthaler, H. Appel, and A. Rubio, “Ab initio optimizedeffective potentials for real molecules in optical cavities: Photon contributions to themolecular ground state,” ACS Photonics , 992–1005 (2018). C. Sch¨afer, M. Ruggenthaler, H. Appel, and A. Rubio, “Modification of excitation andcharge transfer in cavity quantum-electrodynamical chemistry,” Proceedings of the Na-tional Academy of Sciences , 4883–4892 (2019). A. Castro, A first-principles time-dependent density functional theory scheme for the computation of the electromagnetic response of nanostructures,Ph.D. thesis, PhD. Thesis, University of Valladolid, 2004. The thesis can be down-loaded (2004). X. Andrade, “Linear and non-linear response phenomena of molecular systems withintime-dependent density functional theory,” (2010). M. A. Marques, A. Castro, G. F. Bertsch, and A. Rubio, “octopus: a first-principlestool for excited electronion dynamics,” Computer Physics Communications , 60 – 78(2003). A. Castro, H. Appel, M. Oliveira, C. A. Rozzi, X. Andrade, F. Lorenzen, M. A. L.Marques, E. K. U. Gross, and A. Rubio, “octopus: a tool for the application of time-dependent density functional theory,” physica status solidi (b) , 2465–2488 (2006). X. Andrade, D. Strubbe, U. De Giovannini, A. H. Larsen, M. J. T. Oliveira, J. Alberdi-Rodriguez, A. Varas, I. Theophilou, N. Helbig, M. J. Verstraete, L. Stella, F. Nogueira,A. Aspuru-Guzik, A. Castro, M. A. L. Marques, and A. Rubio, “Real-space grids and theOctopus code as tools for the development of new simulation approaches for electronicsystems,” Phys. Chem. Chem. Phys. , 31371–31396 (2015). The Octopus code, tutorials, and examples can be found on its website at https:// ctopus-code.org/ . R. Jest¨adt, M. Ruggenthaler, M. J. T. Oliveira, A. Rubio, and H. Appel, “Light-matter interactions within the Ehrenfest-Maxwell-Pauli-Kohn-Sham framework: funda-mentals, implementation, and nano-optical applications,” Advances in Physics (in press),arXiv:1812.05049. G. F. Bertsch, J.-I. Iwata, A. Rubio, and K. Yabana, “Real-space, real-time method forthe dielectric function,” Phys. Rev. B , 7998–8002 (2000). A. Taflove and S. Hagness, Computational Electrodynamics: The Finite-difference Time-domain Method,Artech House antennas and propagation library (Artech House, 2005). I. Bialynicki-Birula, “On the wave function of the photon,” Acta Physica Polonica-SeriesA General Physics , 97–116 (1994). S. Ludwig, “Elektromagnetische grundgleichungen in bivektorieller behandlung,” Annalender Physik , 579–586 (1907). R. Loudon, The quantum theory of light (Oxford Science Publications, 1988). D. P. Craig and T. Thirunamachandran, Molecular quantum electrodynamics: an introduction to radiation-molecule interactions(Courier Corporation, 1998). J. J. Baumberg, J. Aizpurua, M. H. Mikkelsen, and D. R. Smith, “Extreme nanophotonicsfrom ultrathin metallic gaps,” Nature Materials , 668–678 (2019). C. Pellegrini, J. Flick, I. V. Tokatly, H. Appel, and A. Rubio, “Optimized effectivepotential for quantum electrodynamical time-dependent density functional theory,” Phys.Rev. Lett. , 093001 (2015). S. K¨ummel and L. Kronik, “Orbital-dependent density functionals: Theory and applica-tions,” Rev. Mod. Phys. , 3–60 (2008). S. K¨ummel and J. P. Perdew, “Simple iterative construction of the optimized effectivepotential for orbital functionals, including exact exchange,” Phys. Rev. Lett. , 043004(2003). T. W. Hollins, S. J. Clark, K. Refson, and N. I. Gidopoulos, “Optimized effective potentialusing the hylleraas variational method,” Phys. Rev. B , 235126 (2012). J. Krieger, Y. Li, and G. Iafrate, “Derivation and application of an accurate kohn-shampotential with integer discontinuity,” Physics Letters A , 256 – 260 (1990). J. B. Krieger, Y. Li, and G. J. Iafrate, “Construction and application of an accurate localspin-polarized kohn-sham potential with integer discontinuity: Exchange-only theory,”69hys. Rev. A , 101–126 (1992). J. B. Krieger, Y. Li, and G. J. Iafrate, “Systematic approximations to the optimizedeffective potential: Application to orbital-density-functional theory,” Phys. Rev. A ,5453–5458 (1992). C. Sch¨afer, M. Ruggenthaler, and A. Rubio, “Ab initio nonrelativistic quantum electrody-namics: Bridging quantum chemistry and quantum optics from weak to strong coupling,”Phys. Rev. A , 043801 (2018). V. Rokaj, D. M. Welakuh, M. Ruggenthaler, and A. Rubio, “Lightmatter interaction inthe long-wavelength limit: no ground-state without dipole self-energy,” Journal of PhysicsB: Atomic, Molecular and Optical Physics , 034005 (2018). J. Flick and P. Narang, “Cavity-correlated electron-nuclear dynamics from first princi-ples,” Phys. Rev. Lett. , 113002 (2018). J. Flick, D. M. Welakuh, M. Ruggenthaler, H. Appel, and A. Rubio, “Light-matterresponse functions in quantum-electrodynamical density-functional theory: Modificationsof spectra and of the maxwell equations,” accepted for publication in ACS Photonics(2019), arXiv:1803.02519. N. M. Hoffmann, C. Sch¨afer, A. Rubio, A. Kelly, and H. Appel, “Capturing vacuum fluc-tuations and photon correlations in cavity quantum electrodynamics with multitrajectoryehrenfest dynamics,” Phys. Rev. A , 063819 (2019). J. Galego, C. Climent, F. J. Garcia-Vidal, and J. Feist, “Cavity Casimir-Polder forcesand their effects in ground state chemical reactivity,” Phys. Rev. X , 1–22 (2019). C. Sch¨afer, M. Ruggenthaler, and A. Rubio, “Ab initio nonrelativistic quantum electrody-namics: Bridging quantum chemistry and quantum optics from weak to strong coupling,”Physical Review A , 043801 (2018). Note that Hartree-Fock theory is also included within RDMFT by fixing the orbitaloccupations to 0 and 1. F. Buchholz, I. Theophilou, S. E. B. Nielsen, M. Ruggenthaler, and A. Rubio, “Reduceddensity-matrix approach to strong matter-photon interaction,” ACS Photonics , 2694–2711 (2019). X. Andrade, D. A. Strubbe, U. De Giovannini, A. H. Larsen, M. J. T. Oliveira, J. Alberdi-Rodriguez, A. Varas, I. Theophilou, N. Helbig, M. Verstraete, L. Stella, F. Nogueira,A. Aspuru-Guzik, A. Castro, M. A. L. Marques, and ´A. Rubio, “Real-space grids and70he Octopus code as tools for the development of new simulation approaches for electronicsystems,” Phys. Chem. Chem. Phys. , 31371–31396 (2015). At the current state of the code, we have only implemented the so-called M¨uller func-tional, but one could potentially use any RDMFT functional to approximate two-bodyexpression. A. J. Coleman, “Structure of fermion density matrices,” Rev. Mod. Phys. , 668–686(1963). The respective results for the equilibrium distance b = 1 .
61 bohr are shown in Ref. 44,Sec. 7. Performed with the many-body routine of Octopus, see Ref. 19, Sec. 13 for details. V. I. Anisimov, J. Zaanen, and O. K. Andersen, “Band theory and Mott insulators:Hubbard U instead of Stoner I,” Physical Review B , 943 (1991). V. I. Anisimov, I. V. Solovyev, M. A. Korotin, M. T. Czy˙zyk, and G. A. Sawatzky,“Density-functional theory and NiO photoemission spectra,” Physical Review B , 16929(1993). A. I. Liechtenstein, V. I. Anisimov, and J. Zaanen, “Density-functional theory and stronginteractions: Orbital ordering in Mott-Hubbard insulators,” Physical Review B , R5467(1995). V. I. Anisimov, F. Aryasetiawan, and A. Lichtenstein, “First-principles calculations ofthe electronic structure and spectra of strongly correlated systems: the lda+ u method,”Journal of Physics: Condensed Matter , 767 (1997). K. Haule, “Exact double counting in combining the dynamical mean field theory and thedensity functional theory,” Phys. Rev. Lett. , 196403 (2015). A. G. Petukhov, I. I. Mazin, L. Chioncel, and A. I. Lichtenstein, “Correlated metals andthe LDA + u method,” Phys. Rev. B , 153106 (2003). S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys, and A. P. Sutton,“Electron-energy-loss spectra and the structural stability of nickel oxide: An lsda+ustudy,” Phys. Rev. B , 1505–1509 (1998). M. T. Czy˙zyk and G. A. Sawatzky, “Local-density functional and on-site correlations:The electronic structure of la cuo and lacuo ,” Phys. Rev. B , 14211–14228 (1994). N. Tancogne-Dejean, M. J. T. Oliveira, and A. Rubio, “Self-consistent DFT + u methodfor real-space time-dependent density functional theory calculations,” Phys. Rev. B ,7145133 (2017). L. Xian, D. M. Kennes, N. Tancogne-Dejean, M. Altarelli, and A. Rubio, “Multiflatbands and strong correlations in twisted bilayer boron nitride: Doping-induced correlatedinsulator and superconductor,” Nano letters , 4934–4940 (2019). N. Tancogne-Dejean and A. Rubio, “Parameter-free hybrid functional based on an ex-tended hubbard model: Dft+u+v,” arXiv:1911.10813 (2019). L. A. Agapito, S. Curtarolo, and M. Buongiorno Nardelli, “Reformulation of DFT + u as a pseudohybrid hubbard density functional for accelerated materials discovery,” Phys.Rev. X , 011006 (2015). G. E. Topp, N. Tancogne-Dejean, A. F. Kemper, A. Rubio, and M. A. Sentef, “All-opticalnonequilibrium pathway to stabilising magnetic weyl semimetals in pyrochlore iridates,”Nature communications , 4452 (2018). N. Tancogne-Dejean, M. A. Sentef, and A. Rubio, “Ultrafast modification of hubbard u in a strongly correlated material: Ab initio high-harmonic generation in nio,” Phys. Rev.Lett. , 097402 (2018). K. Berland, V. R. Cooper, K. Lee, E. Schr¨oder, T. Thonhauser, P. Hyldgaard, and B. I.Lundqvist, “van der waals forces in density functional theory: a review of the vdw-dfmethod,” Reports on Progress in Physics , 066501 (2015). K. Berland, V. R. Cooper, K. Lee, E. Schr¨oder, T. Thonhauser, P. Hyldgaard, and B. I.Lundqvist, “van der Waals forces in density functional theory: a review of the vdW-DFmethod,” Reports on Progress in Physics , 066501 (2015). S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, “A consistent and accurate ab initioparametrization of density functional dispersion correction (dft-d) for the 94 elementsh-pu,” The Journal of Chemical Physics , 154104 (2010). A. Tkatchenko and M. Scheffler, “Accurate molecular van der waals interactions fromground-state electron density and free-atom reference data,” Phys. Rev. Lett. , 073005(2009). A. H. Larsen, M. Kuisma, J. Lfgren, Y. Pouillon, P. Erhart, and P. Hyldgaard, “libvdwxc:a library for exchange–correlation functionals in the vdW-DF family,” Modelling andSimulation in Materials Science and Engineering , 065004 (2017). J. P. Perdew and Y. Wang, “Accurate and simple analytic representation of the electron-gas correlation energy,” Physical Review B , 13244–13249 (1992).72 G. Rom´an-P´erez and J. M. Soler, “Efficient implementation of a van der waals densityfunctional: Application to double-wall carbon nanotubes,” Physical Review Letters ,096102 (2009). M. Frigo and S. G. Johnson, “The design and implementation of FFTW3,” Proceedingsof the IEEE , 216–231 (2005), special issue on “Program Generation, Optimization,and Platform Adaptation”. M. Dion, H. Rydberg, E. Schr¨oder, D. C. Langreth, and B. I. Lundqvist, “Van der Waalsdensity functional for general geometries,” Physical Review Letters , 246401 (2004). M. Dion, H. Rydberg, E. Schr¨oder, D. C. Langreth, and B. I. Lundqvist, “Erratum: Vander Waals density functional for general geometries [phys. rev. lett. , 246401 (2004)],”Physical Review Letters , 109902 (2005). K. Lee, E. D. Murray, L. Kong, B. I. Lundqvist, and D. C. Langreth, “Higher-accuracyvan der Waals density functional,” Physical Review B , 081101 (2010). K. Berland and P. Hyldgaard, “Exchange functional that tests the robustness of theplasmon description of the van der Waals density functional,” Physical Review B ,035412 (2014). J. Klimeˇs, D. R. Bowler, and A. Michaelides, “Chemical accuracy for the van der Waalsdensity functional,” Journal of Physics: Condensed Matter , 022201 (2010). V. R. Cooper, “van der Waals density functional: An appropriate exchange functional,”Physical Review B , 161104 (2010). J. Tomasi, B. Mennucci, and R. Cammi, “Quantum mechanical continuum solvationmodels,” Chem. Rev. , 2999–3094 (2005). J. Tomasi and M. Persico, “Molecular interactions in solution: An overview of methodsbased on continuous distributions of the solvent,” Chem. Rev. , 2027–2094 (1994). J. Tomasi, “The physical model,” in Continuum Solvation Models in Chemical Physics: From Theory to Applications,edited by B. Mennucci and R. Cammi (Wiley & Sons, Chichester, UK, 2007) Chap. Sect.1.1, p. 16. A. Delgado, S. Corni, S. Pittalis, and C. A. Rozzi, “Modeling solvation effects in real-spaceand real-time within density functional approaches,” The Journal of Chemical Physics , 144111 (2015). E. Canc´es, B. Mennucci, and J. Tomasi, “A new integral equation formalism for thepolarizable continuum model: Theoretical background and applications to isotropic and73nisotropic dielectrics,” J. Chem. Phys. , 3032 (1997). J.-L. Pascual-ahuir, E. Silla, and I. Tunon, “Gepol: An improved description of molecularsurfaces. iii. a new algorithm for the computation of a solvent-excluding surface,” J.Comput. Chem. , 1127–1138 (1994). M. W. Schmidt, K. K. Baldridge, J. A. Boatz, S. T. Elbert, M. S. Gordon, J. H. Jensen,S. Koseki, N. Matsunaga, K. A. Nguyen, S. Su, et al., “General atomic and molecularelectronic structure system,” Journal of Computational Chemistry , 1347–1363 (1993). J. Perdew, K. Burke, and M. Ernzerhof, “Perdew, burke, and ernzerhof reply,” PhysicalReview Letters , 891 (1998). M. Caricato, F. Ingrosso, B. Mennucci, and J. Tomasi, “A time-dependent polarizablecontinuum model: theory and application.” J Chem Phys , 154501 (2005). R. Cammi and J. Tomasi, “Nonequilibrium solvation theory for the polarizable continuummodel: A new formulation at the SCF level with application to the case of the frequency-dependent linear electric response function,” Int. J. Quantum Chem. , 465–474 (1995). S. Corni, S. Pipolo, and R. Cammi, “Equation of motion for the solvent polarizationapparent charges in the polarizable continuum model: Application to real-time tddft.” J.Phys. Chem. A , 5405–5416 (2014). G. Gil, S. Pipolo, A. Delgado, C. A. Rozzi, and S. Corni, “Nonequilibrium solventpolarization effects in real-time electronic dynamics of solute molecules subject to time-dependent electric fields: A new feature of the polarizable continuum model,” Journal ofChemical Theory and Computation , 2306–2319 (2019). L. Onsager, “Electric moments of molecules in liquids,” J. Am. Chem. Soc. , 1486–1493(1936). S. Corni, R. Cammi, B. Mennucci, and J. Tomasi, “Electronic excitation energies ofmolecules in solution within continuum solvation models: Investigating the discrepancybetween state-specific and linear-response methods,” J. Chem. Phys. (2005). P. Buczek, A. Ernst, and L. M. Sandratskii, “Different dimensionality trends in the lan-dau damping of magnons in iron, cobalt, and nickel: Time-dependent density functionalstudy,” Phys. Rev. B , 174418 (2011). M. Niesert, Ab initio Calculations of Spin-Wave Excitation Spectra from Time-Dependent Density-Functional Theory,Dr. (univ.), RWTH Aachen (2011), record converted from VDB: 12.11.2012; Aachen,RWTH, Diss., 2011. 74 B. Rousseau, A. Eiguren, and A. Bergara, “Efficient computation of magnon disper-sions within time-dependent density functional theory using maximally localized wannierfunctions,” Phys. Rev. B , 054305 (2012). T. Gorni, I. Timrov, and S. Baroni, “Spin dynamics from time-dependent density func-tional perturbation theory,” The European Physical Journal B , 249 (2018). K. Cao, H. Lambert, P. G. Radaelli, and F. Giustino, “Ab initio calculation of spinfluctuation spectra using time-dependent density functional perturbation theory, planewaves, and pseudopotentials,” Phys. Rev. B , 024420 (2018). N. Singh, P. Elliott, T. Nautiyal, J. K. Dewhurst, and S. Sharma, “Adiabatic generalizedgradient approximation kernel in time-dependent density functional theory,” Phys. Rev.B , 035151 (2019). S. Y. Savrasov, “Linear response calculations of spin fluctuations,” Phys. Rev. Lett. ,2570–2573 (1998). N. Tancogne-Dejean, F. G. Eich, and A. Rubio, “Time-dependent magnons from firstprinciples,” Submitted (2019).
L. M. Sandratskii, “Energy band structure calculations for crystals with spiral magneticstructure,” physica status solidi (b) , 167–180 (1986).
H. A. Mook and D. M. Paul, “Neutron-scattering measurement of the spin-wave spectrafor nickel,” Phys. Rev. Lett. , 227–229 (1985). U. De Giovannini, G. Brunetto, A. Castro, J. Walkenhorst, and A. Rubio, “Simulatingpump–probe photoelectron and absorption spectroscopy on the attosecond timescale withtime-dependent density functional theory,” ChemPhysChem , 1363–1376 (2013). J. Walkenhorst, U. De Giovannini, A. Castro, and A. Rubio, “Tailored pump-probe tran-sient spectroscopy with time-dependent density-functional theory: controlling absorptionspectra,” The European Physical Journal B , 128 (2016). U. De Giovannini, “Pump-probe photoelectron spectra,” Handbook of Materials Model-ing: Methods: Theory and Modeling , 1–19 (2018).
S. Sato, H. H¨ubener, U. De Giovannini, and A. Rubio, “Ab initio simulation of attosecondtransient absorption spectroscopy in two-dimensional materials,” Applied Sciences , 1777(2018). L. D. Barron, Molecular Light Scattering and Optical Activity, 2nd ed. (Cambridge Uni-versity Press, Cambridge, 2004). 75 F. R. Keßler and J. Metzdorf, “Landau level spectroscopy: Interband effects and Fara-day rotation,” in Modern Problems in Condensed Matter Sciences, Vol. 27.1, edited byV. M. Agranovich and A. A. Maradudin (Elsevier Science Publishers, Amsterdam, 1991)Chap. 11.
H. Solheim, K. Ruud, S. Coriani, and P. Norman, “Complex polarization propagatorcalculations of magnetic circular dichroism spectra,” J. Chem. Phys. , 094103 (2008).
H. Solheim, K. Ruud, S. Coriani, and P. Norman, “The A and B terms of magneticcircular dichroism revisited,” J. Phys. Chem. A , 9615–9618 (2008).
K.-M. Lee, K. Yabana, and G. F. Bertsch, “Magnetic circular dichroism in real-timetime-dependent density functional theory,” J. Chem. Phys. , 144106 (2011).
M. Seth, M. Krykunov, T. Ziegler, J. Autschbach, and A. Banerjee, “Application ofmagnetically perturbed time-dependent density functional theory to magnetic circulardichroism: Calculation of B terms,” J. Chem. Phys. , 144105 (2008).
M. Seth, M. Krykunov, T. Ziegler, and J. Autschbach, “Application of magneticallyperturbed time-dependent density functional theory to magnetic circular dichroism. II.Calculation of A terms,” J. Chem. Phys. , 234102 (2008).
I. V. Lebedeva, D. A. Strubbe, I. V. Tokatly, and A. Rubio, “Orbital magneto-opticalresponse of periodic insulators from first principles,” npj Computational Materials , 32(2019). R. D. King-Smith and D. Vanderbilt, “Theory of polarization of crystalline solids,” Phys.Rev. B , 1651–1654 (1993). D. Vanderbilt and R. D. King-Smith, “Electric polarization as a bulk quantity and itsrelation to surface charge,” Phys. Rev. B , 4442–4455 (1993). R. Resta, “Macroscopic polarization in crystalline dielectrics: the geometric phase ap-proach,” Rev. Mod. Phys. , 899–915 (1994). A. M. Essin, A. M. Turner, J. E. Moore, and D. Vanderbilt, “Orbital magnetoelectriccoupling in band insulators,” Phys. Rev. B , 205104 (2010). K.-T. Chen and P. A. Lee, “Unified formalism for calculating polarization, magnetization,and more in a periodic insulator,” Phys. Rev. B , 205137 (2011). X. Gonze and J. W. Zwanziger, “Density-operator theory of orbital magnetic susceptibilityin periodic insulators,” Phys. Rev. B , 064445 (2011). X. Andrade, S. Botti, M. A. L. Marques, and A. Rubio, “Time-dependent density76unctional theory scheme for efficient calculations of dynamic (hyper)polarizabilities,”J. Chem. Phys. , 184106 (2007).
D. A. Strubbe, Optical and transport properties of organic molecules: Methods and applications(PhD thesis, University of California, Berkeley, USA, 2012).
D. A. Strubbe, L. Lehtovaara, A. Rubio, M. A. L. Marques, and S. G.Louie, “Response functions in TDDFT: Concepts and implementation,” inFundamentals of Time-Dependent Density Functional Theory (Springer Berlin Hei-delberg, Berlin, Heidelberg, 2012) pp. 139–166.
M. Lazzeri and F. Mauri, “High-order density-matrix perturbation theory,” Phys. Rev. B , 161101 (2003). X. Gonze and J.-P. Vigneron, “Density-functional approach to nonlinear-response coeffi-cients of solids,” Phys. Rev. B , 13120–13128 (1989). A. D. Corso, F. Mauri, and A. Rubio, “Density-functional theory of the nonlinear opti-cal susceptibility: Application to cubic semiconductors,” Phys. Rev. B , 15638–15642(1996). J. A. Berger, “Fully parameter-free calculation of optical spectra for insulators, semicon-ductors, and metals from a simple polarization functional,” Phys. Rev. Lett. , 137402(2015).
S. Albrecht, L. Reining, R. Del Sole, and G. Onida, “Ab initio calculation of excitoniceffects in the optical spectra of semiconductors,” Phys. Rev. Lett. , 4510–4513 (1998). S. Botti, F. Sottile, N. Vast, V. Olevano, L. Reining, H.-C. Weissker, A. Rubio, G. Onida,R. Del Sole, and R. W. Godby, “Long-range contribution to the exchange-correlationkernel of time-dependent density functional theory,” Phys. Rev. B , 155112 (2004). P. Lautenschlager, M. Garriga, L. Vi˜na, and M. Cardona, “Temperature dependence ofthe dielectric function and interband critical points in silicon,” Phys. Rev. B , 4821–4830 (1987). U. De Giovannini and A. Castro, “Real-time and real-space time-dependent density-functional theory approach to attosecond dynamics,” in Attosecond Molecular Dynamics(Royal Society of Chemistry, Cambridge, 2018) pp. 424–461.
U. De Giovannini, H. H¨ubener, and A. Rubio, “A First-Principles Time-Dependent Den-sity Functional Theory Framework for Spin and Time-Resolved Angular-Resolved Photo-electron Spectroscopy in Periodic Systems,” Journal of Chemical Theory and Computa-77ion , 265–273 (2017). L. Tao and A. Scrinzi, “Photo-electron momentum spectra from minimal volumes: thetime-dependent surface flux method,” New Journal of Physics , 013021 (2012). P. Wopperer, U. De Giovannini, and A. Rubio, “Efficient and accurate modeling ofelectron photoemission in nanostructures with TDDFT,” The European Physical JournalB , 1307 (2017). U. De Giovannini, A. H. Larsen, A. Rubio, and A. Rubio, “Modeling electron dynamicscoupled to continuum states in finite volumes with absorbing boundaries,” The EuropeanPhysical Journal B , 1–12 (2015). S. A. Sato, H. H¨ubener, A. Rubio, and U. De Giovannini, “First-principles simulationsfor attosecond photoelectron spectroscopy based on time-dependent density functionaltheory,” The European Physical Journal B , 126 (2018). X. Andrade, A. Castro, D. Zueco, J. L. Alonso, P. Echenique, F. Falceto, and A. Rubio,“Modified Ehrenfest Formalism for Efficient Large-Scale ab initio Molecular Dynamics,”Journal of Chemical Theory and Computation , 728–742 (2009). H. H¨ubener, U. De Giovannini, and A. Rubio, “Phonon Driven Floquet Matter,” NanoLett , 1535–1542 (2018). U. De Giovannini, H. H¨ubener, and A. Rubio, “Monitoring electron-photon dressing inWSe2,” Nano Letters , 7993–7998 (2016). S. K. Takashi Oka, “Floquet Engineering of Quantum Materials,” Annual Review ofCondensed Matter Physics , annurev–conmatphys–031218–013423 (2018). H. H. Umberto De Giovannini, “Floquet analysis of excitations in materials,” Journal ofPhysics: Materials (2019), 10.1088/2515-7639/ab387b.
C. I. Blaga, J. Xu, A. D. DiChiara, E. Sistrunk, K. Zhang, P. Agostini, T. A. Miller, L. F.DiMauro, and C. D. Lin, “Imaging ultrafast molecular dynamics with laser-inducedelectron diffraction,” Nature , 194–197 (2012).
F. Kreˇcini´c, P. Wopperer, B. Frusteri, F. Brauße, J.-G. Brisset, U. De Giovannini, A. Ru-bio, A. Rouz´ee, and M. J. J. Vrakking, “Multiple-orbital effects in laser-induced electrondiffraction of aligned molecules,” Physical Review A , 041401 (2018). A. Trabattoni, S. Trippel, U. De Giovannini, J. F. Olivieri, J. Wiese, T. Mullins, J. Onvlee,S.-K. Son, B. Frusteri, A. Rubio, and J. K¨upper, “Setting the clock of photoelectronemission through molecular alignment,” e-print (2018), arXiv:1802.06622.78 V. Recoules, F. Lambert, A. Decoster, B. Canaud, and J. Cl´erouin, “Ab initio determi-nation of thermal conductivity of dense hydrogen plasmas,” Physical review letters ,075002 (2009).
F. Lambert, V. Recoules, A. Decoster, J. Clerouin, and M. Desjarlais, “On the transportcoefficients of hydrogen in the inertial confinement fusion regime,” Physics of Plasmas ,056306 (2011). X. Andrade, S. Hamel, and A. A. Correa, “Negative differential conductivity in liquidaluminum from real-time quantum simulations,” The European Physical Journal B ,229 (2018). F. Eich, M. Di Ventra, and G. Vignale, “Functional theories of thermoelectric phenom-ena,” Journal of Physics: Condensed Matter , 063001 (2016). P. Hohenberg and W. Kohn, “Inhomogeneous electron gas,” Phys. Rev. , B864–B871(1964).
E. Runge and E. K. U. Gross, “Density-functional theory for time-dependent systems,”Phys. Rev. Lett. , 997–1000 (1984). J. Neugebauer, “Chromophore-specific theoretical spectroscopy: From subsystem densityfunctional theory to mode-specific vibrational spectroscopy,” Physics Reports , 1 – 87(2010).
R. F. W. Bader, Atoms in Molecules. A Quantum Theory (Clarendon Press, 1994).
A. Krishtal, D. Ceresoli, and M. Pavanello, “Subsystem real-time time dependent densityfunctional theory,” The Journal of Chemical Physics , 154116 (2015).
J. P. Perdew, K. Burke, and M. Ernzerhof, “Generalized gradient approximation madesimple,” Phys. Rev. Lett. , 3865–3868 (1996). J. P. Perdew, K. Burke, and M. Ernzerhof, “Generalized Gradient Approximation MadeSimple [Phys. Rev. Lett. 77, 3865 (1996)],” Phys. Rev. Lett. , 1396 (1997). M. Schlipf and F. Gygi, “Optimization algorithm for the generation of oncv pseudopo-tentials,” Computer Physics Communications , 36 – 44 (2015).
E. F. Pettersen, T. D. Goddard, C. C. Huang, G. S. Couch, D. M. Greenblatt, E. C.Meng, and T. E. Ferrin, “Ucsf chimeraa visualization system for exploratory researchand analysis,” Journal of Computational Chemistry , 1605–1612 (2004). J. Jornet-Somoza, J. Alberdi-Rodriguez, B. F. Milne, X. Andrade, M. A. L. Marques,F. Nogueira, M. J. T. Oliveira, J. J. P. Stewart, and A. Rubio, “Insights into colour-79uning of chlorophyll optical response in green plants,” Phys. Chem. Chem. Phys. ,26599–26606 (2015). J. Jornet-Somoza and I. Lebedeva, “Real-time propagation tddft and density analysis forexciton coupling calculations in large systems,” Journal of Chemical Theory and Compu-tation , 3743–3754 (2019). A. G´omez Pueyo, M. A. L. Marques, A. Rubio, and A. Castro, “Propagators for thetime-dependent kohnsham equations: Multistep, rungekutta, exponential rungekutta, andcommutator free magnus methods,” Journal of Chemical Theory and Computation ,3040–3052 (2018). E. Hairer, C. Lubich, and G. Wanner, Geometric Numerical Integration (Springer Verlag,Berlin Heidelberg, 2006).
A. Castro, M. A. L. Marques, and A. Rubio, “Propagators for the time-dependent kohn–sham equations,” The Journal of Chemical Physics , 3425–3433 (2004).
S. Blanes and P. Moan, “Fourth- and sixth-order commutator-free magnus integrators forlinear and non-linear dynamical systems,” Applied Numerical Mathematics , 1519 –1537 (2006). M. Piris and J. M. Ugalde, “Iterative diagonalization for orbital optimization in naturalorbital functional theory,” Journal of Computational Chemistry , 2078–2086 (2009). A. M. K. M¨uller, “Explicit approximate relation between reduced two- and one-particledensity matrices,” Phys. Rev. A , 446–452 (1984).
N. Helbig, J. I. Fuks, M. Casula, M. J. Verstraete, M. A. L. Marques, I. V. Tokatly, andA. Rubio, “Density functional theory beyond the linear regime: Validating an adiabaticlocal density approximation,” Phys. Rev. A , 032503 (2011). M. M. Morrell, R. G. Parr, and M. Levy, “Calculation of ionization potentials fromdensity matrices and natural functions, and the longrange behavior of natural orbitalsand electron density,” The Journal of Chemical Physics , 549–554 (1975). M. C. Payne, M. P. Teter, D. C. Allan, T. A. Arias, and J. D. Joannopoulos, “Iterativeminimization techniques for ab initio total-energy calculations: molecular dynamics andconjugate gradients,” Rev. Mod. Phys. , 1045–1097 (1992). R. L. Frank, E. H. Lieb, R. Seiringer, and H. Siedentop, “M¨uller’s exchange-correlationenergy in density-matrix-functional theory,” Phys. Rev. A , 052517 (2007). J. R. Chelikowsky, N. Troullier, and Y. Saad, “Finite-difference-pseudopotential method:80lectronic structure calculations without a basis,” Phys. Rev. Lett. , 1240–1243 (1994). A. Natan, A. Benjamini, D. Naveh, L. Kronik, M. L. Tiago, S. P. Beckman, and J. R.Chelikowsky, “Real-space pseudopotential method for first principles calculations of gen-eral periodic and partially periodic systems,” Phys. Rev. B , 075109 (2008). A. Togo and I. Tanaka, “
Spglib : a software library for crystal symmetry search,” e-print(2018), arXiv:1808.01590.
K. Lejaeghere, G. Bihlmayer, T. Bj¨orkman, P. Blaha, S. Bl¨ugel, V. Blum, D. Caliste, I. E.Castelli, S. J. Clark, A. Dal Corso, et al., “Reproducibility in density functional theorycalculations of solids,” Science , aad3000 (2016).
J. L. Alonso, X. Andrade, P. Echenique, F. Falceto, D. Prada-Gracia, and A. Rubio, “Effi-cient Formalism for Large-Scale
Ab Initio
Molecular Dynamics based on Time-DependentDensity Functional Theory,” Physical Review Letters , 096403 (2008).
M. Parrinello and A. Rahman, “Crystal structure and pair potentials: A molecular-dynamics study,” Phys. Rev. Lett. , 1196–1199 (1980). M. Parrinello and A. Rahman, “Polymorphic transitions in single crystals: A new molec-ular dynamics method,” Journal of Applied Physics , 7182–7190 (1981). R. M. Martin, Electronic structure: basic theory and practical methods (Cambridge uni-versity press, 2004). http://buildbot.net . Y. Saad, A. Stathopoulos, J. Chelikowsky, K. Wu, and S. ¨O˘g¨ut, “Solution of largeeigenvalue problems in electronic structure calculations,” BIT Numerical Mathematics , 563–578 (1996). A. Meister and C. V¨omel, Numerik linearer Gleichungssysteme (Vieweg+Teubner Verlag,2011).
J. J. Mortensen, L. B. Hansen, and K. W. Jacobsen, “Real-space grid implementation ofthe projector augmented wave method,” Phys. Rev. B , 035109 (2005). A. P. Seitsonen, M. J. Puska, and R. M. Nieminen, “Real-space electronic-structurecalculations: Combination of the finite-difference and conjugate-gradient methods,” Phys.Rev. B , 14057–14061 (1995). X. Andrade and A. Aspuru-Guzik, “Real-space density functional theory on graphical pro-cessing units: Computational approach and comparison to gaussian basis set methods,”Journal of Chemical Theory and Computation9