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 strategy is based on a boundary element method approach, where the generalized Poisson equation is rewritten as an integral 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, differentiable definitions of the solvation energy. These can be used to compute analytical derivatives of the molecular energy, e.g. to compute response properties or derivatives with respect to the nuclear positions (forces and harmonic frequencies). In quantum chemistry, while more sophisticated models like the linearized Poisson-Boltzmann exist , the two most popular models are based either on a dielectric continuum (Polarizable Continuum Model, PCM) or on a conductor continuum (Conductor-like Screening Model, COSMO, also known as C-PCM) . Such models have been implemented in virtually every quantum chemistry package and are available in conjunction with most levels of theory, from semiempirical methods to sophisticated post-Hartree-Fock ones and also with a linear scaling implementation .