Quantum Simulation of Molecules in Solution

Quantum chemical calculations on quantum computers have been focused mostly on simulating molecules in the gas phase. Molecules in liquid solution are, however, most relevant for chemistry. Continuum solvation models represent a good compromise between computational affordability and accuracy in describing solvation effects within a quantum chemical description of solute molecules. In this work, we extend the variational quantum eigensolver to simulate solvated systems using the polarizable continuum model. To account for the state dependent solute–solvent interaction we generalize the variational quantum eigensolver algorithm to treat non-linear molecular Hamiltonians. We show that including solvation effects does not impact the algorithmic efficiency. Numerical results of noiseless simulations for molecular systems with up to 12 spin-orbitals (qubits) are presented. Furthermore, calculations performed on a simulated noisy quantum hardware (IBM Q, Mumbai) yield computed solvation free energies in fair agreement with the classical calculations.


UCCSD Ansatz
Here we report additional results with a noiseless simulation concerning the calculation of the ground state energy of the HeH + molecule at the equilibrium distance in dimethyl sulfoxide (DMSO) (Fig. S1) using the PCM-VQE.
We have run these calculations to test the hybrid algorithm also on a well established ansatz such as the Unitary Coupled Cluster with Singles and Doubles (UCCSD). We recall that in this approach the optimization parameters coincide with the cluster amplitudes of the single and double excitation operators: T 1 (θ) = i>j θ ij (a † j a i − a † i a j ) (S1)  H + 3 calculations with PCM-VQE/6-31G Here we report additional results on the H + 3 molecule using the 6-31G basis set in DMSO (see Fig. S2). To perform the calculations we have adopted an adaptive ansatz including only the relevant excitations that contribute to the ground state wavefunction as described in the computational details' section of the main text. We have also reported the circuit structure in Fig. S6. As we can see in Fig. S2 the adaptive ansatz is able to reproduce exactly the energy in vacuum. Indeed we recall that for a two electron system the CCSD calculation is equivalent to Full Configuration Interaction.
Moving to the results for the solvated molecule (right panel) again we can see a very good agreement between the PCM-CCSD and the PCM-VQE calculations but in contrast with the gas phase result the numerical values are not exactly matching. Particularly, the PCM-VQE free energy is 8 meV lower. This is due to the fact that for the solvated system the PCM-CCSD calculation is performed approximating the reaction field at the HF level (PTE approximation in the original reference 5 ) as implemented in Psi4. 6 Interestingly this effect is not noticeable in the STO-3G results of the main text because adopting a smaller basis set the error due to the PTE approximation is negligible.
Finally, it is worth to comment the convergence rate of this calculations w.r.t. the results obtained with a 6-31G basis set in the previous subsection (Fig. S1). As expected, the number of iterations and the convergence rate of the optimization procedure is dramatically improved using the adaptive ansatz instead of a less chemically-aware UCCSD circuit. This is

Additional results in aqueous solution
In this section we show additional calculations performed to test the PCM-VQE on the same systems studied in the main text but with a different solvent, i.e. water. These results, in accordance with what already discussed for the case of a DMSO solution, allow us to discuss the performance of the algorithm on another very common solvent in the chemical practice.

Single point calculations
In Tab. S1 we report the results obtained running a PCM-VQE calculation in water. Analogously to what discussed in the main text for the case of DMSO, with the PCM-VQE we are able to reproduce the exact results reported for the CCSD-PCM calculation of the H + 3 molecule and to recover a value of free energy in solution for the larger molecules (BeH 2 and H 2 O) close to the total FCI energy comprehensive of the polarization energy contribution.
All the calculations have been carried out according to the computational protocol described in the computational details' section referring to the results shown in the main text.

