Octopus, a computational framework for exploring light-driven phenomena and quantum dynamics in extended and finite systems

I. INTRODUCTION

• Octopus’s recent additions span microscopic light–matter coupling, solvent and van der Waals treatments, real-time spectroscopy and response methods, and numerical and periodic-system improvements.The additions include methods for magnons, conductivities, photoelectron spectra, and orbital magneto-optical responses, as well as new propagators and iterative eigensolvers for reduced density matrix functional theory.Source: §I · p. 5
• The paper covers new physical theories and algorithms for nonequilibrium phenomena, followed by technical developments aimed at code performance and algorithm stability.Source: §I · p. 5

II. COUPLED MAXWELL-KOHN-SHAM EQUATIONS

• Octopus’s prior forward-coupling approach prescribed a classical laser pulse and evolved Kohn–Sham orbitals while neglecting the matter’s back-action on the pulse.The approximation is highly accurate when the induced total current is comparatively small, such as in atoms or smaller molecules.Source: §II · pp. 5–6
• Matter-only and field-only approaches break down when induced currents are large and material-scale electromagnetic near-field effects are non-negligible.Source: §II · p. 6
• Starting from a generalized Pauli–Fierz theory of electrons, nuclei, and photons, the authors derive self-consistent coupling levels and implement a Maxwell solver coupled to electron and nuclear dynamics.In the classical limit, the framework yields coupled Ehrenfest–Maxwell–Pauli–Kohn–Sham equations.Source: §II · p. 6
• Unlike Yee-based FDTD, Octopus uses one grid for the Riemann-Silberstein vector, simplifying matter coupling and enabling higher-order finite differences, larger time steps, and improved stability.The Maxwell Hamiltonian acts on six orbitals corresponding to the electric- and magnetic-field components.Source: §II · pp. 7–8
• Maxwell’s equations couple directly to microscopic matter currents, while the Power-Zienau-Woolley multipole expansion includes dipole, electrical-quadrupole, and magnetic-dipole terms.The Kohn–Sham orbitals and Maxwell fields evolve side by side and are coupled self-consistently at each time step.Source: §II · p. 8
• Because mask absorption proved insufficient for outgoing electromagnetic waves, Octopus implements a perfectly matched layer to reduce boundary reflections and backscattering.Source: §II · p. 8
• For optical-wavelength fields, Octopus supports multiscale space-time propagation, including a larger Maxwell box with interpolation and frozen-current electromagnetic steps between matter updates.The number of intermediate steps is a convergence parameter that depends on the physical situation.Source: §II · p. 9
• Explicitly propagating classical electromagnetic fields gives access to outgoing radiation, which Octopus detectors accumulate at box boundaries for direct observation and spectra on the Maxwell grid.This avoids Fourier-transforming the matter-dipole time signal to obtain optical spectra.Source: §II · p. 10
• Back-reaction reduces the field enhancement between the sodium nanoparticles but increases it farther away, with frequency shifts more pronounced in the near field.The reported field enhancement is also sensitive to the included multipole-coupling terms; the quantitative difference between LDA and PBE is smaller than that from including back-reaction.Source: §II · p. 10, Figure 2 · p. 11
• The coupled real-time, real-space Ehrenfest-Maxwell-Pauli-Kohn-Sham implementation enables fully self-consistent forward-backward light–matter coupling and leads to the classical limit of QEDFT.The authors describe it as a tool for applications in close analogy to experiments.Source: §II · p. 10
• The system is driven by a z field of 5.142 × 10^7 V/m (10^-4 a.u.), with laser period T1 = 1.36 fs.Source: §II · p. 11
• The panels show current density at rcp and field enhancement at rcp, rocpx, and averaged over a detector surface near the box boundary.In panel c, a bright-gray curve identical to the laser pulse is included for comparison.Source: §II · p. 11

III. STRONG ELECTRON-PHOTON INTERACTIONS IN REAL SPACE: QUANTUM-ELECTRODYNAMICAL DENSITY-FUNCTIONAL THEORY (QEDFT)

