S O F TWARE F OCU S ddX: Polarizable continuum solvation from small molecules to proteins Michele Nottoli1 | Michael F. Herbst2 | Aleksandr Mikhalev3 | Abhinav Jha1 | Filippo Lipparini4 | Benjamin Stamm1 1Universität Stuttgart, Institute of Applied Analysis and Numerical Simulation, Stuttgart, Germany 2Mathematics for Materials Modelling, Institute of Mathematics & Institute of Materials, �Ecole Polytechnique Fédérale de Lausanne, Lausanne, Switzerland 3Center for Artificial Intelligence Technology, Skolkovo Institute of Science and Technology, Moscow, Russia 4Dipartimento di Chimica e Chimica Industriale, Universit�a di Pisa, Pisa, Italy Correspondence Michael F. Herbst, Mathematics for Materials Modelling, Institute of Mathematics & Institute of Materials, �Ecole Polytechnique Fédérale de Lausanne, CH-1015 Lausanne, Switzerland. Email: michael.herbst@epfl.ch Filippo Lipparini, Dipartimento di Chimica e Chimica Industriale, Universit�a di Pisa, 56124 Pisa, Italy. Email: filippo.lipparini@unipi.it Benjamin Stamm, Universität Stuttgart, Institute of Applied Analysis and Numerical Simulation, Pfaffenwaldring 57, 70569 Stuttgart, Germany. Email: benjamin.stamm@ians.uni- stuttgart.de Funding information Deutsche Forschungsgemeinschaft, Grant/Award Numbers: 440641818, EXC2075-390740016, INST 40/575-1 FUGG Edited by: Peter R. Schreiner, Editor-in- Chief Abstract Polarizable continuum solvation models are popular in both, quantum chemistry and in biophysics, though typically with different requirements for the numerical methods. However, the recent trend of multiscale modeling can be expected to blur field-specific differences. In this regard, numerical methods based on domain decomposition (dd) have been demonstrated to be sufficiently flexible to be applied all across these levels of theory while remaining systematically accurate and effi- cient. In this contribution, we present ddX, an open-source implementation of dd- methods for various solvation models, which features a uniform interface with clas- sical as well as quantum descriptions of the solute, or any hybrid versions thereof. We explain the key concepts of the library design and its application program inter- face, and demonstrate the use of ddX for integrating into standard chemistry pack- ages. Numerical tests illustrate the performance of ddX and its interfaces. This article is categorized under: Software > Quantum Chemistry Software > Simulation Methods KEYWORD S biophysics, continuum solvation, COSMO, linearized Poisson–Boltzmann, polarizable continuum model, quantum chemistry 1 | INTRODUCTION Polarizable continuum solvation models (PCSM) are powerful, simple and effective tools that are used to describe solva- tion effects on molecules and larger biological systems. Most commonly, such models are used by two large Received: 13 March 2024 Revised: 18 June 2024 Accepted: 8 July 2024 DOI: 10.1002/wcms.1726 This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. © 2024 The Author(s). WIREs Computational Molecular Science published by Wiley Periodicals LLC. WIREs Comput Mol Sci. 2024;14:e1726. wires.wiley.com/compmolsci 1 of 24 https://doi.org/10.1002/wcms.1726 https://orcid.org/0000-0002-6544-0897 https://orcid.org/0000-0003-0378-7921 https://orcid.org/0000-0002-4947-3912 https://orcid.org/0000-0003-3375-483X mailto:michael.herbst@epfl.ch mailto:filippo.lipparini@unipi.it mailto:benjamin.stamm@ians.uni-stuttgart.de mailto:benjamin.stamm@ians.uni-stuttgart.de http://creativecommons.org/licenses/by/4.0/ http://wires.wiley.com/compmolsci https://doi.org/10.1002/wcms.1726 communities: quantum chemists employ them as a cost-effective model to include environmental effects into the pre- diction of molecular geometries and spectroscopic observables.1–6 Biophysicists use them to study the electrostatics of large solvated biomolecules, and hence their binding properties.7–9 Both communities have developed PCSM indepen- dently and coded in different software and with different numerical strategies; the models all share a common ingredi- ent, that is, they treat the solvent as a continuum, infinite, structureless dielectric that is polarized by the solvent and, possibly, back-polarizes it. The electrostatic response of the dielectric is obtained by solving either the generalized Poisson equation or, for (weakly) ionic solutions, the (linearized) Poisson–Boltzmann equation. Despite this common theoretical background, the numerical realizations of PCSM implemented by the two communities are indeed quite different. In the context of quantum chemistry, PCSM have been developed to achieve quantum mechanical calculations of molecules in solution at very little additional cost compared to the in vacuo setting. The most common numerical strat- egy is based on a boundary element method approach, where the generalized Poisson equation is rewritten as an inte- gral equation for an apparent surface charge supported at the boundary of the molecular cavity that accommodates the solute. The resulting linear equations are usually solved by dense linear algebra techniques, such as LU or Cholesky decomposition. This approach scales cubically and thus in principle expensive. However, its cost is still negligible with respect to the cost associated with the quantum mechanical treatment. On the other hand, as PCSM are typically used to compute structures and properties, a lot of work has been done over the past decades to achieve smooth, differentia- ble definitions of the solvation energy. These can be used to compute analytical derivatives of the molecular energy, for example, to compute response properties or derivatives with respect to the nuclear positions (forces and harmonic fre- quencies). In quantum chemistry, while more sophisticated models like the linearized Poisson–Boltzmann exist,10,11 the two most popular models are based either on a dielectric continuum (polarizable continuum model, PCM)2,4,12,13 or on a conductor continuum (Conductor-like Screening Model, COSMO, also known as C-PCM).14–16 Such models have been implemented in virtually every quantum chemistry package and are available in conjunction with most levels of theory, from semiempirical methods17–20 to sophisticated post-Hartree–Fock ones21–25 and also with a linear scaling implementation.26 In the context of biophysics, on the other hand, modeling ionic solutions is paramount, which makes the (linear- ized) Poisson–Boltzmann ((L)PB) model27–29 the de facto standard. Many numerical realizations of such a model have been implemented and are available to the community in popular software that solves the (L)PB equation using a vari- ety of methods that include finite differences,30–33 finite elements,34–36 and boundary elements.37–47 As the modeled sys- tems are typically large biomolecules described with a classical force field (i.e., point charges, or sometimes more sophisticated treatments including distributed multipoles and polarizabilities) and not with quantum mechanics, the treatment of solvation can become a cost-dominating factor in such computations. Therefore, much attention has been given to the numerical efficiency of the implementations, and linear scaling techniques are commonly used to afford the treatment of very large systems. Furthermore, due to the complexity of the systems studied, the use of accurate molecular surfaces, such as the Connolly also known as Solvent Excluded Surface (SES) to define the molecular cavity is common, to be compared with the simpler Van der Waals or Solvent Accessible Surfaces (SAS) used in quantum chemistry applications. On the other hand, the calculation of analytical gradients and the coupling to quantum mechanical codes have not received the same amount of attention as the models used in quantum chemistry. Despite the differences, there are several regions of overlap, from PCSM being applied to large systems treated with multiscale QM/MM methods,48–57 to the possibility of performing Molecular Dynamics simulations in a polarizable continuum solvent,58,59 that makes the possibility of bridging the gap between the two communities interesting. In the last decade, a new numerical paradigm for continuum solvation models has been proposed. The new algo- rithm, based on Schwarz's alternating domain decomposition (dd), has originally been proposed by some of us for the COSMO model. The resulting algorithm, which we have named ddCOSMO,60 exhibits all the properties that make it an ideal candidate to bridge the gap between the two communities. ddCOSMO achieves linear scaling in computational cost and memory requirements, can be used to compute analytical gradients61 and other response properties with arbi- trary accuracy, it can be used to describe very large systems,62 and can be easily interfaced with quantum mechanical56 and polarizable models.49,55,58 After the initial work on ddCOSMO, the domain decomposition paradigm has been extended to the PCM63–65 and, more recently, to the LPB model.66 For the latter models, a straightforward dd imple- mentation exhibits quadratic scaling computational cost. To overcome this difficulty, a linear scaling implementation based on the Fast Multiple Method has further been achieved.67,68 All three models have analytical gradients and are numerically robust and accurate. 2 of 24 NOTTOLI ET AL. In this contribution, we present ddX, an open-source implementation of ddCOSMO, ddPCM, and ddLPB that comes with clear application program interfaces (APIs) and is designed to be easily interfaced with both classical and quantum descriptions. We believe that the ddX library, besides being the ddX coronation of more than 10 years of methodologi- cal developments, can provide both quantum chemists and biophysicists with a unified tool to treat solvation effects with a polarizable model. It can be used consistently for any application, ranging from the simulation of electronic spec- tra to molecular dynamics, from quantum mechanical geometry optimizations to the calculation of the electrostatic sol- vation energy of large biomolecules. The library is distributed on GitHub (https://github.com/ddsolvation/ddX) under the permissive LGPL v3 license and has been provided with an automatically generated documentation, available at https://ddsolvation.github.io/ddX/. So far we have interfaced ddX both with codes from quantum chemistry (Psi469 and adcc70) as well as classical MD simulation (Tinker71), which will be demonstrated below. This paper is organized as follows: in Section 2 we introduce the models available in ddX and give a brief overview on the computation of energy and the analytical derivatives. In Section 3, we present the layout of the implementation of the ddX library along with some details about the fast multipole method (FMM) and our code design. In Section 4, we present our interface to Tinker, Psi4 and adcc in detail. Lastly, in Section 5, we test the ddX library on molecules with atoms up to 105 atoms. The supporting information (Data S1) provides the finer details about the matrix formula- tion for the models. 2 | METHODS 2.1 | Model input: Solute cavity and density The ddX-library considers solute cavities that can be written as a union of balls, thus including van der Waals (vdW)- and Solvent Accessible Surface (SAS)-cavities.72,73 The domain-decomposition approach for Solvent-Excluded Surfaces (SES)-cavities has also been developed,74 but this does not (yet) fit the ddX-framework. We therefore consider a mole- cule consisting of N atoms, for each atom we assign a ball Ωj ≔Brj xj � � , centered at the nucleus position xj �ℝ3 with radius rj �ℝ. The region occupied by the solute molecule is then given by Ω¼ [N j¼1 Ωj: The region not occupied by the solute, that is, Ω¼ℝ3 ∖Ω, is considered to be occupied by the bulk solvent. Note that without many changes, the method can also be applied if xj do not necessarily represent nuclear coordinates or just a subset of those (such as is needed in coarse-grained models). The charge density of the solute is denoted by ρ. Within the framework of ddX, the solute molecule is entirely described by the cavity Ω and the charge distribution ρ. As it is customary in continuum solvation models, it is assumed that the charge ρ is supported in Ω, that is ρ xð Þ¼ 0 for x �ΩC. As we will see in the following, different models and thus different codes have different representations for the charge distribution ρ. Therefore, ρ can be represented either by a sum of (partial) point-charges or, more generally, a sum of multipoles for classical force-field models, or a sum of point-charges (nuclear charges) in combination with the electronic density ρel for QM-based models. In the latter case, ρel is a distribution that is assumed to be supported in Ω and is represented in terms of basis func- tions arising from the interfaced QM-code. In practical calculations, however, ρel is not supported in Ω due to the Gaussian- or Slater-type basis functions having tails that extend beyond the molecular cavity and any electronic charge on the complement of Ω is lost, a problem known as the escaped charge problem.75 This problem has been extensively studied75–77 and has been found not to be an issue in practical applications. 2.2 | Models for electrostatic interaction The ddX-library offers three levels of the theory that model the electrostatic interaction between the solute molecule and the bulk solvent. These models differ by the physical interaction between the solute and the bulk (implicit) solvent. NOTTOLI ET AL. 3 of 24 https://github.com/ddsolvation/ddX https://ddsolvation.github.io/ddX/ We start the presentation with the physically most complex one of the three and present the other two as its simplification: 1. LPB: Linearized Poisson–Boltzmann model.27–29 In the LPB model, the (implicit) solvent is assumed to be a homoge- neous ionic continuum with relative dielectric permittivity εs <∞ and Debye Hückel screening constant κs. Then, the LPB equations are given by: find the potential V that satisfies �r� εrVð Þþ κ2V ¼ 4πρ, inℝ3, ð1Þ where ε xð Þ¼ 1 if x �Ω, εs otherwise, � and κ xð Þ¼ 0 if x �Ω,ffiffiffiffi εs p κs otherwise: � In this simple model, we do not account for the Steric effects78 and hence do not include a Stern layer. 2. PCM: Polarizable Conducting Model.2,4,12,13 The PCM is the particular case of LPB with κs ¼ 0 and εs <∞. The ddPCM method relies on the following equivalent Integral Equation formulation: Find the charge density σ and the intermediate potential Φε that satisfy RεΦε ¼R∞Φ, on ∂Ω, ð2Þ S0σ¼�Φε, on ∂Ω, ð3Þ where Rε ¼ 2π εsþ1 εs�1 I�D0, and S0 and D0 denote the single and the double-layer boundary operators, respectively, with respect to the Green's kernel G0 rð Þ¼ 1 4πr. 12 3. COSMO: COnductor-like Screening Model.5,14–16 The COSMO is the particular case of LPB with εs ¼þ∞ and κs ¼ 0, that is, the solvent is assumed to be a perfect conductor. In consequence, V ¼ 0 in ΩC and the problem reduces to �ΔV ¼ 4πρ, inΩ, V ¼ 0, on ∂Ω, ð4Þ which, for W ¼V �Φ, can be reformulated as �ΔW ¼ 0, inΩ, W ¼�Φ, on ∂Ω: We are interested in the computation of the electrostatic solvation energy and the electrostatic forces of these models. In all three models, the solvation energy Es is defined as Es ¼ f εsð Þ1 2 Z Ω ρ xð ÞW xð Þdx, where W ¼V �Φ denotes the reaction potential and 4 of 24 NOTTOLI ET AL. f εsð Þ¼ 1 for PCM and LPB, εs�1 εsþ0:5 for COSMO, 8< : see reference 60. The force-vector Fi acting on the ith atom corresponding to the solvation energy is given by Fi ¼�rxiEs: ð5Þ 2.3 | Discrete equations in a general framework After discretization of the equations using spherical harmonics and domain decomposition techniques, the resulting linear system for all three methods can be written in the general form MX¼F, ð6Þ where F is a given vector, X is the vector of unknown degrees of freedom, and M is the given solution matrix. In all three cases, the entries of X allows to represent the (approximate) reaction potential W and can be used to compute the (approximate) electrostatic solvation energy that is of the form Es ¼ 1 2 Ψ,Xh i, ð7Þ where � , �h i denotes the scalar product of the two enclosed vectors and Ψ is related to the discretization of ρ. For brevity, we do not provide the matrix and vector entries of M and Ψ for each method here, and refer to the supporting information (Data S1) for a brief overview. For a detailed derivation and discussion on this topic, we refer to references 61, section 2, 63, section III, and 66, section 6.2.2 for the ddCOSMO, ddPCM, and the ddLPB method, respectively. 2.4 | Force computations and derivatives For the computation of the electrostatic forces, one needs the computation of gradients as seen by (5). The derivative of the energy with respect to any parameter λ, such as the position of xk of the k th atom, is given by Eλ s ≔ ∂Es ∂λ ¼ 1 2 Ψλ,X � �þ1 2 Ψ,Xλ � �¼ 1 2 Ψλ,X � �þ1 2 Xadj,Fλ�MλX � � , ð8Þ where Xadj is solution to the adjoint linear system M�Xadj ¼Ψ, ð9Þ and M� is the adjoint of M. Due to the role of Ψ in Equation (9), we will refer to it also as “adjoint RHS” and the deriva- tion of Equation (8) is given in the supporting information (Data S1). For the three models, the RHS F is written as the product of a term which is specific to the solute and a term which is merely geometric. This allows the factorization (due to the product rule) of its derivative in a term that does not depend on the nature of the solute (Fλ non�spec) and a term that has different expressions depending on the nature of the solute (Fλ spec). Hence the derivative is written as Fλ ¼Fλ non�specþFλ spec, and we can further write Eλ s as NOTTOLI ET AL. 5 of 24 Eλ s ¼ 1 2 Ψλ,X � �þ1 2 Xadj,Fλ spec D E þ1 2 Xadj,Fλ non�spec D E �1 2 Xadj,MλX � �¼Eλ s,specþEλ s,non�spec, ð10Þ where Eλ s,spec collects the first two derivative terms and Eλ s,non�spec collects the remaining two. Detailed matrix entries for the derivatives can be found in references 61, section 3, 64, appendix C, and 68, section A for the ddCOSMO, ddPCM, and ddLPB method, respectively. 3 | IMPLEMENTATION With the ddX library, we present a uniform implementation of domain decomposition solvation methods. While the main parts of the library are written in Fortran 2008 making use of the high-performance features of the language for scientific computation, convenient interfaces to C as well as Python are available as well. Based on this multi- language API a seamless interfacing to a variety of common quantum chemistry software is possible, which will be demonstrated in Section 4 for both the case of classical and QM codes and employing both our Fortran and our Python API. In this section, we will provide an overview of the design of ddX and its API, following along the steps of a typical ddX calculation. Afterward we provide some notes on the specific fast multipole (FMM) acceleration, which enabled us to achieve linear scaling for all implemented solvation models, as well as some notes on our means of distributing ddX. 3.1 | Steps of a calculation The typical steps of a ddX calculation for computing the energy and generic analytical derivatives is as follows: 1. Initialize the library. 2. Compute the primal and adjoint RHSs F and Ψ. 3. Solve the problem a. Setup the problem by loading the RHSs F and Ψ. b. Generate a guess for X and then solve the primal linear system to get X (Equation 6). c. Compute the energy as Es ¼ 1 2 Ψ,Xh i (Equation 7). d. If the calculation of an analytical derivative is requested, generate a guess for Xadj and then solve the adjoint lin- ear system to get Xadj (Equation 9). e. Compute the non-specific contribution to the derivative using X, Xadj, Fλ non�spec and Mλ (last two terms Equation 10). 4. Compute the solute specific contributions to the derivative using Ψλ, Fλ spec and Xadj (first two terms Equation 10). 5. Finalize the library. Notably, the Steps 3d, 3e, and 4 are specific for the computation of analytical derivatives. Moreover, Step 3e is non- zero only if λ is a parameter on which the continuum model explicitly depends, and specifically the most relevant case is when λ are the nuclear positions. Let us illustrate the various steps in three cases: • In a plain energy calculation for a classical solute, no analytical derivative needs to be computed, thus none of the Steps 3d, 3e and 4 are required. • A QM/ddX calculation requires the Fock matrix to solve the self-consistent field equations. This is the derivative of the energy with respect to the density matrix elements Pμν, thus Steps 3d and 4 are required. However, Step 3e is not required, as the density matrix elements are not a parameter of the continuum solvation model, and EPμν s,non�spec ¼ ∂Es,non�spec ∂Pμν ¼ 0. • Finally, an energy and force calculation for a classical solute requires all the Steps 3d, 3e and 4. 6 of 24 NOTTOLI ET AL. Depending on the nature of the solute density, the details of the computation in Steps 2 and 4 differ; for example, for solutes consisting of point charges, point dipoles or a quantum density each have different expressions to obtain F or Ψ in Step 2 or the respective derivative contributions from the solute. In ddX we provide implementations for these two steps in the case of solute densities based on point multipoles. This choice was mainly motivated by three factors: (i) point multipoles cover a wide range of use cases, for example classical solutes described by a force field, or the nuclear point charges of QM solutes, (ii) point multipoles of arbitrary order can be treated in a common mathematical formalism, (iii) in ddX we have a tailored implementation of the FMM which can deal with spherical harmonics or multipoles of any order, which can thus be seamlessly employed to acceler- ate the computation in these steps. We remark that providing general implementations for other solute densities is more involved, especially for the case of quantum mechanical solutes, the expressions require the evaluation of elec- tronic integrals and are basis set dependent, see the discussion in Section 4.2. As a result, these are strongly dependent on the host code, and we have opted to not treat this step inside ddX. The following three sections provide details on how the main steps of ddX calculations are handled in the implementations. The discussion will refer to the high level API of ddX, which are reported in Table 1. See Figure 1 for a graphical representation of external codes using ddX. Both the case of a classical solute in a calculation of energy and force, and of a QM solute in a calculation of energy and Fock matrix is reported. To make the following discussion more concise, we will only mention the function names in Fortran. Table 1 can be used to infer the correspondent function names in C and in Python. 3.2 | Step 1: Initialization During the initialization step, the two high-level data structures which are used by ddX are generated. ddX works with a “model” object and a “state” object: • The model object (ddx_type in Fortran) collects all the parameters which control the calculation, like the chosen continuum solvation model, the dielectric permittivity, the information about the discretization, cavity information (the number, centers, and radii of the spheres) and similar. Furthermore, it contains all the precomputed constants which depend on the geometry of the system, and contains the preallocated workspaces, such that most of the mem- ory allocation is done during the initialization phase. This object does not contain variables which depend on the sol- ute's density, like RHSs or solutions. • The state object (ddx_state_type in Fortran) contains all the quantities that depend on a given solute density, like RHSs, solutions, and various intermediate quantities for computing analytical derivatives. TABLE 1 API functions in Fortran, C, and Python. Purpose Fortran C Python Housekeeping ddinit ddx_allocate_model Model allocate_state ddx_allocate_state State deallocate_model ddx_deallocate_model – deallocate_state ddx_deallocate_state – Functionality for multipoles allocate_electrostatics ddx_allocate_electrostatics – multipole_electrostatics ddx_multipole_electrostatics Model. multipole_electrostatics deallocate_electrostatics ddx_deallocate_electrostatics – multipole_psi ddx_multipole_psi Model.multipole_psi multipole_force_terms ddx_multipole_force_terms State. multipole_force_terms Solving the problem ddrun ddx_ddrun State.ddrun Note: We use “–” to indicate operations, which are automatic in the Python interface, such that no explicit function call is needed. NOTTOLI ET AL. 7 of 24 Using separate objects for the state and for the model makes it possible to concurrently solve different electrostatic problems characterized by the same geometry. This could be helpful within certain QM calculations of molecular response properties, which require solving a set of response equation for each perturbation, and thus to treat simulta- neously multiple perturbed densities. Note, however, that if the geometry is updated, it is necessary to recompute a new model object. The model object allocation is done by calling the subroutine ddinit in Fortran, or the corresponding functions in C and Python. This function performs three main steps: (i) a consistency check on the input parameters, (ii) the preallocation of all the workspaces, and (iii) the precomputation of the constants. The host code needs to provide all the parameters that control the calculation, of which some are mandatory, like the chosen continuum solvation model, cavity information, and the dielectric permittivity, whereas the others are optional and can be used to finely tune the method. In Fortran and Python, we used keyword arguments for the optional parameters; in this way, it is possible to implement additional functionalities that require new parameters in the future without breaking backwards compati- bility with existing external software. The state object allocation is done by calling the subroutine allocate_state, or the corresponding functions in C and Python. Allocating the state object requires a model object, as the quantities to be allocated and their size depend on the specific nature of the calculation. At this point, the state object is simply an empty container: the infor- mation about the RHSs has yet to be loaded, and the solutions have to be computed. 3.3 | Step 2: Building the RHSs for multipolar solutes As anticipated, Step 2 depends on the nature of the solute, and in ddX itself we provide support for densities written as point multipolar distributions. For this reason, we cover here only such a case. For other densities, and in particular for QM ones, some functionalities need to be added to the host code. Further information will be given in the description of the Psi4-interface, Section 4.2. The primal RHS F consists of the electrostatic potential and possibly its higher-order derivatives (field and field gra- dient, depending on the PCSM and on whether the forces are requested). To encompass all possible cases without extensive modifications to the host code, the various electrostatic properties are collected into a high level data type. In Fortran and C we provide subroutines for allocating and populating the electrostatics object (ddx_allocate_electrostatics, ddx_multipole_electrostatics). In Python, the same operations are done by the method Model.multipole_electrostatics. FIGURE 1 Interactions between a host code and the ddX library. Left: Case of a point multipolar solute density in the calculation of energy Es and forces rEs,spec and rEs,non�spec. Right: Case of a quantum mechanical density in the calculation of energy Es and Fock matrix elements Eμν s,spec. In this case, the total density ρ is split into a QM part (ρel) and a classical nuclear part (ρnu). Each quantity in the host code represents a variable, and each box in ddX represents a call to an API function. Arrows show the arguments passed to the functions and the return values. The variables model and state are passed to all subsequent functions, but this is not represented graphically for clarity. Note also that the solutions of the linear systems and many intermediate quantities are stored in the state variables, so an explicit return is missing. 8 of 24 NOTTOLI ET AL. The adjoint RHS Ψ is simpler, as it is the same regardless of the PCSM and kind of computation. In Fortran it is computed, respectively, by calling multipole_psi, or in C and Python with the corresponding functions. 3.4 | Step 3: Solving the problem Solving the electrostatic problem requires loading the RHSs in the state object, solving the primal and adjoint linear sys- tems, and possibly computing the energy and forces. We provide a high-level procedure that combines the steps together and returns the energy and, if requested, the solvent contribution to the forces. The subroutine for this step is ddrun in Fortran. This subroutine takes as input the model and state objects, the RHSs, and the tolerance for the linear systems and gives in output the energy and, if requested, the forces. An optional argument controls if the guesses for the linear sys- tems should be generated from scratch or read from the state object. This might help in the case of repeated calculations for similar RHSs, as it might be the case during the iterative solution to the self-consistent field problem. Finally, if requested, the subroutine computes the non-specific contribution to the nuclear gradients (Step 3e). 3.4.1 | Finer control Some particular applications may require a finer control over the workflow. For this reason, we decided to expose a lower-level API alongside the high-level API described so far. In particular, it is possible to run separately the Steps 3a, 3b, 3c, 3d and 3e by calling the associated subroutines. The same functionality is exposed in Fortran, C, and Python as shown in Table 2. The setup subroutine takes care of loading the RHSs in the state and setting up the calculation; this is not required in Python, as this step is done implicitly by the constructor. The subroutines fill_guess and fill_guess_adjoint can be used to generate a guess for the linear system solver; this step can be skipped if a better guess is already present in the state. The solve and solve_adjoint subroutines solve the primal (Step 3b) and adjoint linear systems (Step 3d). The energy subroutine computes and returns the energy (Step 3c), and, finally, the solvation_force_terms subroutine assembles and returns the solvent contributions to the forces (Step 3e). 3.5 | Step 4: Solute specific contributions to the derivatives The expressions that need to be evaluated in this step are quite general, as they depend on both the nature of the solute, and the specific differentiation variable λ. For this reason, in the most general case, the step should be carried out in the host code, for example, in case of the Fock matrix. However, ddX provides an implementation for this step for multipolar solutes and for λ being the nuclear coordi- nates. This is an interesting case, as the forces are often needed, they have more cumbersome expressions with respect to other analytical derivatives, and in ddX this step benefits from the internal FMM library. The solvent contribution to the forces is computed by calling multipole_force_terms in Fortran or the corresponding functions in C and Python. This subroutine takes as arguments the multipolar distribution, the model TABLE 2 Low-level API functions in Fortran, C and Python. Purpose Fortran C Python Solving the problem setup ddx_setup (State) fill_guess ddx_fill_guess State.fill_guess solve ddx_solve State.solve fill_guess_adjoint ddx_fill_guess_adjoint State.fill_guess_adjoint solve_adjoint ddx_solve_adjoint State.solve_adjoint energy ddx_energy State.energy solvation_force_terms ddx_solvation_force_terms State.solvation_force_terms NOTTOLI ET AL. 9 of 24 object and the state object. The latter contains information about the adjoint solution Xadj which is required to assemble the scalar products in Equation (10). 3.6 | Step 5: Finalization The finalization of the library is done by deallocating the memory which is occupied by the model object, the state object, and the electrostatics container. These deallocations are done in Fortran by, respectively, calling deallocate_model, deallocate_state and deallocate_electrostatics. In C, the library provides the corresponding functions, whereas in Python the deallocation is handled by the garbage collector. 3.7 | FMM acceleration Many of the steps reported in Section 3.1 require the computation of electrostatic interactions and are natively quadrati- cally scaling. Such steps are computing the RHS F, solving the linear systems for ddPCM and ddLPB, computing the solvation contribution to the forces for ddPCM and ddLPB, and finally computing the solute specific contribution to the forces. All of these steps are accelerated with the FMM79 and are linear scaling in N , except for the initialization, which scales like N log Nð Þ. The FMM engine implemented in ddX is tailored to accelerating all the calculations involved in the various domain decomposition algorithms. Here, we do not present our FMM implementation in detail, as it is already covered by sev- eral papers, and an in-depth explanation of the ddX internal FMM library is given in reference 67. However, we provide an explanation of the main differences with respect to traditional FMM schemes. • Clusterization: In ddX, we use a non-standard adaptive tree instead of the octree approaches used by standard librar- ies. The spheres building the cavity are usually arranged according to non-uniform distributions; hence, a standard uniform boxification scheme would produce an unbalanced octree, leading to the suboptimal efficiency of the FMM algorithm. This problem is solved using a different clusterization scheme, the recursive inertial bisection, at the cost of having non-standard (and thus non-precomputable) M2M and L2L operators from one tree level to another one. In the inertial bisection step, the distribution of points is split in two by finding the hyperplane that cuts along the moment of inertia. By applying this step recursively, we obtain the recursive inertial bisection algorithm.80,81 This guarantees a balanced tree and thus optimal number operations. The construction of the tree scales like N log Nð Þ as the standard octree generation. • Spherical harmonics: In usual libraries, the sources and targets are usually charges or at the lowest order Cartesian multipoles. In ddX, we allow the sources and targets to be arbitrary multipolar distributions in real spherical harmonics. • Different nature of sources and targets: Traditional libraries are often meant for symmetric interactions in which the nature of the sources and targets is similar. In ddX, the sources are the multipolar distributions at the center of the spheres, while the targets are the integration points on the surface of the spheres. This has two main conse- quences: (i) a particle-to-multipole step (P2M) is missing, whereas a local-to-particle (L2P) step takes care of evaluat- ing the electrostatic properties at the integration points; (ii) the whole FMM scheme is not symmetric, as there is an asymmetry in the targets and sources despite the kernel being symmetric. 3.8 | Software distribution and build system The source code of ddX is freely available on GitHub at https://github.com/ddsolvation/ddX under the LGPL v3 license with an automatically generated documentation available at https://ddsolvation.github.io/ddX/. In ddX we make use of the CMake build system to detect our dependencies and orchestrate the build of the library. By default, only the Fortran and C parts of the library are built and linked into a shared object libddx.so as well as the program ddx_driver. The former can be linked to third-party codes using the headers (ddx.h for C or compiler- specific .mod files for Fortran), while the latter provides a standalone runner for ddX solvation models based on a 10 of 24 NOTTOLI ET AL. https://github.com/ddsolvation/ddX https://ddsolvation.github.io/ddX/ plain text input file. In this form the dependencies of the code are minimal, being them limited to BLAS and LAPACK, as well as an OpenMP-compatible compiler. Alternatively, the library can be installed as a Python package using the pip utility and the Python setuptools module. In this form we provide the code through conda-forge and PiPy, allowing for an even more user-friendly installation—simply by issuing pip install pyddx or conda install conda-forge::pyddx. This also takes care of automatically providing the additional Python-specific dependencies (Pybind11, numpy, scipy) the pyddx Python package requires. 4 | INTERFACING In this section, we illustrate how an interface between ddX and an existing code can be achieved by illustrating in detail two examples: one with a QM code, namely, Psi4,69 and one with a classical molecular dynamics code, namely, Tin- ker.71 The first uses the Python bindings and the second uses the native Fortran API. 4.1 | Tinker-interface Tinker is a collection of different programs and tools that perform various kinds of calculations on systems described using classical molecular mechanics (MM).71 By coupling Tinker with ddX we want to allow Tinker to include solvation contributions to energy and forces computed with ddX. This requires four main modifications to Tinker: (i) as custom- ary in Tinker, implementing a global data structure that contains all the ddX intermediates and parameters; (ii) changing the input parsing to read also ddX parameters; (iii) changing the way the energy is computed; (iv) changing the way the forces are computed. Tinker mostly works using global variables, all the variables related to the same problem are collected in a For- tran module, and whenever they need to be used, the corresponding module is imported. In our case, we followed the same pattern and we created a ddx_interface module. This contains pointers to the ddX quantities (model, state, error, and electrostatics) and temporary storage for the parameters. We will not discuss the input parsing, as it is technical and unrelated to how to use ddX. It is sufficient to know that, when a Tinker input file is read, all the relevant ddX parameters are read and stored in the ddx_interface module. The initialization of ddX is not done at this step, as many important tasks of Tinker—such as molecular dynamics—require updating the geometry, and as the geometry changes the ddX initialization needs to be done from scratch. So it is more convenient to move the ddX initialization to where the energy and forces are computed. In Tinker, all the modules that deal with certain interaction terms are split into one subroutine that computes only the energy, and one that computes energy and forces. In our case, we slightly deviate from the standard pattern and we added a single subroutine called ddx_run which has a logical argument for doing the forces. In turn, this is called by the Tinker subroutines for solvation, if a ddX solvation is set. Algorithm 1 shows the inner workings of the ddx_run subroutine, the subroutines called reflect the example case reported in Figure 1. All the inputs are read from global variables in Fortran modules, and the outputs (energy and optionally forces) are also written to global variables in modules. The first three subroutines called in Lines 2, 3, and 4 are novel additions to Tinker and mostly take care of conver- ting Tinker quantities to ddX quantities, both in terms of format and units of measure. ddx_build_multipoles takes care of reading the multipoles from the Tinker arrays and converting them to real spherical harmonics in atomic units. ddx_build_cavity uses the atom types and positions to create a list of sphere centers and radii in atomic units. ddx_add_sph adds optional user-specified spheres to the list. Then, ddx_init is a simple wrapper to ddinit which takes also care of allocating Ψ. The calls at Lines 6 and 7 are direct to the ddX API and compute the primal and adjoint RHSs. The problem is solved by calling ddrun, which is part of the ddX API, and finally, if the forces are requested, the solute specific contribution is computed by calling multipole_force_terms which is again part of the ddX API. The finalization of the library is done by calling ddx_free, this is a wrapper to deallocate_model and deallocate_state which deallocates also Ψ and the arrays containing the multipoles, centers, and radii. To conclude, coupling Tinker with ddX was possible by adding a new module ddx_interface, an input parsing subroutine, three helper subroutines (ddx_build_multipoles, ddx_build_cavity, and ddx_add_sph), a NOTTOLI ET AL. 11 of 24 wrapper to ddinit, a finalizing subroutine (ddx_free) and a main driver (ddx_run) of only 166 lines of code including comments. 4.2 | Psi4-interface This section discusses the integration of ddX with Psi4,69 a software suite for performing a wide range of first-principle quantum-chemistry calculations. Currently the implicit solvation models of ddX can be used both to compute self- consistent ground states (e.g., HF or DFT methods) as well as subsequent linear response properties including excited states computations based on TDDFT. For more details on using ddX with Psi4 see the ddX section of the Psi4 user manual at https://psicode.org. Following the previous detailed discussion82 about how to couple domain-decomposition solvation methods to QM host codes, three main modifications were needed inside Psi4: (1) the handling of the additional input parameters for ddX, (2) the computation of solvation energy and Fock matrix contributions as well as (3) the addition of integration routines for electrostatic quantities, which are based on a Becke-type quadrature scheme, which Psi4 uses for DFT com- putations. This quadrature is used for ddX, since it is able to properly deal with the overlapping atom-centered balls Ωj of the ddX cavity.82 (1) and (3) are strongly dependent on the internals of Psi4, thus of minor interest to implement interfaces with other host codes; we will thus focus our discussion on point (2). 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 contri- bution 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, for example, in the first SCF iteration. This step processes the user-defined parame- ters 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 specific 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.61,63,66 Next, the nuclear contributions are added by calling model.multi- pole_electrostatics() and model.multipole_psi(). If this is the first iteration of the self-consistent field Algorithm 1 Using ddX from Tinker: Note that, for simplicity, most arguments are omitted from subroutine calls 12 of 24 NOTTOLI ET AL. https://psicode.org loop, the resulting F and Ψ are used to initialize the pyddx.State. This object in particular stores the RHS of the lin- ear systems (6) and (9) as well as iterated solution vectors X and Xadj. A following call to state.fill_guess() and state.fill_guess_adjoint() generates a crude initial guess for X and Xadj. 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 (6) 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 Xadj employed to compute the Fock matrix contributions—using again the host-specific inte- gration routines. 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 NOTTOLI ET AL. 13 of 24 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 ener- gies. For this, we follow the standard approach83 to only consider the fast (electronic) components in the solvent's response, which in Algorithm 2 reflects in two notable changes. First, we employ the optical permittivity ε∞ instead of the static permittivity εs when setting up the pyddx. Model and second, we skip over Lines 14 and 15, that is, we do not add the nuclear contribution to F and Ψ. 4.3 | Employing ddX-solvated ground states in adcc ADC-connect (adcc)70 is a software suite for performing excited states calculations based on the algebraic-diagrammatic construction scheme for the polarization propagator (ADC).84,85 A key focus of the code is a modular design, such that multiple host codes can be employed to supply the Hartree–Fock reference. Additional capabilities, such as interfaces to solvation models, can be passed from the host code to adcc as well. In combination with Psi4 this includes the treat- ment of solvation effects when computing excitation energies. In the past this strategy has already been employed to perform ADC calculations, which consider environmental effects via a polarizable embedding (PE)70,86 (using cppe87) or using the integral-equation formulation of PCM (using PCMSolver88). In a similar fashion, we extended adcc in this work to make use of ddX-based solvation models. Both a full linear response treatment89 as well as an approach based on perturbative linear-response corrections83 is now available in adcc in combination with ddX solvation models. In both cases Algorithm 2 is called directly from adcc to compute the effect of the fast solvent components on the excita- tion (i.e., we employ the optical permittivity and no nuclear contributions). This demonstrates the seamless integration of ddX within a larger existing Python-based ecosystem for quantum chemistry involving the third-party codes Psi4 and adcc. 4.4 | Standalone driver Alternatively, the library can also be run as a standalone using the previously mentioned ddx_driver. The program ddx_driver is a very simple Fortran code that performs two main functions: it performs the initialization from the file and then performs the steps given in Section 3.1. In this case, the initialization from the input file is done with a separate Fortran function called ddfromfile, which reads the text-based input file and then calls ddinit with appropri- ate parameters. 5 | EXAMPLE APPLICATIONS ddX was tested through the Tinker/ddX and the Psi4/ddX interfaces. All the calculations have been performed on the JUSTUS2 computer cluster,90 on nodes equipped with 2 Intel Xeon 6252 Gold and up to 384 GB of memory. For what concerns Tinker/ddX, we used our fork of Tinker 8 at branch ddx, commit d53669b, which is available on GitHub at ddsolvation/tinker-ddX. Tinker was then linked to the ddX shared library obtained by compiling version 0.6.0. Both Tinker and ddX were compiled using the Intel compiler version 19.1.2 and MKL libraries version 2020.2. For what con- cerns Psi4/ddX, we used the version of Psi4 at branch master, commit 990b0e8. The source code was built in a conda environment containing adcc version 0.15.17 and ddX version 0.6.0 provided through PyPI. Note that in this case, ddX was compiled using the GNU compiler available in the conda environment. Data for reproducibility is available at refer- ence 91. The archive contains all the input and output files of the simulations, the conda environment used for Psi4/ ddX and Psi4/adcc/ddX calculations, a copy of the employed versions of ddX, Tinker and Psi4, and the Python scripts used for generating the tables and plots of this paper. 5.1 | Benchmarks We investigated the scaling of the three models by running calculations on systems of increasing size. To do so, we first converted the protein data bank (PDB) files into Tinker input files using the Tinker tool pdb2xyz. We assigned the 14 of 24 NOTTOLI ET AL. parameters using the Amber 99 force field,92 which is bundled within Tinker. Further information about the input structures is given in Table 3. Then, for each input structure and for each model among vacuum, ddCOSMO, ddPCM, and ddLPB, we used the Tinker tool testgrad to compute the energy and the analytical forces. The structures were solvated using water ε¼ 80, and in case of ddLPB we used a Debye Hückel screening constant κs ¼ 0:104 A o �1. The execution time for the evaluation of the energy and analytical forces was measured using the time command, which also returns the maximum memory usage. The solvation energies computed with the three models are reported in Table 3. The time required by the various calculation steps is reported in Figure 2. The time required to assemble the adjoint RHS is negligible with respect to the rest; the next steps, in terms of cost, are the computation of the RHS and of the forces, which also show a perfectly linear scaling regime. Finally, the most expensive steps are the two linear systems; these also show a slightly super-linear scaling regime due to an increasing number of iterations for the largest systems. For the primal linear system of ddCOSMO, the analysis shows that the number of iterations depends on the globularity of the molecule; see reference 104. As shown in a few examples in reference 105, section 5, biological molecules like proteins and DNA often have a chain-like structure and are thus not globular, and a constant number of iterations can be expected, but with some natural variation as some molecules are more globular than others that explains this performance. TABLE 3 PDB codes and number of atoms for the structures used in the benchmark (upper part) and in the MD simulations (bottom part). PDB code Number of atoms Reference Es ddCOSMO Es ddPCM Es ddLPB 2p7r 70 93 �9.92 �9.74 �9.85 1etn 143 94 �91.1 �89.9 �91.5 1du9 381 95 �270 �267 �270 1caa 975 96 �960 �955 �966 1d3w 2051 97 �2651 �2641 �2673 1hde 9734 98 �4300 �4271 �4328 1ju2 20,288 99 �1989 �1870 �1978 6ftl 39,151 100 �10,057 �9655 �10,009 1stm 138,408 101 �80,489 �95,512 �80,916 1cnl 169 102 – – – 1ucs 997 103 – – – Note: The table reports also the solvation energies computed with the three different models, all the energies are in kcal/mol. FIGURE 2 Scaling of the various steps of a ddX calculation for systems of increasing size for the three models (ddCOSMO, ddPCM, ddLPB). The various steps are reported using different colors and different models using different panels. The gray dotted lines report an indication of the linear scaling regime. NOTTOLI ET AL. 15 of 24 The overall performance of ddX is remarkably good. For the largest system, which consists of more than 138,000 atoms, the total time required to compute the energy ranges between 3 min for ddCOSMO, 4 min for ddPCM and 56 min for ddLPB, the timings increasing to 9, 10 and 204 min, respectively, if also the forces are computed. We remark that the timings were obtained on a single computer node equipped with two Intel Xeon 6252 Gold CPUs, which is commonly available hardware, and that the results are obtained using conservative discretization parameters, that ensure no compromise on the accuracy of the compute gradients. Various previous studies can further illustrate the good performance: ddCOSMO was used to solvate viral capsides of up to 7 million atoms in less than 1 h,49 ddPCM was used to solvate a system of �625k atoms in �5 h on a single core,67 and finally ddLPB was tested against the TABI-PB and the APBS finite difference method LPB solvers, finding that it performs similarly for loose accuracies and that outperforms both when a higher accuracy is requested.68 Figure 3 reports the peak memory usage for the three models, which is only an estimate as the values are computed as the difference between the peak memory usage using ddX and the peak memory usage of a vacuum calculation. The three models present a perfectly linear scaling memory consumption, and as expected, the memory consumption increases with the complexity of the model, so it increases going from ddCOSMO to ddPCM and finally to ddLPB. 5.2 | Molecular dynamics The interface with Tinker and the rigorous analytical gradients implemented in ddX allow one to perform MD simula- tions. We tested the MD simulations on two systems (reported in the bottom of Table 3). For each PDB structure, the input files were prepared using the pdbxyz tool from Tinker and the Amber 99 sb force field,106 which is bundled in Tinker. The input files were used to run the first 2 ps of NVT MD simulation at a temperature of 300 K, and then the restart files were used to run 2 ps long NVE simulations. Also in this case, the structures were solvated using water ε¼ 80, with a Debye Hückel screening constant κs ¼ 0:104 A ∘ �1 in case of ddLPB. The simulations were repeated using two different threshold for the linear solvers of ddX, 10�4 as a loose threshold and 10�8 as a tight threshold. These two simple-minded examples show that the forces computed by ddX are indeed accurate, making thus rigorously energy conserving NVE simulations with ddX possible. Figure 4 reports an analysis of the energy conservation during the NVE simulations for the two systems and for the four models (including vacuum). The energy conservation is overall good, as it remains of the same order of magnitude as the one observed for the vacuum trajectory. In terms of relative energy conservation, the energy fluctuations in the antifreeze protein MD remain between �2% and þ4%, whereas the ones of the peptide remain in the �2% range. The usage of a tighter threshold mostly affects the short-time fluctuations, with a more pronounced effect on ddLPB due to the inner–outer linear system solver. 5.3 | Psi4 calculations The ddX interface with Psi4 was tested by comparing between using ddX and using PCMSolver,88 which we use as a reference. As a test molecule, we choose the dye nile red. This molecule presents a large solvatochromic red-shift and is of medium size (with 42 atoms) which makes it an interesting test case for continuum solvation models. FIGURE 3 Memory consumption of ddX in an energy and forces calculation done with the Tinker/ddX implementation. The three models (ddCOSMO, ddPCM, and ddLPB) are reported with different colors. The gray dotted lines report an indication of the linear scaling regime. 16 of 24 NOTTOLI ET AL. First, we optimized the ground state geometry of nile red using the MP2 method with a cc-pVDZ basis set. The opti- mized ground state geometry was then used to run both MP2 calculations and TDDFT calculations. In this case, contin- uum solvation was enabled and the calculations were repeated in three different solvents, water (wat, εs ¼ 80:1, ε∞ ¼ 1:777), acetonitrile (acn, εs ¼ 37:5, ε∞ ¼ 1:807), and cyclohexane (cyhex, εs ¼ 2:02, ε∞ ¼ 2:034). The MP2 results are reported in Table 4, whereas the TDDFT results are reported in Table 5. The ddX results agree perfectly with the PCMSolver results, confirming the validity of the implementation. It was not possible to perform the PCMSolver calculations in cyclohexane due to numerical errors, and in this case, ddX showed its better stability. We did not report the timings, as all the investigated models performed similarly, suggesting that the bottleneck of the cal- culations is not solving the linear systems but rather computing the electronic integrals at the grid points, an operation which is similar for both ddX and PCMSolver. However, it should be noted that since ddX is entirely linear scaling, we can expect a molecule size at which the PCMSolver calculation becomes more expensive as the quadratic contribution for solving the linear system becomes dominant. FIGURE 4 Energy conservation during the NVE trajectories, the values reported are the total energy minus the starting energy of the NVE simulation. The different models (vacuum, ddCOSMO, ddPCM, and ddLPB) are reported using different colors, the combination of the two systems, and the two convergence thresholds (loose 10�4, tight 10�8) using different panels. The dark curves report a moving average of the energies over 100 steps to improve the readability, and the light curves report the energies and show the short-term fluctuations. TABLE 4 Results of MP2 calculations on nile red, done in three solvents (water, acetonitrile and cyclohexane), with methods from PCMSolver (PCM, COSMO) and from ddX (ddPCM, ddCOSMO). Solvent Model Energy (Hartree) wat ddPCM �1026.5343 wat PCM �1026.5343 wat ddCOSMO �1026.5343 wat COSMO �1026.5344 acn ddPCM �1026.5340 acn PCM �1026.5341 acn ddCOSMO �1026.5342 acn COSMO �1026.5342 cyhex ddPCM �1026.5279 cyhex PCM – cyhex ddCOSMO �1026.5290 cyhex COSMO �1026.5290 vac – �1026.5241 Note: The results of calculation done in vacuum are added as a reference. NOTTOLI ET AL. 17 of 24 To conclude, the ddX implementation in Psi4 shows similar performance and improved stability with respect to tra- ditional implementations while generating the same results. Further, for this size of the solute, there is no significant difference in the timings between ddCOSMO and ddPCM as opposed to PCMsolver. Overall, the Psi4-ddX interface is ready to be used in production calculations and its linear scaling allows us to consider much larger molecules than the nile red molecule. 5.4 | ADC(2) energies The adcc support for ddX was tested by performing ADC(2) calculations on the p-nitroaniline molecule. This molecule is a push–pull system characterized by a bright absorption due to a charge transfer transition. The band is subject to a large solvatochromic red-shift when going from low to high polarity.107 The initial structures were obtained by performing MP2 geometry optimizations with the cc-pVDZ basis set. The ini- tial structure was first optimized in vacuum, and then the resulting structure was used as a starting point for three geometry optimizations in the investigated solvents. The solvation effects were included using ddPCM. As the analytical forces are not yet supported in Psi4/ddX, we used the automatic numerical differentiation provided by Psi4 to carry out the optimizations. On the final four structures, we performed excited state calculations using ADC(2). The solvation effect on the excited states was included through a perturbative linear response correction.86 The results are reported in Table 6. The ground state energy (Hartree–Fock) progressively decreases with an increas- ing dielectric constant of the environment. The first and third excited states are dark and are only slightly influenced by the different dielectric constants of the environment, whereas the second excited state is bright and undergoes the TABLE 5 Results of TDDFT calculations on nile red, done in three solvents (water, acetonitrile and cyclohexane), with methods from PCMSolver (PCM, COSMO) and from ddX (ddPCM, ddCOSMO). Solvent Model Energy (Hartree) Exc (eV) wat ddPCM �1032.4381 2.80 3.56 3.63 4.11 4.33 wat PCM �1032.4382 2.79 3.56 3.63 4.11 4.33 wat ddCOSMO �1032.4382 2.78 3.56 3.63 4.11 4.33 wat COSMO �1032.4382 2.78 3.56 3.63 4.11 4.33 acn ddPCM �1032.4379 2.80 3.56 3.64 4.11 4.33 acn PCM �1032.4380 2.79 3.56 3.63 4.11 4.33 acn ddCOSMO �1032.4380 2.78 3.56 3.63 4.11 4.33 acn COSMO �1032.4381 2.78 3.56 3.63 4.11 4.33 cyhex ddPCM �1032.4326 2.85 3.45 3.68 4.12 4.31 cyhex PCM – – – – – – cyhex ddCOSMO �1032.4334 2.82 3.46 3.68 4.12 4.30 cyhex COSMO �1032.4334 2.82 3.46 3.68 4.12 4.30 vac – �1032.4292 2.99 3.37 3.72 4.12 4.30 Note: The results of calculation in vacuum are added as a reference. TABLE 6 Results of ADC(2) calculations on p-nitroaniline, done in three solvents (water, acetonitrile, and cyclohexane) and in vacuum. Solvent Energy (Hartree) Exc (eV) wat �489.2550 3.80 3.97 4.54 acn �489.2547 3.80 3.98 4.54 cyhex �489.2483 3.81 4.26 4.69 vac �489.2436 3.76 4.53 4.72 18 of 24 NOTTOLI ET AL. expected red-shift when going from a low polarity solvent to a high polarity solvent. The general trend is in agreement with experimental results108–110; however, a quantitative agreement is missing mostly due to limitations of the contin- uum solvation models.107 6 | CONCLUSION AND OUTLOOK The ddX library represents the culmination of over 10 years of dedicated effort in the field of polarizable continuum sol- vation models (PCSM). With the library, we provide three PCSMs based on the domain decomposition paradigm, namely ddCOSMO, ddPCM, and ddLPB, under a single easy-to-use API available via multiple programming languages. Notably, the implemented methods are robust and exhibit a linear scaling regime in both time and memory. Being interfaced to multiple host codes ddX is already able to tackle solvation effects in molecular systems ranging from mole- cules to proteins, thus bridging for the first time the application regimes of quantum chemistry and biophysics. Based on our standardized API we expect interfacing with further quantum-chemistry and molecular mechanics software packages to be easily feasible. We believe ddX to be a versatile ingredient in the toolbox of modular open- source atomistic modeling packages. In the future, we aspire to include GPU acceleration and more involved para- llelization strategies in ddX, developments which will directly be available to all interfaced packages. Also, we would like to implement the SES-based domain-decomposition method74 to support SES-based cavities. AUTHOR CONTRIBUTIONS Michele Nottoli: Investigation (equal); methodology (equal); software (equal); writing – original draft (equal); writing – review and editing (equal). Michael F. Herbst: Investigation (equal); methodology (equal); software (equal); writing – original draft (equal); writing – review and editing (equal). Aleksandr Mikhalev: Investigation (equal); methodology (equal); software (equal). Abhinav Jha: Investigation (equal); methodology (equal); software (equal); writing – original draft (equal); writing – review and editing (equal). Filippo Lipparini: Investigation (equal); methodology (equal); software (equal); supervision (equal); writing – original draft (equal); writing – review and editing (equal). Benjamin Stamm: Funding acquisition (equal); investigation (equal); methodology (equal); software (equal); supervision (equal); writing – original draft (equal); writing – review and editing (equal). FUNDING INFORMATION AJ and BS acknowledge support from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under project 440641818. We thank the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) for supporting this work by funding—EXC2075-390740016 under Germany's Excellence Strategy. We acknowledge the sup- port by the Stuttgart Center for Simulation Science (SimTech). The authors acknowledge support by the state of Baden- Württemberg through bwHPC and the German Research Foundation (DFG) through grant no INST 40/575-1 FUGG (JUSTUS 2 cluster). CONFLICT OF INTEREST STATEMENT The authors declare no conflicts of interest. DATA AVAILABILITY STATEMENT The data required to reproduce all the calculations presented in this software focus are available at https://doi.org/10. 18419/darus-4030. ORCID Michele Nottoli https://orcid.org/0000-0002-6544-0897 Michael F. Herbst https://orcid.org/0000-0003-0378-7921 Filippo Lipparini https://orcid.org/0000-0002-4947-3912 Benjamin Stamm https://orcid.org/0000-0003-3375-483X RELATED WIREs ARTICLES Dielectric continuum methods for quantum chemistry NOTTOLI ET AL. 19 of 24 https://doi.org/10.18419/darus%E2%80%904030 https://doi.org/10.18419/darus%E2%80%904030 https://orcid.org/0000-0002-6544-0897 https://orcid.org/0000-0002-6544-0897 https://orcid.org/0000-0003-0378-7921 https://orcid.org/0000-0003-0378-7921 https://orcid.org/0000-0002-4947-3912 https://orcid.org/0000-0002-4947-3912 https://orcid.org/0000-0003-3375-483X https://orcid.org/0000-0003-3375-483X https://doi.org/10.1002/wcms.1519 REFERENCES 1. Cramer CJ, Truhlar DG. Implicit solvation models: equilibria, structure, spectra, and dynamics. Chem Rev. 1999;99(8):2161–200. https://doi.org/10.1021/cr960149m 2. Tomasi J, Mennucci B, Cammi R. Quantum mechanical continuum solvation models. Chem Rev. 2005;105(8):2999–3094. https://doi. org/10.1021/cr9904009 3. Lange AW, Herbert JM. Polarizable continuum reaction-field solvation models affording smooth potential energy surfaces. J Phys Chem Lett. 2009;1(2):556–61. https://doi.org/10.1021/jz900282c 4. Mennucci B. Polarizable continuum model. WIREs Comput Mol Sci. 2012;2(3):386–404. https://doi.org/10.1002/wcms.1086 5. Klamt A. The COSMO and COSMO-RS solvation models. WIREs Comput Mol Sci. 2018;8:e1338. https://doi.org/10.1002/wcms.1338 6. Herbert JM. Dielectric continuum methods for quantum chemistry. WIREs Comput Mol Sci. 2021;11(4):e1519. https://doi.org/10.1002/ wcms.1519 7. Roux B, Simonson T. Implicit solvent models. Biophys Chem. 1999;78(1–2):1–20. https://doi.org/10.1016/s0301-4622(98)00226-9 8. Demerdash O, Yap E-H, Head-Gordon T. Advanced potential energy surfaces for condensed phase simulation. Annu Rev Phys Chem. 2014;65(1):149–74. https://doi.org/10.1146/annurev-physchem-040412-110040 9. Decherchi S, Masetti M, Vyalov I, Rocchia W. Implicit solvent methods for free energy estimation. Eur J Med Chem. 2015;91:27–42. https://doi.org/10.1016/j.ejmech.2014.08.064 10. Mennucci B, Cances E, Tomasi J. Evaluation of solvent effects in isotropic and anisotropic dielectrics and in ionic solutions with a uni- fied integral equation method: theoretical bases, computational implementation, and numerical applications. J Phys Chem B. 1997; 101(49):10506–17. https://doi.org/10.1021/jp971959k 11. Tomasi J, Mennucci B, Cances E. The IEF version of the PCM solvation method: an overview of a new method addressed to study molecular solutes at the QM ab initio level. J Mol Struct (THEOCHEM). 1999;464(1–3):211–26. https://doi.org/10.1016/s0166-1280(98) 00553-3 12. Cances E, Mennucci B, Tomasi J. A new integral equation formalism for the polarizable continuum model: theoretical background and applications to isotropic and anisotropic dielectrics. J Chem Phys. 1997;107(8):3032–41. https://doi.org/10.1063/1.474659 13. Klamt A. The COSMO and COSMO-RS solvation models. WIREs Comput Mol Sci. 2011;1(5):699–709. https://doi.org/10.1002/wcms.56 14. Klamt A, Schuurmann G. COSMO: a new approach to dielectric screening in solvents with explicit expressions for the screening energy and its gradient. J Chem Soc Perkin Trans. 1993;2(5):799–805. https://doi.org/10.1039/p29930000799 15. Truong TN, Stefanovich EV. A new method for incorporating solvent effect into the classical, ab initio molecular orbital and density functional theory frameworks for arbitrary shape cavity. Chem Phys Lett. 1995;240(4):253–60. https://doi.org/10.1016/0009-2614(95) 00541-b 16. Barone V, Cossi M. Quantum calculation of molecular energies and energy gradients in solution by a conductor solvent model. Chem A Eur J. 1998;102(11):1995–2001. https://doi.org/10.1021/jp9716997 17. Gregersen BA, Khandogin J, Thiel W, York DM. Smooth solvation method for d-orbital semiempirical calculations of biological reac- tions. 2. Application to transphosphorylation thio effects in solution. J Phys Chem B. 2005;109(19):9810–7. https://doi.org/10.1021/ jp044061l 18. Benighaus T, Thiel W. Efficiency and accuracy of the generalized solvent boundary potential for hybrid QM/MM simulations: imple- mentation for semiempirical Hamiltonians. J Chem Theory Comput. 2008;4(10):1600–9. https://doi.org/10.1021/ct800193a 19. Benighaus T, Thiel W. A general boundary potential for hybrid QM/MM simulations of solvated biomolecular systems. J Chem Theory Comput. 2009;5(11):3114–28. https://doi.org/10.1021/ct900437b 20. Kříž K, Řez�ač J. Reparametrization of the COSMO solvent model for semiempirical methods PM6 and PM7. J Chem Inf Model. 2019; 59(1):229–35. https://doi.org/10.1021/acs.jcim.8b00681 21. Aguilar MA, Olivares del Valle FJ, Tomasi J. Electron correlation and solvation effects. II. The description of the vibrational properties of a water molecule in a dielectric given by the application of the polarizable continuum model with inclusion of correlation effects. Chem Phys. 1991;150(2):151–61. https://doi.org/10.1016/0301-0104(91)80125-2 22. Olivares del Valle FJ, Tomasi J. Electron correlation and solvation effects. I. Basic formulation and preliminary attempt to include the electron correlation in the quantum mechanical polarizable continuum model so as to study solvation phenomena. Chem Phys. 1991; 150(2):139–50. https://doi.org/10.1016/0301–0104(91)80124-z 23. Lipparini F, Scalmani G, Mennucci B. Non covalent interactions in RNA and DNA base pairs: a quantum-mechanical study of the cou- pling between solvent and electronic density. Phys Chem Chem Phys. 2009;11(48):11617. https://doi.org/10.1039/b915898g 24. Cammi R. Quantum cluster theory for the polarizable continuum model. I. The CCSD level with analytical first and second derivatives. J Chem Phys. 2009;131(16):164104. https://doi.org/10.1063/1.3245400 25. Caricato M. Coupled cluster theory with the polarizable continuum model of solvation. Int J Quantum Chem. 2018;119(1):e25710. https://doi.org/10.1002/qua.25710 26. Scalmani G, Barone V, Kudin KN, Pomelli CS, Scuseria GE, Frisch MJ. Achieving linear-scaling computational cost for the polarizable continuum model of solvation. Theor Chem Acc. 2004;111(2–6):90–100. https://doi.org/10.1007/s00214-003-0527-2 27. Gouy M. Sur la constitution de la charge electrique a la surface d'un electrolyte. J Phys Theor Appl. 1910;9(1):457–68. https://doi.org/ 10.1051/jphystap:019100090045700 28. Chapman DL. LI. A contribution to the theory of electrocapillarity. Lond Edinb Dublin Philos Mag J Sci. 1913;25(148):475–81. https:// doi.org/10.1080/14786440408634187 20 of 24 NOTTOLI ET AL. https://doi.org/10.1021/cr960149m https://doi.org/10.1021/cr9904009 https://doi.org/10.1021/cr9904009 https://doi.org/10.1021/jz900282c https://doi.org/10.1002/wcms.1086 https://doi.org/10.1002/wcms.1338 https://doi.org/10.1002/wcms.1519 https://doi.org/10.1002/wcms.1519 https://doi.org/10.1016/s0301-4622(98)00226-9 https://doi.org/10.1146/annurev-physchem-040412-110040 https://doi.org/10.1016/j.ejmech.2014.08.064 https://doi.org/10.1021/jp971959k https://doi.org/10.1016/s0166-1280(98)00553-3 https://doi.org/10.1016/s0166-1280(98)00553-3 https://doi.org/10.1063/1.474659 https://doi.org/10.1002/wcms.56 https://doi.org/10.1039/p29930000799 https://doi.org/10.1016/0009-2614(95)00541-b https://doi.org/10.1016/0009-2614(95)00541-b https://doi.org/10.1021/jp9716997 https://doi.org/10.1021/jp044061l https://doi.org/10.1021/jp044061l https://doi.org/10.1021/ct800193a https://doi.org/10.1021/ct900437b https://doi.org/10.1021/acs.jcim.8b00681 https://doi.org/10.1016/0301-0104(91)80125-2 https://doi.org/10.1016/0301%E2%80%930104(91)80124-z https://doi.org/10.1039/b915898g https://doi.org/10.1063/1.3245400 https://doi.org/10.1002/qua.25710 https://doi.org/10.1007/s00214-003-0527-2 https://doi.org/10.1051/jphystap:019100090045700 https://doi.org/10.1051/jphystap:019100090045700 https://doi.org/10.1080/14786440408634187 https://doi.org/10.1080/14786440408634187 29. Debye P, Huckel E. Zur Theorie der Elektrolyte. I Gefrierpunktserniedrigung und verwandte Erscheinungen. Phys Z. 1923;24(185):305. 30. Gilson MK, Rashin A, Fine R, Honig B. On the calculation of electrostatic interactions in proteins. J Mol Biol. 1985;184(3):503–16. https://doi.org/10.1016/0022-2836(85)90297-9 31. Baker NA, Sept D, Joseph S, Holst MJ, McCammon JA. Electrostatics of nanosystems: application to microtubules and the ribosome. Proc Natl Acad Sci. 2001;98(18):10037–41. https://doi.org/10.1073/pnas.181342398 32. Jurrus E, Engel D, Star K, Monson K, Brandi J, Felberg LE, et al. Improvements to the APBS biomolecular solvation software suite. Pro- tein Sci. 2017;27(1):112–28. https://doi.org/10.1002/pro.3280 33. Li C, Jia Z, Chakravorty A, Pahari S, Peng Y, Basu S, et al. DelPhi suite: new developments and review of functionalities. J Comput Chem. 2019;40(28):2502–8. https://doi.org/10.1002/jcc.26006 34. Bond SD, Chaudhry JH, Cyr EC, Olson LN. A first-order system least-squares finite element method for the Poisson–Boltzmann equa- tion. J Comput Chem. 2009;31(8):1625–35. https://doi.org/10.1002/jcc.21446 35. Holst M, McCammon JA, Yu Z, Zhou Y, Zhu Y. Adaptive finite element modeling techniques for the Poisson–Boltzmann equation. Commun Comput Phys. 2012;11(1):179–214. https://doi.org/10.4208/cicp.081009.130611a 36. Xie D. New solution decomposition and minimization schemes for Poisson–Boltzmann equation in calculation of biomolecular electro- statics. J Comput Phys. 2014;275:294–309. https://doi.org/10.1016/j.jcp.2014.07.012 37. Lu B, Cheng X, Huang J, McCammon JA. Order N algorithm for computation of electrostatic interactions in biomolecular systems. Proc Natl Acad Sci. 2006;103(51):19314–9. https://doi.org/10.1073/pnas.0605166103 38. Lotan I, Head-Gordon T. An analytical electrostatic model for salt screened interactions between multiple proteins. J Chem Theory Comput. 2006;2(3):541–55. https://doi.org/10.1021/ct050263p 39. Yap E-H, Head-Gordon T. Calculating the bimolecular rate of protein–protein association with interacting crowders. J Chem Theory Comput. 2013;9(5):2481–9. https://doi.org/10.1021/ct400048q 40. Bajaj C, Chen S-C, Rand A. An efficient higher-order fast multipole boundary element solution for Poisson–Boltzmann-based molecu- lar electrostatics. SIAM J Sci Comput. 2011;33(2):826–48. https://doi.org/10.1137/090764645 41. Cooper CD, Bardhan JP, Barba LA. A biomolecular electrostatics solver using Python, GPUs and boundary elements that can handle solvent-filled cavities and Stern layers. Comput Phys Commun. 2014;185(3):720–9. https://doi.org/10.1016/j.cpc.2013.10.028 42. Cooper CD. A boundary-integral approach for the Poisson–Boltzmann equation with polarizable force fields. J Comput Chem. 2019; 40(18):1680–92. https://doi.org/10.1002/jcc.25820 43. Search SD, Cooper CD, Van't Wout E. Towards optimal boundary integral formulations of the Poisson–Boltzmann equation for molec- ular electrostatics. J Comput Chem. 2022;43(10):674–91. https://doi.org/10.1002/jcc.26825 44. Addison-Smith I, Guzman HV, Cooper CD. Accurate boundary integral formulations for the calculation of electrostatic forces with an implicit-solvent model. J Chem Theory Comput. 2023;19(10):2996–3006. https://doi.org/10.1021/acs.jctc.3c00021 45. Geng W, Krasny R. A treecode-accelerated boundary integral Poisson–Boltzmann solver for electrostatics of solvated biomolecules. J Comput Phys. 2013;247:62–78. https://doi.org/10.1016/j.jcp.2013.03.056 46. Krasny R, Wang L. A treecode based on barycentric Hermite interpolation for electrostatic particle interactions. Comput Math Biophys. 2019;7(1):73–84. https://doi.org/10.1515/cmb-2019-0006 47. Wilson L, Geng W, Krasny R. TABI-PB 2.0: an improved version of the Treecode-accelerated boundary integral Poisson–Boltzmann solver. J Phys Chem B. 2022;126(37):7104–13. https://doi.org/10.1021/acs.jpcb.2c04604 48. Rega N, Brancato G, Barone V. Non-periodic boundary conditions for ab initio molecular dynamics in condensed phase using localized basis functions. Chem Phys Lett. 2006;422(4–6):367–71. https://doi.org/10.1016/j.cplett.2006.02.051 49. Nottoli M, Nifosì R, Mennucci B, Lipparini F. Energy, structures, and response properties with a fully coupled QM/- AMOEBA/ddCOSMO implementation. J Chem Theory Comput. 2021;17(9):5661–72. https://doi.org/10.1021/acs.jctc.1c00555 50. Brancato G, Barone V, Rega N. Theoretical modeling of spectroscopic properties of molecules in solution: toward an effective dynami- cal discrete/continuum approach. Theor Chem Acc. 2007;117(5–6):1001–15. https://doi.org/10.1007/s00214-006-0216-z 51. Brancato G, Rega N, Barone V. A hybrid explicit/implicit solvation method for first-principle molecular dynamics simulations. J Chem Phys. 2008;128(14):144501. https://doi.org/10.1063/1.2897759 52. Mancini G, Fusè M, Lipparini F, Nottoli M, Scalmani G, Barone V. Molecular dynamics simulations enforcing nonperiodic boundary conditions: new developments and application to the solvent shifts of Nitroxide magnetic parameters. J Chem Theory Comput. 2022; 18(4):2479–93. https://doi.org/10.1021/acs.jctc.2c00046 53. Steindal AH, Ruud K, Frediani L, Aidas K, Kongsted J. Excitation energies in solution: the fully polarizable QM/MM/PCM method. J Phys Chem B. 2011;115(12):3027–37. https://doi.org/10.1021/jp1101913 54. Caprasecca S, Curutchet C, Mennucci B. Toward a unified modeling of environment and bridge-mediated contributions to electronic energy transfer: a fully polarizable QM/MM/PCM approach. J Chem Theory Comput. 2012;8(11):4462–73. https://doi.org/10.1021/ ct300620w 55. Caprasecca S, Jurinovich S, Lagardère L, Stamm B, Lipparini F. Achieving linear scaling in computational cost for a fully polarizable MM/continuum embedding. J Chem Theory Comput. 2015;11(2):694–704. https://doi.org/10.1021/ct501087m 56. Lipparini F, Lagardère L, Scalmani G, Stamm B, Cancès E, Maday Y, et al. Quantum calculations in solution for large to very large mol- ecules: a new linear scaling QM/continuum approach. J Phys Chem Lett. 2014;5(6):953–8. https://doi.org/10.1021/jz5002506 NOTTOLI ET AL. 21 of 24 https://doi.org/10.1016/0022-2836(85)90297-9 https://doi.org/10.1073/pnas.181342398 https://doi.org/10.1002/pro.3280 https://doi.org/10.1002/jcc.26006 https://doi.org/10.1002/jcc.21446 https://doi.org/10.4208/cicp.081009.130611a https://doi.org/10.1016/j.jcp.2014.07.012 https://doi.org/10.1073/pnas.0605166103 https://doi.org/10.1021/ct050263p https://doi.org/10.1021/ct400048q https://doi.org/10.1137/090764645 https://doi.org/10.1016/j.cpc.2013.10.028 https://doi.org/10.1002/jcc.25820 https://doi.org/10.1002/jcc.26825 https://doi.org/10.1021/acs.jctc.3c00021 https://doi.org/10.1016/j.jcp.2013.03.056 https://doi.org/10.1515/cmb-2019-0006 https://doi.org/10.1021/acs.jpcb.2c04604 https://doi.org/10.1016/j.cplett.2006.02.051 https://doi.org/10.1021/acs.jctc.1c00555 https://doi.org/10.1007/s00214-006-0216-z https://doi.org/10.1063/1.2897759 https://doi.org/10.1021/acs.jctc.2c00046 https://doi.org/10.1021/jp1101913 https://doi.org/10.1021/ct300620w https://doi.org/10.1021/ct300620w https://doi.org/10.1021/ct501087m https://doi.org/10.1021/jz5002506 57. Egidi F, Carnimeo I, Cappelli C. Optical rotatory dispersion of methyloxirane in aqueous solution: assessing the performance of density functional theory in combination with a fully polarizable QM/MM/PCM approach. Opt Mater Express. 2014;5(1):196. https://doi.org/ 10.1364/ome.5.000196 58. Lipparini F, Lagardère L, Raynaud C, Stamm B, Cancès E, Mennucci B, et al. Polarizable molecular dynamics in a polarizable contin- uum solvent. J Chem Theory Comput. 2015;11(2):623–34. https://doi.org/10.1021/ct500998q 59. Schnieders MJ, Baker NA, Ren P, Ponder JW. Polarizable atomic multipole solutes in a Poisson–Boltzmann continuum. J Chem Phys. 2007;126(12):124114. https://doi.org/10.1063/1.2714528 60. Cances E, Maday Y, Stamm B. Domain decomposition for implicit solvation models. J Chem Phys. 2013;139(5):054111. https://doi.org/ 10.1063/1.4816767 61. Lipparini F, Stamm B, Cancès E, Maday Y, Mennucci B. Fast domain decomposition algorithm for continuum solvation models: energy and first derivatives. J Chem Theory Comput. 2013;9(8):3637–48. https://doi.org/10.1021/ct400280b 62. Nottoli M, Mikhalev A, Stamm B, Lipparini F. Coarse-graining ddCOSMO through an interface between tinker and the ddX library. J Phys Chem B. 2022;126(43):8827–37. https://doi.org/10.1021/acs.jpcb.2c04579 63. Stamm B, Cancès E, Lipparini F, Maday Y. A new discretization for the polarizable continuum model within the domain decomposi- tion paradigm. J Chem Phys. 2016;144(5):054101. https://doi.org/10.1063/1.4940136 64. Gatto P, Lipparini F, Stamm B. Computation of forces arising from the polarizable continuum model within the domain-decomposition paradigm. J Chem Phys. 2017;147(22):224108. https://doi.org/10.1063/1.5008329 65. Nottoli M, Stamm B, Scalmani G, Lipparini F. Quantum calculations in solution of energies, structures, and properties with a domain decomposition polarizable continuum model. J Chem Theory Comput. 2019;15(11):6061–73. https://doi.org/10.1021/acs.jctc.9b00640 66. Quan C, Stamm B, Maday Y. A domain decomposition method for the Poisson–Boltzmann solvation models. SIAM J Sci Comput. 2019;41(2):B320–50. https://doi.org/10.1137/18m119553x 67. Mikhalev A, Nottoli M, Stamm B. Linearly scaling computation of ddPCM solvation energy and forces using the fast multipole method. J Chem Phys. 2022;157(11):114103. https://doi.org/10.1063/5.0104536 68. Jha A, Nottoli M, Mikhalev A, Quan C, Stamm B. Linear scaling computation of forces for the domain-decomposition linear Poisson– Boltzmann method. J Chem Phys. 2023;158(10):104105. https://doi.org/10.1063/5.0141025 69. Smith DGA, Burns LA, Simmonett AC, Parrish RM, Schieber MC, Galvelis R, et al. PSI4 1.4: open-source software for high-throughput quantum chemistry. J Chem Phys. 2020;152(18):184108. https://doi.org/10.1063/5.0006002 70. Herbst MF, Scheurer M, Fransson T, Rehn DR, Dreuw A. adcc: a versatile toolkit for rapid development of algebraic-diagrammatic con- struction methods. WIREs Comput Mol Sci. 2020;10(6):e1462. https://doi.org/10.1002/wcms.1462 71. Rackers JA, Wang Z, Lu C, Laury ML, Lagardère L, Schnieders MJ, et al. Tinker 8: software tools for molecular design. J Chem Theory Comput. 2018;14(10):5273–89. https://doi.org/10.1021/acs.jctc.8b00529 72. Lee B, Richards FM. The interpretation of protein structures: estimation of static accessibility. J Mol Biol. 1971;55(3):379-IN4. https:// doi.org/10.1016/0022-2836(71)90324-x 73. Richards FM. Areas, volumes, packing, and protein structure. Annu Rev Biophys Bioeng. 1977;6(1):151–76. https://doi.org/10.1146/ annurev.bb.06.060177.001055 74. Quan C, Stamm B, Maday Y. A domain decomposition method for the polarizable continuum model based on the solvent excluded sur- face. Math Models Methods Appl Sci. 2018;28(7):1233–66. https://doi.org/10.1142/s0218202518500331 75. Cances E, Mennucci B. The escaped charge problem in solvation continuum models. J Chem Phys. 2001;115(13):6130–5. https://doi. org/10.1063/1.1401157 76. Chipman DM. Simulation of volume polarization in reaction field theory. J Chem Phys. 1999;110(16):8012–8. https://doi.org/10.1063/1. 478729 77. Chipman DM. New formulation and implementation for volume polarization in dielectric continuum theory. J Chem Phys. 2006; 124(22):224111. https://doi.org/10.1063/1.2203068 78. Stern O. Zur theorie der elektrolytischen doppelschicht. Z Elektrochem Angew Phys Chem. 1924;30(21–22):508–16. https://doi.org/10. 1002/bbpc.192400182 79. Greengard L, Rokhlin V. A fast algorithm for particle simulations. J Comput Phys. 1987;73(2):325–48. https://doi.org/10.1016/0021- 9991(87)90140-9 80. Pothen A. Graph partitioning algorithms with applications to scientific computing. Parallel numerical algorithms. Dordrecht: Springer Netherlands; 1997. p. 323–68. https://doi.org/10.1007/978-94-011-5412-3_12 81. Leland R, Hendrickson B. An empirical study of static load balancing algorithms. Proceedings of IEEE scalable high performance com- puting conference. SHPCC-94. Knoxville, TN: IEEE Computer Society Press; 1994. https://doi.org/10.1109/shpcc.1994.296707 82. Stamm B, Lagardère L, Scalmani G, Gatto P, Cancès E, Piquemal J-P, et al. How to make continuum solvation incredibly fast in a few simple steps: a practical guide to the domain decomposition paradigm for the conductor-like screening model. Int J Quantum Chem. 2018;119(1):e25669. https://doi.org/10.1002/qua.25669 83. Cammi R, Corni S, Mennucci B, Tomasi J. Electronic excitation energies of molecules in solution: state specific and linear response methods for nonequilibrium continuum solvation models. J Chem Phys. 2005;122(10):104513. https://doi.org/10.1063/1.1867373 84. Schirmer J. Beyond the random-phase approximation: a new approximation scheme for the polarization propagator. Phys Rev A. 1982; 26(5):2395–416. https://doi.org/10.1103/physreva.26.2395 22 of 24 NOTTOLI ET AL. https://doi.org/10.1364/ome.5.000196 https://doi.org/10.1364/ome.5.000196 https://doi.org/10.1021/ct500998q https://doi.org/10.1063/1.2714528 https://doi.org/10.1063/1.4816767 https://doi.org/10.1063/1.4816767 https://doi.org/10.1021/ct400280b https://doi.org/10.1021/acs.jpcb.2c04579 https://doi.org/10.1063/1.4940136 https://doi.org/10.1063/1.5008329 https://doi.org/10.1021/acs.jctc.9b00640 https://doi.org/10.1137/18m119553x https://doi.org/10.1063/5.0104536 https://doi.org/10.1063/5.0141025 https://doi.org/10.1063/5.0006002 https://doi.org/10.1002/wcms.1462 https://doi.org/10.1021/acs.jctc.8b00529 https://doi.org/10.1016/0022-2836(71)90324-x https://doi.org/10.1016/0022-2836(71)90324-x https://doi.org/10.1146/annurev.bb.06.060177.001055 https://doi.org/10.1146/annurev.bb.06.060177.001055 https://doi.org/10.1142/s0218202518500331 https://doi.org/10.1063/1.1401157 https://doi.org/10.1063/1.1401157 https://doi.org/10.1063/1.478729 https://doi.org/10.1063/1.478729 https://doi.org/10.1063/1.2203068 https://doi.org/10.1002/bbpc.192400182 https://doi.org/10.1002/bbpc.192400182 https://doi.org/10.1016/0021-9991(87)90140-9 https://doi.org/10.1016/0021-9991(87)90140-9 https://doi.org/10.1007/978-94-011-5412-3_12 https://doi.org/10.1109/shpcc.1994.296707 https://doi.org/10.1002/qua.25669 https://doi.org/10.1063/1.1867373 https://doi.org/10.1103/physreva.26.2395 85. Dreuw A, Papapostolou A, Dempwolff AL. Algebraic diagrammatic construction schemes employing the intermediate state formalism: theory, capabilities, and interpretation. Chem A Eur J. 2023;127(32):6635–46. https://doi.org/10.1021/acs.jpca.3c02761 86. Scheurer M, Herbst MF, Reinholdt P, Olsen JMH, Dreuw A, Kongsted J. Polarizable embedding combined with the algebraic diagram- matic construction: tackling excited states in biomolecular systems. J Chem Theory Comput. 2018;14(9):4870–83. https://doi.org/10. 1021/acs.jctc.8b00576 87. Scheurer M, Reinholdt P, Kjellgren ER, Haugaard Olsen JM, Dreuw A, Kongsted J. CPPE: an open-source C++ and python library for polarizable embedding. J Chem Theory Comput. 2019;15(11):6154–63. https://doi.org/10.1021/acs.jctc.9b00758 88. Di Remigio R, Steindal HA, Mozgawa K, Weijo V, Cao H, Frediani L. PCMSolver/pcmsolver: Fluctuating charges MM. PCMSolver, an open-source library for the polarizable continuum modelelectrostatic problem, written by R. Di Remigio, L. Frediani and contributors; 2020. https://doi.org/10.5281/ZENODO.1156166; [cited 2020 Nov 30]. Available from: http://pcmsolver.readthedocs.io/ 89. Lunkenheimer B, Kohn A. Solvent effects on electronically excited states using the conductor-like screening model and the second- order correlated method ADC(2). J Chem Theory Comput. 2012;9(2):977–94. https://doi.org/10.1021/ct300763v 90. JUSTUS2 bwHPC [cited 2020 Nov 30]. https://wiki.bwhpc.de/e/JUSTUS2 91. Nottoli M, Herbst MF, Mikhalev A, Jha A, Lipparini F, Stamm B. Replication data for: “ddX: polarizable continuum solvation from small molecules to proteins.” Version V1. 2024 https://doi.org/10.18419/darus-4030 92. Wang J, Cieplak P, Kollman PA. How well does a restrained electrostatic potential (RESP) model perform in calculating conforma- tional energies of organic and biological molecules? J Comput Chem. 2000;21(12):1049–74. https://doi.org/10.1002/1096-987x(200009) 21:12<1049::aid-jcc3>3.0.co;2-f 93. Hall PR, Malone LD, Sillerud LO, Ye C, Hjelle BL, Larson RS. Characterization and NMR solution structure of a novel cyclic pentapep- tide inhibitor of pathogenic hantaviruses. Chem Biol Drug Des. 2007;69(3):180–90. https://doi.org/10.1111/j.1747-0285.2007.00489.x 94. Ozaki H, Sato T, Kubota H, Hata Y, Katsube Y, Shimonishi Y. Molecular structure of the toxin domain of heat-stable enterotoxin pro- duced by a pathogenic strain of Escherichia coli. A putative binding site for a binding protein on rat intestinal epithelial cell mem- branes. J Biol Chem. 1991;266(9):5934–41. https://doi.org/10.1016/s0021-9258(19)67688-x 95. Yingqi X, Wu J, Pei J, Shi Y, Ji Y, Tong Q. Solution structure of BmP02, a new potassium channel blocker from the venom of the Chi- nese ScorpionButhus martensiKarsch. Biochemistry. 2000;39(45):13669–75. https://doi.org/10.1021/bi000860s 96. Day MW, Hsu BT, Joshua-Tor L, Park J-B, Zhou ZH, Adams MWW, et al. X-ray crystal structures of the oxidized and reduced forms of the rubredoxin from the marine hyperthermophilic archaebacterium pyrococcus furiosus. Protein Sci. 1992;1(11):1494–507. https://doi. org/10.1002/pro.5560011111 97. Chen K, Hirst J, Camba R, Bonagura CA, Stout CD, Burgess BK, et al. Atomically defined mechanism for proton transfer to a buried redox centre in a protein. Nature. 2000;405(6788):814–7. https://doi.org/10.1038/35015610 98. Schanstra JP, Ridder IS, Heimeriks GJ, Rink R, Poelarends GJ, Kalk KH, et al. Kinetic characterization and x-ray structure of a mutant of Haloalkane dehalogenase with higher catalytic activity and modified substrate range. Biochemistry. 1996;35(40):13186–95. https:// doi.org/10.1021/bi961151a 99. Dreveny I, Gruber K, Glieder A, Thompson A, Kratky C. The hydroxynitrile lyase from almond. Structure. 2001;9(9):803–15. https:// doi.org/10.1016/s0969-2126(01)00639-6 100. Valegard K, Andralojc PJ, Haslam RP, Pearce FG, Eriksen GK, Madgwick PJ, et al. Structural and functional analyses of Rubisco from arctic diatom species reveal unusual posttranslational modifications. J Biol Chem. 2018;293(34):13033–43. https://doi.org/10.1074/jbc. ra118.003518 101. Ban N, McPherson A. The structure of satellite panicum mosaic virus at 1.9 A resolution. Nat Struct Biol. 1995;2(10):882–90. https:// doi.org/10.1038/nsb1095-882 102. Gehrmann J, Daly NL, Alewood PF, Craik DJ. Solution structure of α-conotoxin ImI by 1H nuclear magnetic resonance. J Med Chem. 1999;42(13):2364–72. https://doi.org/10.1021/jm990114p 103. Ko T-P, Robinson H, Gao Y-G, Cheng C-HC, DeVries AL, Wang AH-J. The refined crystal structure of an eel pout type III antifreeze protein RD1 at 0.62-A resolution reveals structural microheterogeneity of protein and solvation. Biophys J. 2003;84(2):1228–37. https:// doi.org/10.1016/s0006-3495(03)74938-8 104. Reusken A, Stamm B. Analysis of the Schwarz domain decomposition method for the conductor-like screening continuum model. SIAM J Numer Anal. 2021;59(2):769–96. https://doi.org/10.1137/20m1342872 105. Ciaramella G, Hassan M, Stamm B. On the scalability of the Schwarz method. SMAI J Comput Math. 2020;6:33–68. https://doi.org/10. 5802/smai-jcm.61 106. Hornak V, Abel R, Okur A, Strockbine B, Roitberg A, Simmerling C. Comparison of multiple Amber force fields and development of improved protein backbone parameters. Proteins Struct Funct Bioinform. 2006;65(3):712–25. https://doi.org/10.1002/prot.21123 107. Ambrosetti M, Skoko S, Giovannini T, Cappelli C. Quantum mechanics/fluctuating charge protocol to compute solvatochromic shifts. J Chem Theory Comput. 2021;17(11):7146–56. https://doi.org/10.1021/acs.jctc.1c00763 108. Kovalenko SA, Schanz R, Farztdinov VM, Hennig H, Ernsting NP. Femtosecond relaxation of photoexcited para-nitroaniline: solvation, charge transfer, internal conversion and cooling. Chem Phys Lett. 2000;323(3–4):312–22. https://doi.org/10.1016/s0009-2614(00) 00432-2 109. Stahelin M, Burland DM, Rice JE. Solvent dependence of the second order hyperpolarizability in p-nitroaniline. Chem Phys Lett. 1992; 191(3–4):245–50. https://doi.org/10.1016/0009-2614(92)85295-l NOTTOLI ET AL. 23 of 24 https://doi.org/10.1021/acs.jpca.3c02761 https://doi.org/10.1021/acs.jctc.8b00576 https://doi.org/10.1021/acs.jctc.8b00576 https://doi.org/10.1021/acs.jctc.9b00758 https://doi.org/10.5281/ZENODO.1156166 http://pcmsolver.readthedocs.io/ https://doi.org/10.1021/ct300763v https://wiki.bwhpc.de/e/JUSTUS2 https://doi.org/10.18419/darus-4030 https://doi.org/10.1002/1096-987x(200009)21:12%3C1049::aid-jcc3%3E3.0.co;2-f https://doi.org/10.1002/1096-987x(200009)21:12%3C1049::aid-jcc3%3E3.0.co;2-f https://doi.org/10.1111/j.1747-0285.2007.00489.x https://doi.org/10.1016/s0021-9258(19)67688-x https://doi.org/10.1021/bi000860s https://doi.org/10.1002/pro.5560011111 https://doi.org/10.1002/pro.5560011111 https://doi.org/10.1038/35015610 https://doi.org/10.1021/bi961151a https://doi.org/10.1021/bi961151a https://doi.org/10.1016/s0969-2126(01)00639-6 https://doi.org/10.1016/s0969-2126(01)00639-6 https://doi.org/10.1074/jbc.ra118.003518 https://doi.org/10.1074/jbc.ra118.003518 https://doi.org/10.1038/nsb1095-882 https://doi.org/10.1038/nsb1095-882 https://doi.org/10.1021/jm990114p https://doi.org/10.1016/s0006-3495(03)74938-8 https://doi.org/10.1016/s0006-3495(03)74938-8 https://doi.org/10.1137/20m1342872 https://doi.org/10.5802/smai-jcm.61 https://doi.org/10.5802/smai-jcm.61 https://doi.org/10.1002/prot.21123 https://doi.org/10.1021/acs.jctc.1c00763 https://doi.org/10.1016/s0009-2614(00)00432-2 https://doi.org/10.1016/s0009-2614(00)00432-2 https://doi.org/10.1016/0009-2614(92)85295-l 110. Millefiori S, Favini G, Millefiori A, Grasso D. Electronic spectra and structure of nitroanilines. Spectrochim Acta A Mol Spectrosc. 1977;33(1):21–7. https://doi.org/10.1016/0584-8539(77)80143-8 SUPPORTING INFORMATION Additional supporting information can be found online in the Supporting Information section at the end of this article. How to cite this article: Nottoli M, Herbst MF, Mikhalev A, Jha A, Lipparini F, Stamm B. ddX: Polarizable continuum solvation from small molecules to proteins. WIREs Comput Mol Sci. 2024;14(4):e1726. https://doi. org/10.1002/wcms.1726 24 of 24 NOTTOLI ET AL. https://doi.org/10.1016/0584-8539(77)80143-8 https://doi.org/10.1002/wcms.1726 https://doi.org/10.1002/wcms.1726 ddX: Polarizable continuum solvation from small molecules to proteins 1 INTRODUCTION 2 METHODS 2.1 Model input: Solute cavity and density 2.2 Models for electrostatic interaction 2.3 Discrete equations in a general framework 2.4 Force computations and derivatives 3 IMPLEMENTATION 3.1 Steps of a calculation 3.2 Step 1: Initialization 3.3 Step 2: Building the RHSs for multipolar solutes 3.4 Step 3: Solving the problem 3.4.1 Finer control 3.5 Step 4: Solute specific contributions to the derivatives 3.6 Step 5: Finalization 3.7 FMM acceleration 3.8 Software distribution and build system 4 INTERFACING 4.1 Tinker-interface 4.2 Psi4-interface 4.3 Employing ddX-solvated ground states in adcc 4.4 Standalone driver 5 EXAMPLE APPLICATIONS 5.1 Benchmarks 5.2 Molecular dynamics 5.3 Psi4 calculations 5.4 ADC(2) energies 6 CONCLUSION AND OUTLOOK AUTHOR CONTRIBUTIONS FUNDING INFORMATION CONFLICT OF INTEREST STATEMENT DATA AVAILABILITY STATEMENT RELATED WIREs ARTICLES REFERENCES