Energy and free energy plot for the H 2 O double dissociation profile
Here we report additional results regarding our analysis on the polarization energy and hydration free energy discussed in the main text for the double bond dissociation profile of the water molecule. In Fig. S3a-b we report, respectively, the total energy ( Finally, we also notice that the rate of convergence of the optimization slows down considerably as the bond distance increases. This behavior is not surprising and can be explained considering (i) that as the dissociation regime approaches the initial guess state |HF⟩ provides a worse approximation of the optimal result and therefore a greater number of iterations is required to achieve convergence. (ii) The adaptive ansatz here used includes excited configurations among the most relevant to capture electronic correlation at the equilibrium geometry whose electronic structure is increasingly different from the one at the bond dissociation limit.

Measurement budget allocation and statistical errors
In this section we give a brief account of the effects of finite sample size on the accuracy of the expectation value estimates performed during the VQE and PCM-VQE procedure.
The goal of this section is two-fold: on one hand we want to provide a practical recipe to estimate the number of shots N needed to achieve a desired accuracy on solvation quantities (and, in general, on any expectation value calculated from the VQE wavefunction); on the other hand we want to rationalize the magnitude of the errors at fixed number of shots for different systems and solvation quantities.
The analysis we report explicitly refers to an error due to finite sampling for which it is exact to consider the variances of the expectation values. We point out that in presence of quantum computer noise it would be more appropriate to work explicitly with the mean squared error MSE(X) defined as where X is a random variable representing the expectation value, Var(X) is its variance and Bias(X) comprises the error due to the presence of the quantum computer noise.
We first consider the problem of estimating the number of optimal shots N needed to estimate the expectation value of an observableÔ within a desired accuracy ϵ. We follow Rubin et al. 11 who derived an optimal number of measurement for ⟨Ô⟩ = γ h γ P γ given by where the factor K is obtained taking into account the fluctuation of each independently measured Pauli word P α needed to compute the desired expectation value: Here B is a label that runs over the different sets of non-commuting Pauli words and h α , h β are the corresponding matrix elements of the mapped operatorÔ.
We highlight that Eq. S4 is a direct consequence of error propagation theory and formalizes two important statements: (i) the number of shots increases as the number of independently measured Pauli words increases, (ii) the number of shots increases as the magnitude of the matrix elements ofÔ increases.
In order to use Eq. S4 in practice we need to replace Covar(P α , P β ) with a reasonable estimate. A common choice, 12,13 since Var(P α ) = 1 − ⟨P α ⟩ 2 , is to set all the variances to 1 and all the covariances to zero. It is possible to show that this corresponds to assuming that the Pauli strings are sampled from a uniform distribution. 14 With this assumption, K (and in turn the number of shots N ) will depend only on the Frobenius norm of the operator ||O||: To use what reported here in practice, if a given precision ϵ is sought, one has first to estimate the factor K summing all the squaredÔ matrix elements corresponding to noncommuting Pauli strings and then to divide by ϵ 2 to obtain the corresponding number of samples N to collect. Now we shall consider the error estimated on the total energy and total free energy. This is also relevant for the error on solvation free energies, which are obtained as differences of total free energy in solution and total energy in gas-phase. Suppose for example that an error of ϵ = 0.1 Ha on the total energy (or free energy) for water at the STO-3G level is desired. The first step is to calculate K = ||H|| 2 for water from eq. S5, that turns out to be ≈ 6 · 10 3 . Then, eq. S3 gives a number of samples N ≈ 6 · 10 5 .
The presented discussion can also be used to numerically explore how error estimates depend on the specific calculation, via the Hamiltonian norm and the number of non-commuting sets of Pauli words at fixed number of shots, see Tab. S2.
As we can see from such Table, the error increases both as the number of non commuting Pauli terms and the Hamiltonian norm increases. This is in accordance with Eq. S3 and Eq.

S5.
Turning now to the error on the polarization energy, that we identify in the main text as a robust solvation quantity also for noisy quantum computers, we note that the same procedure just described to estimate N to obtain a desired accuracy ϵ can be applied directly to the operator V σ . Moreover, by using Eq. S4 we can also rationalize the different accuracy observed in estimating the solvation free energy and the polarization energy. Ultimately this should be ascribed to the fact that the number of independently measured Pauli strings is much smaller for the polarization energy alone than for the overall solvation free energy since in one case we deal with a one-electron operator while in the latter case we have to include all the measurements coming from the two-electrons part. More specifically, in Tab. S3 we show estimated sampling errors for 1 2 V σ (the operator needed to compute the polarization energy). The comparison of the magnitude of the errors reported in this Table with those in Tab. S2 confirms our qualitative analysis. Interestingly, the error for the smallest molecule H + 3 here is the largest in the Table. This is not surprising, since H + 3 is the only charged molecule in the set, and as a consequence the V σ norm is the largest.

Algorithmic complexity: PCM classical overhead
The proposed PCM-VQE algorithm is an hybrid one, combining a classical portion with a quantum one. In view of realistic simulations on future quantum computers, it is important to assess that the extension to PCM does not hamper potential quantum advantage from the VQE approach. The goal of this section is not to provide a thorough analysis of the algorithmic cost, which depends on the particular flavour of VQE considered, but to highlight the computational overhead due to the inclusion of solvation effects within the PCM framework.
Since we are considering an hybrid algorithm we estimate the overall cost C = C Q + C C where C Q is the cost associated with operations run on the quantum computer and C C the classical cost to compute the molecular integrals .
Concerning C Q , we note that gas-phase VQE requires evaluation of 1-and 2-RDM to obtain the expectation value of the molecular Hamitonian ⟨Ĥ 0 ⟩ (see Figure 1 in the main text). The 1-RDM is the only information from the quantum computer that is also needed to calculate ⟨V σ ⟩ , i.e., the additional part with respect to ⟨Ĥ 0 ⟩ to calculate the total free energy. Since the 1-RDM is calculated for the standard VQE as well, there is no quantum overhead for PCM-VQE with respect to gas-phase VQE (i.e., C gas Q = C PCM Q ).
Now we move our analysis to the classical cost of the algorithm. Here the inclusion of solvation effects may potentially affect the algorithmic cost. Particularly, the standard molecular Hamiltonian is usually built computing the molecular one-and two-electron integrals. 15 On the other hand, inclusion of solvation effects calls for (i) computing the solvent response at each iteration q and (ii) calculating the solute-solvent interaction term ⟨V σ ⟩ at each iteration. (i) Formally, computing the response charges implies the inversion of the PCM response matrix whose computational cost is O(N 3 tess ). In practice it is possible to achieve a linear scaling with N tess for this step adopting techniques such as the Fast Mul-tipole Method (FMM) and parallelization. 16 In turn, N tess is scaling with the size of the molecular surface, that is linear with the number of atoms N atoms in the worst case of linear molecules. This is a negligible scaling w.r.t. the following contribution. (ii) Concerning the calculation of ⟨V σ ⟩, at each step we have to contract the 1-RDM with the v rs array a step costing O(N 2 a N tess ) = O(N 2 a N atoms ) additional operations. Reasonably assuming that N a is scaling linearly with N atoms , this step has a O(N 3 a ) scaling. This is still a better scaling than the cost of the contraction of the 2-RDM with the gas-phase bielectronic integrals (that dominates the gas-phase C C ), where a sum over four index is needed (O(N 4  a )). As such the scaling of C PCM C will be at most equal to that of C gas C . Of course, any strategy to improve such O(N 4 a ) scaling (e.g., use of molecular symmetry, prescreening of the molecular integrals) can also be used to improve the scaling of the classical PCM term.
In conclusion, the overhead of the PCM extension to VQE is null for the quantum part, and does not worsen the scaling of the classical part with respect to gas-phase calculations; a potential quantum advantage of gas-phase VQE is therefore unhampered by the inclusion of the PCM solvent. The bottleneck remains to be the solution of the electronic structure problem.

Quantum circuits
(a) (b) Figure S5: Quantum circuits used for the PCM-VQE simulations for the BeH 2 (a) and H 2 O (b) molecules. Please note the difference w.r.t. the circuit reported in the main text for H + 3 . The space of multiple excitations is realised by applying more Givens rotations (18 and 30 for BeH 2 and H 2 O respectively). We also note the presence of single excitation operators missing in the H + 3 circuit as a result of the adaptive procedure that discards irrelevant excitations. Labels G and G 2 refer to single-and double-excitation unitaries implemented as rotations in the subspace spanned by two or four qubits. Figure

TOC Graphic
Some journals require a graphical entry for the Table of Contents. This should be laid out "print ready" so that the sizing of the text is correct. Inside the tocentry environment, the font used is Helvetica 8 pt, as required by Journal of the American Chemical Society. The surrounding frame is 9 cm by 3.5 cm, which is the maximum permitted for Journal of the American Chemical Society graphical table of content entries. The box will not resize if the content is too big: instead it will overflow the edge of the box. This box and the associated title will always be printed on a separate page at the end of the document.