• QEDFT extends time-dependent density-functional theory to ab initio weak and strong light–matter interactions, enabling exploration of photon-dressed electronic states in real space.The framework is presented for applications to chemistry and materials and for exploring, predicting, and controlling nonequilibrium states of matter [4,10].Source: §III · p. 12
• The light–matter Hamiltonian treats N interacting electrons coupled to Np photonic modes in the length gauge under the long-wavelength dipole approximation.Each mode is characterized by its elongation, frequency, and coupling-strength vector, which includes its polarization.Source: §III · p. 12
• QEDFT expands the internal variables from electronic density to include each photonic mode’s electric-displacement contribution, mapped bijectively to external electronic and photonic variables.The resulting auxiliary system combines electronic Kohn–Sham equations with Maxwell equations.Source: §III · pp. 12–13
• The Kohn–Sham reformulation places light–matter quantumness in the local exchange-correlation potential, which gains an electron–photon component while the photon exchange-correlation contribution vanishes.The passage describes calculating this local potential as the key additional requirement for including quantum photons.Source: §III · p. 13
• Initial electron–photon exchange-correlation potentials were constructed using many-body perturbation theory with exact exchange and the optimized-effective-potential method [13,29,30].Source: §III · p. 13

A. The optimized effective potential (OEP)

• The final OEP equation includes electron-electron and electron-photon effects and defines an exchange-correlation orbital shift for constructing the corresponding potential.The exchange-correlation shift is defined along the lines of the orbital shift used in conventional OEP calculations.Source: A. The optimized effective potential… · p. 14
• The reformulation replaces the all-unoccupied-state OEP calculation with Np+1 self-consistent Sternheimer equations involving occupied orbitals only.The equations are solved self-consistently with Eq. (11).Source: A. The optimized effective potential… · p. 15
• The convergence measure Sσ(r) vanishes at convergence, and stable algorithms were found with either a constant update or the Barzilai-Borwein method.The update function can also use the inverse electron density.Source: A. The optimized effective potential… · p. 16
• KLI is often accurate, efficient, and stable, but its diagonal response approximation misses polarization and can violate translational invariance in light-matter ground states.The authors suggest referencing the dipole operator to the electronic center of charge to reduce permanent-dipole dependence and improve stability.Source: A. The optimized effective potential… · pp. 16–17
• For weakly coupled sodium dimers, KLI is close to OEP; the implementation also compares convergence using a constant update and the Barzilai-Borwein method.Figure 3 presents the dimer density, OEP–KLI comparison, and convergence behavior.Source: Figure 3 · p. 15, A. The optimized effective potential… · p. 17

• Dressed RDMFT maps N electrons in d dimensions plus one mode onto N dressed fermions in d+1 dimensions, preserving the same physics.The auxiliary system uses dressed coordinates z = (r, q).Source: §IV · p. 18
• The auxiliary Hamiltonian has only one- and two-body dressed-coordinate terms, making standard electronic-structure methods applicable; Octopus uses it for dressed RDMFT and HF.Source: §IV · p. 18
• The dressed 1RDM lets the electronic RDMFT total-energy functional be replaced with auxiliary-Hamiltonian terms, transferring common 1RDM-based two-body approximations to the dressed system.Source: §IV · p. 18
• At increasing coupling, dressed RDMFT stays very close to exact total energy while dressed HF deviations grow, but both miss the strongly increasing photon number.The comparison is for stretched H2 in the dressed-orbital description.Source: Figure 4 · p. 19, §IV · p. 20
• Dressed RDMFT improves considerably over dressed HF for energy, while both methods capture photon number similarly.The figure compares differences from exact ground-state energies and photon number as a function of coupling for stretched H2.Source: Figure 4 · p. 19, §IV · p. 19
• The approach currently neglects the auxiliary wave function’s exchange symmetry with respect to auxiliary coordinates.Source: §IV · p. 19

V. TOWARDS DYNAMICS OF STRONGLY CORRELATED SYSTEMS: THE TDDFT+U FUNCTIONAL

• DFT+U adds a mean-field Hubbard model on a localized subspace to improve local electron interactions, but its double-counting term is not known generally.Source: §V · p. 20
• Octopus extends DFT+U to Wannier-like localized states and intersite interactions, and implements self-consistent ab initio U and J through time-dependent ACBN0.Source: §V · p. 21
• TDDFT+U simulations reproduce the absorption spectra of transition-metal oxides including NiO and MnO.The paper reports this result for its TDDFT+U simulations.Source: §V · p. 21
• Figure 5 shows that the effective Hubbard U for Ni 3d orbitals in light-driven NiO varies over time, so it cannot be assumed constant out of equilibrium.The figure plots self-consistent Hubbard U dynamics alongside the driving vector potential.Source: Figure 5 · p. 22, §V · p. 22
• Light-driven tuning of effective electronic parameters could enable phase transitions, including light-induced magnetic Weyl semimetals.Source: §V · p. 22

