Algorithm 2 sketches our main addition to Psi4, the function get_solvation_contributions. In a ground state calculation this function is called in each SCF step to add to the Fock or Kohn-Sham matrix the solvation contribution which is computed from the current density matrix. The first block of Algorithm 2 performs the basic setup, which only needs to be done once, e.g. in the first SCF iteration. This step processes the user-defined parameters discussed in Section 3.2 to initialize the ddX model (pyddx.Model) as well as the Psi4-specific helper classes for the required analytical (MintsHelper) and numerical (NumIntHelper) integrals. Based on these helper classes, the electronic contributions to F and Ψ (Step 3a in Section 3.1) as well as the solute contributions to the derivatives (Step 4) will be computed on the host side. For the detailed expressions, we refer the reader to the aforementioned method-specific literature . Next, the nuclear contributions are added by calling model.multipole_electrostatics() and model.multipole_psi(). If this is the first iteration of the selfconsistent field loop, the resulting F and Ψ are used to initialize the pyddx.State. This object in particular stores the RHS of the linear systems ( ) and ( ) as well as iterated solution vectors X and X adj . A following call to state.fill_guess() and state.fill_guess_adjoint() generates a crude initial guess for X and X adj . As discussed in Section 3.4.1, in any of the subsequent SCF iterations we only update the RHS of the ddX equations (state.update_problem()), and keep the solution from the previous SCF step to be used as a guess for the ddX iterative solver. This reduces the required number of iterations to solve the updated linear systems ( ) and (9) in particular in the later SCF steps, where the density matrix only changes little. Next, the forward and adjoint problems defined inside the state object are solved and their solutions X and X adj employed to compute the Fock matrix contributions -using again the host-specific integration routines. Finally, we consider the computation of solvated excited states in a linear-response TDDFT formalism. Psi4 already offers the method set_external_cpscf_perturbation to add further potential terms to the response equations. We use this function to add a call to get_solvation_contributions when computing vertical excitation energies. For this, we follow the standard approach to only consider the fast (electronic) components Algorithm 2: Pseudocode illustrating the integration of ddX within Psi4. The function takes a density matrix (or transition density matrix) and returns the respective Fock matrix contribution. Psi4 specific code is highlighted. in the solvent's response, which in Algorithm 2 reflects in two notable changes. First, we employ the optical permittivtiy ε ∞ instead of the static permittivity ε s when setting up the pyddx.Model and second, we skip over lines 14 and 15, i.e., we do not add the nuclear contribution to F and Ψ.