VI. VAN DER WAALS INTERACTIONS

• Because van der Waals interactions are nonlocal and long-range, local and semilocal functionals cannot describe them; Octopus implements vdW-DF, vdW-TS, and vdW-D3 for isolated and periodic systems.Source: §VI · pp. 22–23

A. vdW-DF

• Octopus evaluates vdW-DF nonlocal correlation using the Roman-Pérez–Soler convolution approach, reducing the scaling from O(N^2) to O(N log N).The method expresses the integrand as a convolution using helper functions and the Fourier convolution theorem.Source: A. vdW-DF · pp. 23–24
• Octopus supports vdW-DF1, vdW-DF2, and vdW-DF-cx, plus variants including vdW-DF-optPBE, vdW-DF-optB88, and vdW-DF-C09.Source: A. vdW-DF · p. 24

B. vdW-TS

• A small van der Waals-induced red shift appears in the LDA optical spectrum of the hydrogen fluoride dimer, demonstrating TS-vdW correction in TDDFT absorption calculations.The dimer serves as a proof-of-concept application of the modular TS-vdW functional correction.Source: B. vdW-TS · pp. 24–25, Figure 6 · p. 25

C. vdW-D3

• Octopus supports DFT-D3 [66], a van der Waals correction dependent on atomic positions rather than density, validated against library reference data for consistency with other codes.The implementation links to the authors’ library, distributed with Octopus.Source: C. vdW-D3 · pp. 25–26

VII. POLARIZABLE CONTINUUM MODEL

• PCM models solvent as an infinite dielectric continuum around a molecular cavity; apparent surface charges represent the reaction potential, which Octopus computes using IEF-PCM.The dielectric response depends on cavity geometry and the solvent dielectric function.Source: §VII · p. 26
• Octopus builds PCM cavities with GEPOL and smooths discretized polarization charges with normalized spherical Gaussians to regularize Coulomb singularities between surface tesserae and real-space grid points.The independent surface and real-space grids can cause numerical problems when tesserae and grid points are close.Source: §VII · pp. 26–27
• After 10 of 19 Kohn–Sham iterations, nitrobenzene’s solvation free energy stabilizes, while total polarization charge differs from valence-electron count by only 0.03%.The example is nitrobenzene in water with static dielectric constant ϵs = 80.Source: §VII · pp. 27–28, Figure 7 · p. 28
• TD-PCM offers equilibrium, inertial, and EOM approaches: equilibrium assumes instantaneous solvent response, inertial separates fast and slow response, and EOM evolves polarization charges with full history dependence.The inertial approximation is suited to solvents whose relaxation times are large relative to solute electronic excitations.Source: §VII · pp. 28–29
• All TD-PCM schemes shift nitrobenzene absorption features to lower excitation energies than in vacuo, with non-rigid shifts that depend on excited state and solvation scheme.Figure 8 compares the gas-phase and solvated spectra.Source: §VII · p. 29, Figure 8 · p. 30
• For water’s 3.37 ps relaxation time, EOM and inertial TD-PCM spectra are almost coincident, making inertial TD-PCM the computationally preferable choice.Within Debye’s model, inertial TD-PCM is the infinite-relaxation-time limit and EOM interpolates between it and the zero-time equilibrium limit.Source: §VII · p. 29
• For a water-like solvent with 1 fs relaxation, EOM shifts approach equilibrium TD-PCM, while peak intensities lie between inertial and equilibrium results.This illustrates how EOM departs from inertial results toward equilibrium as solvent relaxation becomes faster.Source: §VII · pp. 29–30
• Reaction-field effects mainly account for peak redshifts, while cavity fields directly affect absorption intensities and can change relative peak intensities under EOM non-equilibrium solvation.The cavity field reflects polarization induced by the external field; its effect depends on the effective local field allowed by dielectric properties and cavity geometry.Source: §VII · p. 31
• Octopus can include an implicit dielectric continuum in time-domain quantum-mechanical excited-state calculations and select among TD-PCM variants to capture different solvent relaxation regimes.The three variants provide a range of relaxation treatments.Source: §VII · p. 31

VIII. MAGNONS FROM REAL-TIME TDDFT

• An infinitely short transverse Zeeman kick excites magnetic fluctuations; time-evolving the spin magnetization and analyzing it in Fourier space yields spin susceptibilities under the ground-state linear-response assumption.The method computes χ+−(q; ω) and a corresponding χ−+(q; ω).Source: §VIII · p. 32
• For finite momenta, the generalized Bloch theorem avoids costly large supercells, but the implemented real-time magnon boundary conditions apply only without spin-orbit coupling.Large supercells can be computationally expensive when a few meV energy resolution is needed.Source: §VIII · pp. 32–33
• The generalized Bloch and supercell simulations agree up to numerical precision, while results for cubic Ni, Fe, and Co are in very good agreement with previous works 93.Figure 9 displays the calculated bulk-Ni spin susceptibility along ΓX.Source: §VIII · p. 33, Figure 9 · p. 34
• Beyond weak kicks, the real-time approach can investigate strong-field nonlinear phenomena, excited-state kicks, and coupling to phonons or photons without new theory or code development.Weak magnetic kicks recover linear-response results.Source: §VIII · pp. 33–34
• Real-space time propagation scales linearly with the number of states rather than the quadratic scaling of sum-over-states approaches, and requires an exchange-correlation potential but no kernel.The absence of a kernel is useful for testing functionals and theory levels whose kernel expressions are complicated.Source: §VIII · p. 34

IX. ORBITAL MAGNETO-OPTICAL RESPONSE OF SOLIDS AND MOLECULES FROM A STERNHEIMER APPROACH

• For periodic solids, density-matrix perturbation theory treats uniform magnetic fields with periodic boundary conditions and automatically preserves gauge invariance (113,117,119).The formalism uses the gauge-invariant counterpart of operators and solves the response with the Sternheimer approach; local-field effects require self-consistent solution.Source: §IX · pp. 35–36
• The magnetic-field-dependent polarizability is obtained from the current response; under the “2n + 1” theorem, the second-order derivative ρ(EµBγ) need not be calculated explicitly.124,125The dielectric tensor in the magnetic field is ϵνµ = δνµ + 4πανµ/w, with w the unit-cell volume.Source: §IX · pp. 36–37
• In the large-supercell limit, the periodic formalism becomes equivalent to the standard finite-system formulation, which uses electric and orbital magnetic dipole moments (113).For finite systems with real wavefunctions and time-reversal symmetry, the supplementary electric-field response need not be solved separately.Source: §IX · p. 37
• For silicon, excitonic effects improve agreement with experiment, while calculated magneto-optical peak magnitudes remain about two times smaller than measurements (107).Without excitonic effects, Re/Im ϵxy is already qualitatively similar to experiment at the direct absorption edge; reducing the assumed linewidth δ can correct the peak magnitudes.Source: §IX · p. 37, Figure 10 · p. 38
• With self-consistent local-field effects, magneto-optical spectra for solids take twice as long to calculate as polarizability, while finite-system calculations require the same effort (113).The scheme’s efficiency is otherwise comparable to standard linear-response calculations of simple optical polarizability without a magnetic field.Source: §IX · p. 38

X. TIME-DEPENDENT ANGULAR RESOLVED PHOTOELECTRON SPECTROSCOPY

• The implementation yields momentum-resolved probability P(p), from which energy spectra P(E) and angle-resolved spectra P(p∥, E) can be derived.It supports solids as well as atoms and molecules; the formalism can accommodate non-collinear spin, though the presentation restricts itself to closed-shell collinear systems.Source: §X · p. 39
• t-SURFF computes photoelectron probability from the time-integrated current flux through a surface separating inner Kohn-Sham dynamics from outer Volkov-electron dynamics.The method assumes the outer-region electrons can be treated as independent particles driven by the external field.Source: §X · pp. 39–40
• The flux surface can be chosen as a plane for periodic systems or a sphere for non-periodic systems, but calculations require explicit vacuum and, for solids, a surface slab.The surface position must be converged at increasing distances; this adds cost relative to bulk unit-cell calculations, while including surface-layer signals relevant to experiments.Source: §X · pp. 40, 42
• For strong-field ionization of C60, electrons rescattering at the same time form a final-momentum ring centered at the vector potential at rescattering.The angular distribution is shown for an infrared field polarized along x.Source: §X · p. 41, Figure 11 · p. 41
• Apart from the larger simulation box needed for absorbing boundaries, the method is not computationally costly and supports existing k-point, state, and mesh parallelizations.Absorbing boundaries prevent spurious reflection of electrons from the simulation-box edges.Source: §X · p. 42
• Arbitrary linear combinations of pulses enable pump-probe configurations, and the method also supports complex techniques such as RABBITT.Real-time excitation need not come from an optical pulse, as the vector potential may take an unrestricted functional form.Source: §X · p. 42
• In simulated pump-probe ARPES of monolayer h-BN driven resonantly at K, excited-state population transfer to the conduction band is visible in the spectrum.Figure 11(c) shows the simulated measurement.Source: Figure 11 · p. 41, §X · p. 42
• TDDFT-Ehrenfest dynamics can simulate electron-phonon coupling signatures in ARPES [136,137].Source: §X · p. 43
• Pump-probe simulations can study steady-state modifications of driven electronic states and Floquet physics ab initio, including periodic-force effects on dressed electronic structure [138,139].Source: §X · p. 43
• With t-SURFF, electrons can flow back and forth through the surface without artifacts in the final spectrum.The method accumulates flux through the surface over the propagation time.Source: §X · p. 43
• C60 ionized by a strong IR pulse shows characteristic photoelectron rings centered at the vector potential at rescattering, and Octopus simulations agree excellently with experiments [141,142,143].The pulse induces rescattering dynamics imprinted in the photoelectron angular distribution.Source: §X · p. 43

XI. ELECTRIC AND THERMAL CONDUCTIVITIES

• Octopus tracks current and heat-current densities at each simulation step after an electric-field impulse introduced through a vector potential compatible with periodic boundaries.The field is E(t) = E0δ(t), and the vector-potential approach satisfies the periodic boundary conditions of an extended system.Source: §XI · p. 44
• Fourier transforming the simulated currents yields frequency-dependent electrical and thermal conductivities, with DC conductivity evaluated at ω = 0.The frequency-dependent conductivities are shown for the hydrogen example in Fig. 12.Source: §XI · pp. 44–46
• For liquid metallic hydrogen at 1400 K and 400 GPa, the time-dependent current density for one ionic snapshot of a 128-atom supercell after an applied electric field decays to zero; frequency-dependent conductivities are shown in Fig. 12.Source: §XI · pp. 44–45, Figure 12 · p. 45
• The illustrated heat current includes only the kinetic-energy contribution; the potential-energy and electron-interaction contributions to the energy current are not yet implemented.The authors identify evaluation of these remaining terms as future work to determine whether they are significant or can be neglected for some systems.Source: §XI · p. 45
• Unlike linear-response Kubo–Greenwood theory, the TDDFT method can study nonlinear conduction, as demonstrated for liquid aluminum [146].The passage describes that nonlinear conductivity effects in liquid aluminum could not be studied with standard linear-response theory.Source: §XI · p. 46
• The authors plan to extend the tools with a thermal vector potential to induce heat currents in extended systems during time-dependent simulations.Source: §XI · p. 46

XII. LOCAL DOMAIN CONTRIBUTION TO PHYSICAL OBSERVABLES

• Density-functional theory makes observables functionals of electronic density; post-processing partitions that density among regions, while additive operators such as the time-dependent dipole support local contributions.The dipole is relevant to optical absorption cross sections.Source: §XII · pp. 46–47
• With non-overlapping local domains, density-dependent (semi-)local operators satisfy the additivity condition; Octopus offers geometric and atom-dependent regions and Bader volumes with zero-flux boundaries.Overlapping regions require care when comparing domain results.Source: §XII · p. 47
• For the benzene-fulvene dimer, Bader atomic decomposition defines molecular domains as sums of atomic volumes; Fourier-transforming each domain’s local time-dependent dipole moment recovers its local polarizability tensor.Source: §XII · p. 48
• Local analysis attributes the major spectral change to reduced oscillator strength of the fulvene peak around 6.5 eV, in agreement with FDE-rt-TDDFT [152].The cited comparison is presented as validation of the local-domain strategy.Source: §XII · pp. 48–49
• At 8 Å the dimer shows peaks at isolated-molecule frequencies; closer spacing red-shifts the ~5 eV peak and weakens the 6.5–7 eV peak, transferring oscillator strength higher.Figure 13 presents the spectral decomposition as a function of intermolecular separation.Source: §XII · p. 48, Figure 13 · p. 49
• The method requires no prior fragmentation, basis selection, or localization, can treat large systems such as LHC-II, and enables exciton-coupling calculations with real-time TDDFT transition densities [158].Source: §XII · p. 49

XIII. NEW PROPAGATORS FOR REAL-TIME TDDFT

• For adiabatic-approximation TDKS, symplecticity is a relevant property; symplectic propagators better conserve energy over long times, with energy oscillating around its true value rather than diverging.Source: §XIII · pp. 50–51
• Tests of multistep, exponential Runge–Kutta, and commutator-free Magnus families found CFM suitable for TDDFT based on accuracy, stability, and performance, leading to Octopus implementation.Source: §XIII · p. 51
• Octopus implements fourth-order CFM4 by extrapolating Hamiltonian components from prior steps, yielding an explicit scheme that preserves fourth-order accuracy when extrapolation is also fourth order.The extrapolated components are the Hartree, exchange, and correlation parts; the explicit method requires no linear or nonlinear algebraic solves.Source: §XIII · p. 52
• CFM4 outperforms EMR at every examined error, with clearer gains at higher accuracy; its fourth-order error slope contrasts with EMR’s second-order slope.The benchmark compares propagation cost with accuracy for a benzene molecule after an instantaneous perturbation; CFM4 requires two exponentials versus EMR’s one.Source: §XIII · p. 52, Figure 14 · p. 53

XIV. CONJUGATE GRADIENT IMPLEMENTATION IN RDMFT

• The conjugate-gradient method converges with fewer natural orbitals and to a slightly lower total energy than the tested Piris basis-set implementations, exploiting grid flexibility beyond those bases.The authors attribute the lower energy to contributions to the natural orbitals not covered by the available basis sets.Source: Figure 15 · p. 55, §XIV · p. 57
• Among tested Piris-method starting bases, independent-particle orbitals gave the fastest convergence and lowest energies; extra localization helped only few orbitals and slowed convergence with basis size.Independent-particle orbitals also yielded the largest number of localized orbitals among the tested bases, so the added localization step was removed.Source: §XIV · p. 55
• The conjugate-gradient implementation optimizes orbitals directly on real-space grids, supports systematic refinement by enlarging the grid or reducing spacing, and adapts the DFT procedure to RDMFT orthogonality constraints.Unlike the Piris method, it can start from random states and does not require the initial orbitals to serve as a basis; non-diagonal Lagrange multipliers modify the descent direction and line minimization.Source: §XIV · p. 56
• For the Müller functional, initial-state quality affected iteration count but not the converged result; finite orbital counts and possible local minima mean other functionals may require start-point convergence checks.Although the Müller functional is convex with infinitely many orbitals, practical calculations use finite sets, and not all RDMFT functionals are known to share that property.Source: §XIV · pp. 56–57

XV. PERIODIC SYSTEMS AND SYMMETRIES

• Octopus generates real-space grid points along primitive axes, calculates finite-difference gradient and Laplacian weights, and samples reciprocal space through a k-point grid.Periodic calculations use the primitive cell under Bloch’s theorem and sample the Brillouin zone.Source: §XV · p. 58
• Crystal symmetries obtained with spglib reduce the Brillouin zone to its irreducible zone, but the implementation currently includes only symmorphic symmetries.Non-symmorphic symmetries involving fractional translations are deferred because they are incompatible with arbitrary real-space grids.Source: §XV · p. 58
• 1.50 meV/atom was obtained in the Delta-factor test using Schlipf-Gygi ONCVPSP 2015 pseudopotentials, a value reported as very close to those from other codes.The comparison is described as using the same pseudopotential set.Source: §XV · p. 58
• For symmetry-breaking perturbations, Octopus uses the original space group’s small group, retaining only symmetries that leave the perturbation direction invariant for time-dependent calculations.The described perturbations include vector potentials, finite-momentum kicks, and strain.Source: §XV · p. 58
• Real-space symmetrization of charge and current densities and other observables is important for obtaining good, stable numerical results.It complements Brillouin-zone reduction to the irreducible zone.Source: §XV · p. 58
• Periodic ion dynamics combines Newton equations for ionic reduced coordinates with motion of the primitive cell, and Octopus implements stress-tensor calculations and lattice dynamics.The stress tensor is identified as a key ingredient of the lattice equation of motion.Source: §XV · p. 59

XVI. ADDITIONAL TECHNICAL CODE IMPROVEMENTS

• The section’s technical developments target code reliability and efficiency, easier maintenance, and challenges posed by computer architectures.The authors identify the code’s scale and complexity as reasons for ongoing optimization, validation, and code-quality work.Source: §XVI · p. 59

A. Web application for analyzing regression tests

• The regression suite spans 30 toolchains and about 180 tests, runs Octopus roughly 740 times, makes about 12,000 reference comparisons, and has about 71% code coverage.Buildbot runs it automatically as part of continuous integration.Source: A. Web application for analyzing… · pp. 59–60
• An interactive web application helps diagnose failed regression tests, distinguish toolchain-specific problems from numerical variation, and has already helped identify regression-causing bugs.It also facilitates updating tests and improving the testsuite.Source: A. Web application for analyzing… · p. 60

B. Improving ground-state calculations

• The revised default conjugate-gradients eigensolver improves solver reliability and, by orthogonalizing against lower-energy bands, converges faster in most tested cases.Its implementation follows Ref. 167 closely, with the stated orthogonalization difference.Source: B. Improving ground-state calculations · p. 60
• The real-space low-pass filter preconditioner was most effective across varied systems, but its optimal damping depends on the system and adding Jacobi iterations did not improve it.The filter applies a weighted combination of a grid value and its neighboring values; it was compared across molecules, semiconductors, metals, and surfaces.Source: B. Improving ground-state calculations · p. 61
• The filter preconditioner most effectively reduced total ground-state calculation time by cutting SCF iterations without excessive computational cost.The compared multigrid method reduced SCF iterations but was less effective because it applied the Laplacian several times per iteration.Source: B. Improving ground-state calculations · pp. 61–62
• Some states, especially in calculations with many unoccupied states, may not fully converge with a preconditioner and require restarting without one.The reasons remain under investigation.Source: B. Improving ground-state calculations · p. 62

C. Novel multi-system framework

• Octopus is being refactored to represent the physical system as interacting subsystems, such as electrons, nuclei, and Maxwell fields.Motivated by Maxwell–Kohn–Sham coupling, the generic framework can accommodate existing features and future developments requiring coupled subsystems.Source: C. Novel multi-system framework · p. 62

D. Memory layout

• The new orbital layout stores states in batches with the state index innermost, enabling contiguous memory access and identical grid-loop operations across states.This supports GPU warp parallelization and CPU vectorization; finite-difference kernels such as the Laplacian use SSE, AVX, or AVX512 instructions.Source: D. Memory layout · p. 63

E. GPUs

• Octopus added CUDA support, multiple devices per host, and custom GPU memory management; packed storage can keep states in GPU memory and reduce transfers.Removing frequent temporary-variable allocations and deallocations improved scaling across multiple GPUs and GPU-equipped nodes.Source: E. GPUs · p. 63
• The GPU implementation was 4.8× faster than one full CPU node, reducing time to solution by factor 4.8; despite drawing ca. 950 W versus 450 W, it cut energy to solution by factor 2.3.Source: E. GPUs · p. 64
• The Octopus GPU implementation, scaling from one GPU node to 2, 4, and 8 nodes, achieved speed-ups of 1.26, 2.14, and 3.37, with corresponding parallel efficiencies of 63%, 54%, and 42%, respectively.Source: Figure 16 · p. 64, E. GPUs · p. 64
• For Octopus’s GPU implementation, speed-ups for a simulation using 2, 4, and 8 GPUs were 1.95, 3.70, and 5.97, corresponding to parallel efficiencies of 97%, 93%, and 75%.Source: Figure 16 · p. 64, E. GPUs · p. 64

XVII. CONCLUSIONS

• The conclusions emphasize Octopus’s real-space flexibility and new approaches for non-equilibrium electron, phonon, and photon dynamics, including Maxwell-Kohn-Sham, OEP, and dressed RDMFT methods.The framework is presented as a setting for implementing and testing ideas that may later be adopted by other codes.Source: §XVII · p. 65
• The conclusions also highlight real-time magnons and Sternheimer orbital magneto-optical response, with scalable real-space grids making supercell studies promising for complex correlated magnets.The paper additionally reports efficient treatment of molecules, nanoparticles, solvents, solids, and monolayers, and periodic systems without loss of accuracy in the Delta-factor test.Source: §XVII · pp. 65–66
• Planned extensions include Floquet and cavity-QED materials engineering, multitrajectory electron-ion dynamics, open quantum dissipative systems, and spectroscopy with entangled photons.The paper also emphasizes recent performance and scalability improvements, particularly GPU support.Source: §XVII · p. 66