Quantum Monte Carlo Treatment of the Charge Transfer and Diradical Electronic Character in a Retinal Chromophore Minimal Model.

The penta-2,4-dieniminium cation (PSB3) displays similar ground state and first excited state potential energy features as those of the retinal protonated Schiff base (RPSB) chromophore in rhodopsin. Recently, PSB3 has been used to benchmark several electronic structure methods, including highly correlated multireference wave function approaches, highlighting the necessity to accurately describe the electronic correlation in order to obtain reliable properties even along the ground state (thermal) isomerization paths. In this work, we apply two quantum Monte Carlo approaches, the variational Monte Carlo and the lattice regularized diffusion Monte Carlo, to study the energetics and electronic properties of PSB3 along representative minimum energy paths and scans related to its thermal cis–trans isomerization. Quantum Monte Carlo is used in combination with the Jastrow antisymmetrized geminal power ansatz, which guarantees an accurate and balanced description of the static electronic correlation thanks to the multiconfigurational nature of the antisymmetrized geminal power term, and of the dynamical correlation, due to the presence of the Jastrow factor explicitly depending on electron–electron distances. Along the two ground state isomerization minimum energy paths of PSB3, CASSCF calculations yield wave functions having either charge transfer or diradical character in proximity of the two transition state configurations. Here, we observe that at the quantum Monte Carlo level of theory, only the transition state with charge transfer character can be located. The conical intersection, which becomes highly sloped, is observed only if the path connecting the two original CASSCF transition states is extended beyond the diradical one, namely by increasing the bond-length-alternation (BLA). These findings are in good agreement with the results obtained by MRCISD+Q calculations, and they demonstrate the importance of having an accurate description of the static and dynamical correlation when studying isomerization and transition states of conjugated systems.

Introduction

The retinal protonated Schiff base (RPSB, represented in Figure 1a) is the chromophore responsible for the photochemical properties of a vast family of biological photoreceptors referred to as retinal proteins which are, among other functions, involved in the mechanism of vision of dim light in vertebrates.[1−3] RPSB undergoes a very fast cis–trans isomerization in the protein (opsin) environment (∼200 fs) with high quantum yield (∼0.68) upon photon absorption.[1] This isomerization process has been deeply investigated by different spectroscopic techniques[4] and theoretical calculations.[5,6] In particular, femtosecond spectroscopy[7] and hybrid quantum mechanics/molecular mechanics (QM/MM) molecular dynamics calculations[7,8] highlight the essential role played by the surrounding protein environment and by the hydrogen-out-of-plane motion[9] in the isomerization mechanism, which involves decay via a S1/S0 conical intersection (CI)[10] reached by the selective torsion of the central (C11=C12) double bond.

Figure 1

(a) Lewis representation of the full chromophore, the 11-cis retinal protonated Schiff base (RPSB), covalently bound to a lysine of the opsin protein, and of the cis-,trans-penta-2,4-dieniminium cations (PSB3). (b) Schematic representation of the thermal isomerization from the cis to the trans configuration of PSB3, following the two minimum energy paths (obtained in ref (11) from CASSCF calculations) undergoing respectively a charge transfer (CT) and a diradical (DIR) transition state. The former minimum energy path (MEPCT) is represented in red with a saddle point in the charge transfer transition state (TSCT), and the orange region indicates that the ground state potential energy surface is characterized by a charge transfer character. The latter minimum energy path (MEPDIR) is represented in blue with a saddle point in the diradical transition state (TSDIR), and the cyan region indicates that the ground state potential energy surface is characterized by a diradical character. The energy barriers corresponding to TSCT and TSDIR are comparable, and most of the computational studies agree that TSCT represents the lowest barrier.[11−13] The main difference between TSCT and TSDIR configurations lies in their bond length alternation (BLA), which is positive for the TSDIR and negative for the TSCT (see discussion in the text). A scan of the configurations connecting TSCT with TSDIR (called BLA scan) presents a conical intersection (CI) for a positive BLA, but the actual value varies with the adopted computational method.[11−14]

The penta-2,4-dieniminium cation (PSB3, represented in Figure 1a) is a small conjugated molecular system widely used as reduced computational model of the full RPSB because PSB3 and RPSB exhibit similar ground and excited state features. First, both PSB3 and RPSB have equilibrium cis and trans isomers where the positive charge is localized on the nitrogen-containing side of the conjugated chain. Second, the transition from the ground state (S0) to the first excited state (S1) is characterized by a transfer of the positive charge from the nitrogen region toward the opposite end of the conjugated chain. Third, twisting the conjugated chain along one of the central double bonds leads to a CI structure between the S0 and S1 electronic states. Such a CI plays an important role in the photoisomerization mechanism of RPSB because it mediates population transfer from S1 to S0 along the isomerization coordinate.

Not only do PSB3 and RPSB have similar near-equilibrium and S1 potential energy surface features, but it was recently found that they have similar ground state potential energy features in the vicinity of the CI as well. In a hybrid quantum mechanics/molecular mechanics (QM/MM) study of the thermal isomerization mechanism of bovine rhodopsin, it was found that thermal 11-cis to all-trans isomerization of RPSB may occur via one of two distinct saddle points that are in the vicinity of the S0/S1 CI.[2] Both saddle points (i.e., chemically, transition states) feature an almost orthogonally twisted C11=C12 double bond, similar to the CI but very different bond length alternations (BLA). We define the BLA as the difference between the average bond length of formal single bonds (C1–C2 and C3–C4 for PSB3) and the average bond length of formal double bonds (C1=N, C2=C3, and C4=C5 for PSB3), such that BLA is positive if there is no bond inversion and negative if there is. The situation is analogous in the case of PSB3. Indeed, PSB3 also features two transition states which are almost 90° twisted along the central (C2=C3) double bond but having different BLA patterns. In both PSB3 and RPSB, the two transition states not only have different BLA geometries but also different electronic structures. One transition state (TSCT) is characterized by a transfer of the positive charge from the nitrogen-containing side of the molecule to the other end of the conjugated chain. The other transition state (TSDIR) retains the positive charge on the nitrogen side of the chain and therefore has covalent/diradical character due to the homolytic breaking of the isomerizing double bond. To characterize the regions of the ground state potential energy surface driving the thermal isomerization, Gozem et al.[11] optimized the two transition states in PSB3 and used them to map three pathways at the CASSCF/6-31G* level of theory. The first two pathways are minimum energy paths (MEPs) leading away from each transition state and toward the two equilibrium structures, cis-PSB3 and trans-PSB3. These paths are called MEPCT (for the path passing through the charge transfer transition state TSCT) and MEPDIR (for the path passing through the covalent/diradical transition state TSDIR). The third path is an interpolation/extrapolation of coordinates between the two transition states. Because the transition states essentially have different BLAs but otherwise are almost geometrically identical, this scan follows a BLA coordinate and is therefore called the BLA scan. This scan intercepts a S1/S0 CI structure. For a simplified representation of the three paths, see Figure 1b; further information are reported in ref (11).

The three aforementioned paths served useful as benchmarks to understand the importance of a balanced representation of electronic correlation in correctly describing the topology of the potential energy surface in that region as well as in computing the relative energies of TSCT and TSDIR. If the CASSCF(6,6)/6-31G* can be considered a good compromise between accuracy and computational effort for mapping (small) parts of the potential energy surface of complex molecular systems, using more sophisticated methodologies for single-point energy calculations on the CASSCF paths is seen to be a reliable strategy to explicitly include dynamical correlation in the wave function,[11−13,15] allowing one to correct the original CASSCF findings such as CI location, energy of the π2 and ππ* configurations, and relative stability of the TSCT and TSDIR.

Thanks to the tremendous progress in high performance computing (HPC) in recent years, quantum Monte Carlo (QMC) methods represent a powerful alternative to other ab initio and DFT approaches for the accurate description of systems where electronic correlation plays an essential role. QMC methods[16,17] have been widely applied to several problems of physical and chemical interest such as materials,[18−22] molecular properties,[23−34] reaction pathways,[35,36] and geometry and excited states of biochromophores.[37−41] The good scalability with respect to the system size (N, with 3 < d < 4 and N the number of electrons)[16,39] and the use of massively parallel algorithms make QMC methods particularly suitable for Petascale architectures. All these aspects justify the growing number of applications of QMC in problems of quantum chemistry and molecular physics.

The variational Monte Carlo (VMC) method[42] is the simplest QMC approach: thanks to a combined use of the Monte Carlo integration and the variational principle for the ground state, the many-body trial wave function Ψ can be efficiently optimized and used to extract information other than just the energy of the molecular target. Higher accuracy in the determination of the system properties (energy, geometrical and electronic parameters, etc.) is usually achieved by the fixed node projection Monte Carlo methods such as the diffusion Monte Carlo (DMC)[43,44] or the lattice regularized diffusion Monte Carlo[45,46] (LRDMC).

The choice of the functional form of Ψ represents a crucial step determining the overall quality of QMC calculations, both in the variational and in the fixed node projection schemes. The Jastrow antisymmetrized geminal power (JAGP)[47−52] has been seen to be very efficient in the investigation of chemical systems,[28−30,32,33,35,39−41] with accuracy comparable to that of high-level quantum chemistry methods. Its compactness, coupled to the use of efficient algorithms for the optimization of all parameters, including linear coefficients and exponents of the atomic basis set,[25,53,54] leads to a fast convergence of the variational results for electronic and geometrical properties with the size of the basis set,[28,55,56] with a computational cost comparable to that of a simple wave function defined by the product of a Jastrow factor and a single Slater determinant. Although AGP is not in general size consistent, the presence of a flexible Jastrow factor makes JAGP size consistent for systems which splits in fragments with spin zero or 1/2.[25,49,57] JAGP has already proven to give a good description of the static and dynamical correlation in several crucial benchmarks[32,50,51,58] like in the estimation of the torsional energy of the ethylene and of the singlet–triplet gap of methylene.[32]

In the following, we compute the energy and electronic structure of PSB3 along the MEPCT, MEPDIR, and BLA paths reported previously.[11] As shown below, the combination of QMC and JAGP ansatz allows us to get a proper description of both static and dynamical correlation; for this reason, using the QMC and JAGP wave functions represents an optimal choice to study the intrinsic properties of the PSB3 model.

The paper is organized as follows: in section 2 we report the main features of VMC and LRDMC schemes and a detailed analysis of the JAGP ansatz and of its potential to correctly describe multiconfigurational systems such as diradicals; computational details on our calculations are given in section 3; the current results are shown in section 4, pointing out the importance of the electronic correlation properly introduced by the QMC/JAGP study for the characterization of the conical intersection; conclusions and comments on future perspectives are reported in the last section.

Quantum Monte Carlo

Variational and Lattice Regularized Diffusion Monte Carlo

The accuracy of QMC approaches, both in the simplest VMC scheme and in the fixed-node projection schemes, are strictly related to the wave function ansatz. Typically, the electronic wave function Ψ in QMC[16,17,56] is defined by the productwhere is the antisymmetric function taking into account the Fermionic nature of electrons and is the Jastrow factor depending explicitly on the interparticle (electrons and nuclei) distances; x̅ and R̅ represent the collective electronic (x̅ refers to space r̅ and spin σ̅) and nuclear coordinates, respectively.

The Jastrow factor is a symmetric positive function of the electronic positions; therefore it does not change the nodal surface (determined by the antisymmetric term ), but it introduces the dynamical correlation among electrons and satisfies the electron–electron and electron–nucleus cusp conditions.[16,56,59]

In VMC, the parameters that define Ψ are optimized in order to minimize the electronic energy within the functional freedom of the ansatz. The VMC results can further be improved by using the fixed-node (FN) projection Monte Carlo techniques, which provide the lowest possible energy with the constraint that the wave function ΦFN has the same nodal surface of an appropriately chosen guiding function Ψ,[16,43] which is usually optimized using the VMC method. The fixed-node projection Monte Carlo method that we have adopted is LRDMC,[45,46] which is efficient for systems with a large number of electrons[46] and preserves the variational principle even when used in combination with nonlocal pseudopotentials.[46] Because the LRDMC calculations are much more demanding than the VMC calculations, in terms of computational time, they have been performed only for a few key structures.

The Jastrow Antisymmetrized Geminal Power

The trial wave function ansatz used in the QMC calculations presented in this paper is the Jastrow antisymmetrized geminal power (JAGP),[47,48,56] that is, the productof the antisymmetrized geminal power (AGP) function ΨAGP and the Jastrow factor ΨJ, where the dependence on the nuclear coordinates R̅ is here omitted.

For an unpolarized system (zero total spin S) of N = 2Np electrons and M atoms, the AGP function is defined aswhere is the antisymmetrization operator and the geminal pairing function G is a product of a singlet function and a spatial wave function symmetric with respect to the particle exchange :The spatial function is a linear combination of products of atomic orbitals ϕμ:where the indexes μ and ν run over all the basis in all the atoms in the system, for a total of L atomic orbitals (note that L is determined by the overall basis set size). The coefficients gμν have to be optimized in order to minimize the variational energy of the system (together with the other parameters in the wave function).

In our calculations, we used this Jastrow factorthat involves the one-electron interaction term Uen, the homogeneous two-electron interaction term Uee, and the inhomogeneous two-electron interaction terms Ueen and Ueenn (representing respectively an electron–electron–nucleus function and an electron–electron–nucleus–nucleus function). They are defined as follows:where the vector r = r – R is the difference between the position of the nucleus a and the electron i, r is the corresponding distance, r is the distance between electrons i and j, Z is the electronic charge of the nucleus a, which is described by l atomic orbitals χμ (with index μ = 1,...,l), and b1, b2, fμ, f¯μ,ν, and f̃μ,μ are variational parameters. [Note that the atomic orbitals χ used in the Jastrow term are similar to the atomic orbitals ϕ used for the AGP, although they are not the same orbitals, and in general a reliable description of molecular systems requires a much smaller number of orbitals in the Jastrow than the number of orbitals used in the AGP, see for instance ref (56).] The leading contribution for the description of electronic correlation is given by Uee, but also the inhomogeneous two-electron interaction terms Ueen and Ueenn are particularly important in the JAGP ansatz because they reduce the unphysical charge fluctuations included in the AGP function, as discussed in refs (25, 49).

The pairing spatial function in eq 5 is written in terms of the (localized) atomic orbitals ϕμ, offering an interesting correspondence between the AGP ansatz and the resonating valence bond framework.[60,61] An equivalent way to write the pairing function is obtained by using the molecular orbitals (MOs) ψ. The expansion of the pairing function in terms of MOs is obtained by performing a generalized (the atomic orbitals ϕμ are not necessarily orthonormal, so the overlap matrix Sμν = ⟨ϕμ|ϕν⟩ ≠ δμν) diagonalization of the coupling matrix G, which is the L × L matrix of the gμν coefficients:where Λ = diag(λ1,...,λ)

andThe resulting pairing function iswhere the orthonormal single particle functions are written aswith the Pμ coefficients defining the eigenvectors P.

It is interesting to investigate the connection between the expansion of the pairing function in terms of MOs and the standard configuration interaction expansion of the wave function in multiconfigurational approaches.

By substitution of eq 12 in eq 3, and expanding the summation out of the antisymmetrization operator, the following multideterminant expansion is obtained for the AGP function:where the coefficients are given by|Ψ0⟩ is the leading closed-shell Slater determinant:the determinant |Ψ⟩ is equal to |Ψ0⟩, but with the virtual orbital ψ substituting the valence orbital ψ, etc. From the expression of the coefficients in eq 15 and the ordering of the eigenvalues λ in eq 11, it follows that the leading contribution beyond the determinant |Ψ0⟩ is given by the determinant |Ψ⟩ with i = Np and a = Np + 1. The multideterminant expansion of ΨAGP in eq 14 allows us to directly compare the ΨAGP with wave functions from other quantum chemical frameworks. In ΨAGP, all the odd excited determinants (single, triple, etc.) are excluded, whereas a subset of the even excitations (those with a multiple excitation to the same virtual orbital) are taken into account; only doubly occupied molecular orbitals are present. In other words, ΨAGP is contained in the seniority zero sector of the electronic full configuration interaction and its expansion coefficients are determined by the ratios of the eigenvalues of the Λ matrix.

The seniority number Ω represents an alternative tool to classify singlet wave functions. Ω is defined as the number of unpaired electrons in the Slater determinant, i.e., the number of singly occupied molecular orbitals. In cases when the static correlation plays a major role Ω-based selection of important Slater determinants in the expansion has been seen to be superior than the traditional one, based on the number of excitations with respect to the reference configuration.[62] Wave functions with Ω = 0 for benchmark systems are accurate enough to recover most of the static correlation, but the FCI limit (including dynamical correlation) is achieved only when configurations from Ω = 2, 4, 6, ... sectors are explicitly included.[62] In the case of JAGP wave function, the combination between a Ω = 0 determinantal term and a Jastrow factor allows us to estimate the correlation energy more accurately than Ω = 0 CI wave functions. The set of MOs ψ is optimized within the JAGP framework, i.e., in the presence of the Jastrow factor and of the multiconfigurational character of the wave function: MOs extracted from our optimization procedure represent therefore the optimal choice for the correlated description of the system under study. The way to move from the AGP MOs to the standard ones and vice versa is explained in the next paragraph and in Appendix.

AGP for Diradicals

The multiconfigurational nature of the AGP function, clearly shown by eq 14, has been extensively discussed by Zen et al. in a recent paper on the use of AGP for diradicals.[32] Here we focus the attention on the application of the AGP ansatz to the ππ* state of PSB3 (see the discussion in section 4), starting from the simple but accurate model with two electrons in two orbitals, originally introduced by Salem and Rowland.[63,64] For a generic diradical system, two atomic orbitals ϕA and ϕB are centered on nuclei A and B. Such a model is a representative scheme for molecules undergoing a double bond breaking, like the twisted ethylene or the twisted PSB3: in these cases, the two involved orbitals are of p type and positioned on two (central) carbon atoms.

The main goal of this discussion is to demonstrate that the AGP ansatz, a formally Ω = 0 wave function, contains the ππ* configuration, needed for a proper description of the ground state surface of PSB3 along the paths. [The π and π* orbitals have to be intended here as π-molecular orbitals residing on highly twisted configurations.]

In case of two orbitals, the pairing function term in the AGP formulation is explicitly given bywhere gμν coefficients represent the coupling terms of the G matrix in the expansion of the AGP spatial factor: gAA and gBB are referred to the ionic terms ϕA(r1)ϕA(r2) and ϕB(r1)ϕB(r2) in which the two electrons are localized on the same atom, whereas the elements gAB and gBA are related to the covalent terms ϕA(r1)ϕB(r2) and ϕB(r1)ϕA(r2) .

In terms of molecular orbitals, becomes (L = 2)On the other hand, in the standard delocalized picture[63,64] two ψ+ and ψ– molecular orbitals are defined:assuming zero overlap between ϕA and ϕB (the assumption is easily verified for two orthogonal p orbitals).

The diradical Ψ(r1,r2) = ψ+ψ– wave function (ψ+ψ– is a short notation for ψ+(r1)ψ–(r2) + ψ–(r1)ψ+(r2)) is written in the atomic basis as the followingBecause the AGP formally contains the terms introduced in eq 20, the comprehension of a consistent way to link together eqs 17, 18, and 20 is a mandatory task. Some questions are therefore arising: (i) which relation occurs between ψ1, ψ2, ψ+, and ψ–; (ii) which relation exists between the atomic basis, ϕA and ϕB, and the molecular orbitals, ψ1 and ψ2 obtained by the diagonalization of the matrix (eq 18); (iii) finally, if a diradical ground state ψ+ψ– can be properly described by the AGP ansatz.

A unitary matrix transforms the ψ+ and ψ– molecular orbitals, by rotating them with a certain angle θ, into ψ̃+ and ψ̃–If θ = π/4, the normalized orbitals becomeBy defining ψ1 ≡ ψ̃+ and ψ2 ≡ ψ̃–, which become ψ1 = ϕA and ψ2 = ϕB by taking into account the relations in eqs 19, and by substitution in eq 18, we obtainwhich is equivalent to eq 20 for λ1 = −λ2 = 1/√2. We also observe that eq 17 is equivalent to eq 20 for gAB = gBA = 0 and gAA = −gBB = 1/√2.

Summarizing, it is always possible to transform a ψ+ψ– configuration into a combination of ψ̃+2 and ψ̃–2 by applying an opportunely chosen proper unitary transformation. In this simple model for diradicals, the AGP molecular orbitals are a linear combination of the orbitals deriving from the traditional picture of electronic delocalization. The detailed derivation of the mapping procedure is reported in the Appendix, where the generalization from the (2,2) to the (2,n) active space (where n is the number of orbitals) is shown and the proof of the capability of AGP to represent such active space is reported. The expansion in eq 17 shows that the AGP ansatz contains all the terms reported by the picture in terms of delocalized molecular orbitals and localized atomic orbitals. The gμν coefficients are variational parameters optimized by the stochastic methods mentioned before, and for this reason the AGP optimization is a fundamental step to select the right wave function for the ground state of interest.

As a conclusion of this discussion, the graphical representation of CASSCF and AGP molecular orbitals for highly twisted diradical ππ* state of a certain structure of PSB3 (the last structure of the BLA path, see the discussion below), together with the combinations explained in this paragraph, is reported in Figure 2: the AGP frontier orbitals (ψ1 and ψ2, eq 18) are identical to the normalized sum and difference of the corresponding CASSCF orbitals (π and π* orbitals), according to the rotation given by the unitary transformation .

Figure 2

Comparison between CASSCF frontier (π and π*) and AGP frontier (ψ1 and ψ2, see eq 18) orbitals; corresponding combinations, given the unitary transformation (eq 21), are also shown.

Computational Details

The QMC calculations reported in this paper have been obtained using the TurboRVB package, developed by S. Sorella and co-workers,[65] that includes a complete suite of variational and diffusion Monte Carlo codes for wave function and geometry optimization of molecules and solids. The scalar-relativistic energy consistent pseudopotentials (ECP) of Burkatzki et al.[66] have been used in order to describe the two core electrons of the carbon and nitrogen atoms. In detail, the basis sets we have used for the AGP part are (10s,9p,2d,1f) contracted in {8} hybrid orbitals for the carbon atom, (8s,9p,2d,1f) contracted in {8} hybrid orbitals for the nitrogen atom, and (6s,5p,1d) contracted in {1} hybrid orbitals for the hydrogen atom. As basis sets for the atomic orbitals included in the inhomogeneous terms of the Jastrow factor, namely in Uen, Ueen and Ueenn reported in eqs 6, 8 and 9, we used an uncontracted basis for the Uen and Ueen term and a contracted with hybrid orbitals basis for Ueenn. This allowed us to have an accurate basis set for the Jastrow factor while keeping the number of parameters of the wave function reasonably small. In more details, in Uen and Ueen we used a (4 s,2p,1d) basis set for the carbon or nitrogen atoms, and a (3 s,2p) for hydrogen atom, whereas in Ueenn the orbitals are contracted in {2} hybrid orbitals for the oxygen, nitrogen, or hydrogen atoms.

The parameters of the wave function Ψ, including the values of the exponents of the atomic orbitals, have been optimized by using the already validated and stable optimization schemes discussed in ref (56). In particular, the optimization that we have followed for the singlet π2 and triplet states of PSB3 considered here in the different structures starts from an initial configuration where the AGP matrix is diagonal, the exponents are initialized to values taken from standard Dunning’s basis sets (where too small and too large values are eliminated because they are not not necessary due to the presence of our Jastrow factor, see discussion in ref (56)), and all the Jastrow parameters are set to zero, with the exception of b1 = b2 = 1. Next, the optimization procedure follows the protocol: (i) optimization of the AGP, namely of the matrix elements and the contraction coefficients of the basis set, with fixed exponents and Jastrow parameters b1 = b2 = 1, (ii) optimization of the AGP and relaxation of the values of the exponents of the AGP basis set and of the b1 and b2 parameters, (iii) optimization of the Jastrow terms, keeping the AGP parameters fixed, (iv) optimization of the overall JAGP, keeping fixed the exponents in the basis set, both for the AGP and the Jastrow, (v) optimization of all the parameters, including the exponents of the basis set, with increasing statistical accuracy. For the single diradical state (ππ* configuration) of PSB3 in the proximity of the conical intersection, we have used a slightly different procedure to avoid the possibility to be trapped in a local minimum. Therefore, we forced the wave function to be in the correct electronic configurations by taking the triplet JAGP optimized wave function and obtained from that the corresponding singlet diradical configuration. This wave function has been used as the starting point of an optimization that started from the step (iii) of the previously stated optimization schedule. We have verified a posteriori for every nuclear structure, where we have calculated both the π2 and the ππ* singlet configurations, that the overlap ⟨ΨJAGPππ*|ΨJAGPπ⟩ between the two JAGP wave functions (calculated using the correlated sampling techniques) is almost zero, thus the two wave functions actually correspond to different electronic states.

It is important to note that QMC approaches use stochastic methods both to evaluate an observable and to optimize the variational parameters of the wave function. The stochastic uncertainty due to the former point is easy to calculate, and it has been reported in figures and tables of the present work. The latter point is instead much more difficult to evaluate. We have carefully tested the reliability of the optimization schemes used in this work, and indeed the profiles reported in the following figures, although not perfectly smooth, are pretty regular, both for energy and charge-transfer values. Moreover, the most interesting configurations (cis, trans, and structures close to TSCT and TSDIR) have been optimized with some extra effort, thus the results reported in the tables are fully reliable.

In this work we also report several results computed at the level of the fixed-node projection Monte Carlo scheme that has been realized by performing LRDMC calculations with mesh size a = 0.3 au. Although we have not performed, for computational reasons, the continuous extrapolation of the lattice mesh size a→ 0, we know from previous works[32,56] and preliminary calculations that the bias given by the finite mesh size a = 0.3 au is almost negligible in the evaluations of the considered energy differences.

Results and Discussion

In this section, we present the energetic and electronic features of the MEPCT, MEPDIR, and BLA paths computed using VMC and LRDMC with the JAGP ansatz. A high-level treatment of electron correlation is crucial for the correct description of the energy surface of the ground state isomerization of PSB3. As discussed in ref (11), dynamical electronic correlation modifies the mapped CASSCF potential energy surface in two ways: the TSCT transition state is found lower in energy than the TSDIR, at variance with the CASSCF findings, and the CASSCF CI is seen shifted to larger BLA values when dynamical correlation is included in the calculations. The energies along the three paths have been computed using a number of electronic structure methods, namely multiconfigurational approaches,[11] multireference perturbation theory,[11] DFT schemes,[13,67] EOM-CC,[12] and SORCI.[15] All the most accurate methods qualitatively predict similar changes in the potential energy surface with respect to CASSCF. In particular, in the case of MRCISD+Q (the most accurate method tested previously), TSCT becomes more stable than TSDIR by 4.7 kcal/mol (compared to CASSCF where it is less stable by 1.2 kcal/mol), and the CI gets shifted to a BLA value of ∼0.03 Å (compared to ∼0.00 Å for CASSCF).

Before starting, we validated our computational protocol by looking at the electronic properties of the cis isomer of PSB3 such as the dipole moment μ and the charge transfer of the ground state S0, defined as the partial charge on the allyl moiety H2C=CH–CH= (the net charge of the system is +1). The charge transfer character at cis-PSB3 is 0.313 at the CASSCF/6-31G* level of theory and 0.355 at the MRCISD+Q level of theory, as derived from Mulliken population analyses[11] (see Table 1). Because we cannot define Mulliken charges in our QMC framework, we compute the charge transfer by finding the portion of the electronic density in the region of the allyl moiety up to the plane perpendicular to the C2=C3 bond and cutting it in the middle. This method of obtaining the charge-transfer character is tested on densities extracted from DFT and wave function methods and is shown to produce very similar charge transfer character as Mulliken charges.

Table 1

Singlet Ground State Energy (in Hartree, H), Total Dipole μ (in Debye, D) and Charge-Transfer Value for the cis Isomer of PSB3, Evaluated with Several Computational Methods, As Defined in the First Columna

method	ref	core	energy [H]	μ [D]	charge transfer
CASSCF(6,6)//6-31G*	(11)	AE			0.313
MRCISD+Q//6-31G*	(11)	AE			0.355
PBE//cc-pVTZ	this work	AE		3.784	0.375
PBE//VTZ-ANO	this work	ECP		3.758	0.369
B3LYP//cc-pVTZ	this work	AE		3.718	0.380
B3LYP//VTZ-ANO	this work	ECP		3.639	0.376
HF//cc-pVTZ	this work	AE		3.441	0.405
HF//VTZ-ANO	this work	ECP	–41.7258	3.472	0.396
VMC/SDb	this work	ECP	–41.7048(9)	3.625	0.365
VMC/J1-bodySDc	this work	ECP	–41.7130(7)	3.633	0.365
VMC/JSD-proj	this work	ECP	–42.8361(2)	3.895	0.363
VMC/JSD-opt	this work	ECP	–42.8373(2)	3.900	0.360
VMC/JAGP	this work	ECP	–42.8490(2)	3.983	0.356
LRDMC/JAGP	this work	ECP	–42.9160(3)	4.066	0.352

J1-bodySD, JSD-proj, and JSD-opt are defined in the text, and the basis sets for the QMC calculations are defined in section 3. In the core column, AE stands for all-electron calculation and ECP for energy-conserving pseudopotential calculation. The reported numbers in the last column represent the net charge on the allyl moiety. VMC and LRDMC errors on μ and charge-transfer are of the order of 10–3.

Wave function optimization by DFT/LDA; EDFT/LDA = −42.6663848 H.

Wave function optimization by DFT/LDA; EDFT/LDA = −42.6769218 H.

The single-reference nature of the S0 state is highlighted by the fact that charge-transfer values computed with standard DFT (with PBE and B3LYP functionals) and HF are in good agreement with the MRCISD+Q result. The effect of applying pseudopotentials on the carbon and nitrogen atoms on the charge transfer is found to be negligible when comparing all-electron (AE) and ECP results obtained using similar basis sets, as shown in Table 1. A further evidence of the reliability of our approach is given by the performance of several variants of the single-determinant (SD) wave function (J1-bodySD, where only the 1-body term for the Jastrow is used; JSD-proj, the SD wave function is projected out from the full AGP; JSD-opt, the SD wave function is optimized after projection) in the VMC framework (e.g., taking into account only the first term of the AGP expansion in eq 14) that, using ECP, are in full agreement with the more accurate VMC/JAGP (0.358, the best variational result) and the MRCISD+Q. LRDMC/JAGP only slightly corrects (0.352) the VMC/JAGP result for the charge transfer. The same conclusions are easily extended to the dipole moment.

MEPCT and MEPDIR Paths

As first, we consider the energy difference between the cis and trans isomers of PSB3. VMC/AGP and LRDMC/AGP values (Table 2), −2.9(2) and −3.0(2) kcal/mol respectively, and the VMC/JSD value (−2.9(2) kcal/mol) are fully consistent with the reference MRCISD+Q (−3.1 kcal/mol),[11,12] XMCQDPT2 (−2.8 kcal/mol),[11] and EOM-CCSD (−3.0 kcal/mol).[12] The negative value indicates that trans-PSB3 is more stable than cis-PSB3 because throughout this work the cis-PSB3 energy is taken as the reference.

Table 2

Energy Differences ΔE (in kcal/mol) between the Singlet Ground State of cis PSB3 Isomer and The Singlet Ground State trans Isomer, the TSCT, and TSDIR obtained by Gozem et al.[11] from CASSCF-Based Calculationsa

method	ref	ΔE trans	ΔE TSCT	ΔE TSDIR	ΔES-T cis
VMC/JSD	this work	–2.9(2)	44.7(2)	51.1(2)	62.9(2)
VMC/JAGP	this work	–2.9(2)	45.2(2)	51.7(2)	66.2(2)
LRDMC/JAGP	this work	–3.0(2)	45.5(2)	51.4(2)	63.9(2)
MRCISD+Q	(11, 12)	–3.1	48.7	54.9
XMCQDPT2	(11)	–2.8	46.9	50.5
EOM-CCSD	(12)	–3.0	46.6	52.5

The energy difference ΔE between the first singlet and the first triplet electronic states of the cis isomer is also reported. The reported QMC results are compared with MRCISD+Q, XMCQDPT2, and EOM-CCSD calculations.

At VMC/JAGP level, the TSCT (45.5(2) kcal/mol) is lower in energy than the TSDIR (51.7 kcal/mol), making the CT path energetically favored, similarly to what was reported by the aforementioned correlated approaches.[11,12] LRDMC and VMC findings are equal within the stochastic error to the VMC/JAGP values, evidence that the trial wave function Ψ is fully optimized. The singlet–triplet gap for the cis-PSB3 is also reported, with a difference of 2.3(3) kcal/mol between VMC and LRDMC using the complete JAGP.

The VMC/JAGP energy profile of the MEPCT path (Figure 3a) is characterized by a shallow plateau around the transition state structure, at variance with the shape of the MEPDIR path (Figure 3c); moreover, expensive LRDMC calculations do not alter the picture. Parts b and d of Figure 3 show the ratio between λLUMO and λHOMO, according to the AGP expansion given in eq 14; as already discussed by some of us in the case of the application of the AGP ansatz on the diradical twisted ethylene C2H4,[32] very small values of this ratio correspond to a single-reference wave function, with the lowest molecular orbitals doubly occupied (the ratio is exactly zero in the limit of a pure single Slater determinant), whereas large values of the ratio indicate two (near)-equivalent configurations. Following the analysis in ref (32), the absolute value of λLUMO/λHOMO approaching to unity means that two configurations are contributing with the same weight to the electronic structure. λHOMO and λLUMO correspond to λ1 and λ2 in the simple (2,2) model in eq 18. The wave function along the MEPCT path (Figure 3b) is dominated by a single configuration, with the λLUMO/λHOMO oscillating around the value for the cis isomer (−0.0635), as expected by the previous investigations on the MEPCT path of PSB3. Similar arguments can be found in ref (67). The same behavior has been found along the MEPDIR path (Figure 3b), with no appreciable contribution given by higher-energy configurations; this result unequivocally shows that the introduction by QMC of balanced dynamical correlation strongly alters the CASSCF description of the electronic structure of PSB3 along the MEPDIR path, similarly to what found in MRCISD+Q calculations.

Figure 3

Left column: (a) thermal isomerization energy profile (with respect to the cis PSB3) and (b) ratio between λLUMO and λHOMO (see eq 14) along the MEPCT reaction coordinate. Right column: (c) thermal isomerization energy profile and (d) ratio between λLUMO and λHOMO along the MEPDIR reaction coordinate. SD stands for single determinant, e.g., a single configuration (with ratio λLUMO/λHOMO exactly equal to zero). Error bars are within the symbols. For both paths CASSCF and MRCISD+Q energy profiles are taken from refs (11,12).

The AGP wave function can be formally expanded into a linear combination of Slater determinants (eq 14), so the single-electron molecular orbitals are obtained by the diagonalization of the geminal coupling matrix (eq 10) and are defined within a correlated framework. The AGP spans the seniority number Ω = 0 sector in Hilbert space, with double occupation for each orbital. The terms “closed-shell” and “open-shell” are widely used to indicate systems without and with unpaired electrons, respectively, implying that the molecular orbitals come from single-reference approaches, like Hartree–Fock or DFT. In the case of highly twisted configurations of PSB3, for instance, the π2 state, involving charge transfer with respect to the equilibrium ground state, has closed-shell character, while the diradical ππ* state has open-shell character: for ππ* static correlation plays an important role, and a multiconfigurational approach must be used. Such definitions strictly depend on the choice of the reference for the molecular orbitals, as explained in section 2.3. Even though the AGP wave function is formally characterized by only doubly occupied orbitals, its application is not limited so far to the study of closed-shell systems because the molecular orbitals involved in the AGP expansion are the results of a variational optimization and they can be qualitatively different from the Hartree–Fock ones, as explicitly shown in section 2.3. Data show that the MEPDIR path does not have a diradical character anymore; the VMC/JAGP description of MEPCT and MEPDIR paths produces a ground state of closed-shell π2 character, similarly to the MRCISD+Q.[13]

The results above are again reflected in the charge transfer profiles along both MEPCT and MEPDIR paths shown in Figure 4. Indeed, it is clear that in both MEPCT and MEPDIR, the wave function gains charge transfer character, with TSCT and TSDIR both having ca. 65% of the positive charge on the allyl moiety. The picture presented at the VMC/JAGP and LRDMC/JAGP levels is again at variance with what is found at the CASSCF level of theory. With CASSCF, MEPDIR passes through a transition state, TSDIR, which has diradical character and has almost no charge transfer at all, with all the charge localized on the nitrogen-containing moiety. This difference is due to the change in CI position on the energy surface after introducing a balanced description of electronic correlation in the calculations. At the CASSCF level of theory, the CI is peaked and lies in between the TSCT and TSDIR transition states, causing them to have wave functions with different character. At correlated levels of theory, the CI is shifted to a larger BLA value than that corresponding to TSDIR. As a result, both TSCT and TSDIR lie on the same side of the CI (which is now intermediate/sloped) and therefore both have the same wave function character (charge transfer). This change in local topology causes only one of the transition states to remain as a saddle point, TSCT, while TSDIR is no longer a transition state on the S0 potential energy surface. The charge transfer character in both MEPCT and MEPDIR decreases as PSB3 moves toward the cis or trans isomers of PSB3, converging to the values corresponding to the two minima.

Figure 4

Charge transfer along the MEPCT and MEPDIR paths. CASSCF and MRCISD charge-transfer profiles are taken from refs (11,12).

Furthermore, the MEPCT maximum (i.e., the TSCT) region is flatter for the two QMC levels than for MRCISD+Q level; this may be due to a better treatment of the electronic dynamical correlation in this region.

BLA Path

In the CASSCF(6,6)/6-31G* landscape, the CI is located between the TSCT and TSDIR. Consistently with the other approaches, PSB3 assumes a closed-shell, charge transfer character (π2) at smaller BLA values than the CI, while its wave function becomes covalent and diradical (ππ*) at larger BLA values (upper panel of Figure 5). As anticipated above, the crossing between the two states does not produce a peaked CI, like in the CASSCF case, but instead a sloped and intermediate CI. Even though we are not able to identify the exact position of the CI, VMC/JAGP calculations clearly show that CI moves toward large values of BLA (∼0.075 Å), beyond the TSDIR. The present result differs from the collection of data obtained by other correlated approaches, locating the CI around 0.03–0.04 Å,[11−13,15] with the exception of the QD-NEVPT2/CAS(6,6) analysis (extrapolated value at 0.05 Å).[11] When using the LRDMC/JAGP approach, the CI position (∼0.06 Å) comes closer to MRCISD+Q data. We observe that the charge transfer (π2) curves obtained using VMC/JAGP and LRDMC/JAGP are actually more stable, with respect to cis-PSB3 energy, than that obtained with MRCISD+Q, consistently with the lower energy TSCT of VMC/JAGP and LRDMC/JAGP. Meanwhile, the diradical (ππ*) curve from VMC/JAGP or LRDMC/JAGP are less stable that those from MRCISD+Q. This is what causes the CI from VMC/JAGP and LRDMC/JAGP to shift to higher BLA than in MRCISD+Q. Moreover, it is important to point out that relaxing the structures and the minimum energy paths of PSB3 at the QMC correlated level may produce slightly different topologies for the BLA and MEP scans close to the CI.

Figure 5

Energy profile along the BLA path (with respect to the cis PSB3) and the related λLUMO/λHOMO ratio for the π2 and ππ* configurations, calculated by VMC and LRDMC methods on the variationally optimized JAGP wave functions. The triplet energy profile is also reported. The MRCISD+Q, MRCISD, and CASSCF profiles are shown, for comparison, and the conical intersections (CI) obtained with the different approaches are marked in the plot.

In the lower panel of Figure 5, the VMC/JAGP points corresponding to the ππ* state are characterized by a high ratio (in absolute value) between λHOMO and λLUMO, close to 1. As already mentioned, this is due to the multiconfigurational nature of the diradical state, where the frontier orbitals are near-degenerate and singly occupied. From the AGP analysis, near-degeneracy for S0 is therefore found in the large-BLA portion of the path. This reinforces the results of the MEPCT and MEPDIR paths above because it is clear here that TSCT and TSDIR both lie on the same side of the CI after introduction of the dynamical electronic correlation. Also, one can see here that while TSCT remains a saddle point on the S0 potential energy surface, TSDIR actually becomes a minimum on the S1 surface and is no longer a transition state as in CASSCF. The π2 state, also extending in the diradical portion of space previously defined by the CASSCF study, is clearly single-reference, and it dominates the wave function of PSB3 up to a BLA of 0.05 Å.

The reliability of the present results is further confirmed looking at the triplet energy profile along the BLA path (blue triangles in the upper panel of Figure 5, filled and open symbols for VMC and LRDMC, respectively). Because the comparison of the π2 and ππ* energies along the BLA scan is a fundamental step in order to give an accurate representation of the ground state potential energy surface of PSB3 surrounding the CI, the triplet energy profile can be considered a lower bound for the ππ* energy. For diradicals, triplet should be the ground state spin multiplicity, according to the molecular version of Hund’s rule:[63,64] this is the case, for instance, of the orthogonally twisted ethylene molecule[32] where, due to the homolytic cleavage of the double bond, the wave function is dominated by two configurations with the same weight and the two involved p atomic orbitals have zero overlap. The structures of the BLA path with diradical character are characterized by a torsion of about 90 deg around the formal central double bond of PSB3, similar to the prototypical example given by C2H4 system: the central double bond is broken and the two p orbitals are mutually (almost) perpendicular. This manifest similarity between the twisted diradical PSB3 and the prototypical C2H4 system yields to reasonably expect that the ground state of the twisted diradical PSB3 is also a triplet, as we observe a posteriori with VMC/JAGP and LRDMC/JAGP calculations. [However, it is not true in general for all diradicals that the triplet state has a a lower energy that the singlet state, for example, disjointed and non-Kekule molecules, like the tetramethyleneethane,[31,68] have a very small singlet–triplet gap with the singlet lower in energy in some specific geometries, while the ground state of oligacenes larger than hexacene[69] is also a singlet.]

As one can see from the upper panel of Figure 5, the triplet energy is higher than the π2 singlet energy for BLA values smaller than 0.055 Å for VMC and 0.05 Å for LRDMC, and consequently the same certainly occurs for the ππ* state. A similar gap of ∼10–12 kcal/mol is also found at CASSCF(6,6)/6-31G* level in the diradical portion of the BLA path.

These simple arguments, combined with the fact that the convergence of the wave function optimization for a triplet state is easier to achieve in the VMC/AGP framework, allow us to be extremely confident with the robustness of the present results for the location of CI along the BLA path.

The ground state wave function is dominated by the π2 state before the CI (i.e., at lower BLA values), while it assumes ππ* character after the crossing; it is interesting to understand how the charge distribution changes in the two configurations. Figure 6 reports, as an example, the charge density difference between the π2 and ππ* states for VMC/JAGP and LRDMC/JAGP CI points along the BLA path, more precisely, for the closest structures to CI found at VMC and LRDMC level. The shape of the isosurfaces obviously resembles the p atomic orbitals involved in the πelectronic structure of the twisted PSB3. The most evident changes in the charged distribution are observed for the two central carbon atoms, where the double bond has been broken because of the torsion: the yellow part of the plot (online color version) corresponds to an accumulation of electronic charge in the π2 state with respect to the ππ*, whereas the green isosurfaces indicate the opposite situation.

Figure 6

Graphical representation of the charge density difference between the π2 and the ππ* configurations for VMC/JAGP and LRMDC/JAGP CI points: the yellow (green) isosurface indicates an excess (reduction) of 0.01 in terms of electronic charge in the π2 state with respect to ππ*.

Such analysis is consistent with the charge transfer nature of π2 configuration and with the covalent and diradical character of ππ*. We know from Figure 7 that the charge transfer of π2 along the path oscillates between 0.63 and 0.72, and that the charge transfer of ππ* is much smaller (∼0.10). A large (small) value of charge transfer means an excess (reduction) of electronic charge on the protonated imine heteroallyl moiety, according to the plot reported in Figure 7 for VMC and LRDMC calculations, with the charge transfer of the triplet coinciding with the ππ* values.

Figure 7

Charge transfer along the BLA path for the π2, ππ*, and triplet configurations, obtained by VMC and LRDMC approaches (stochastic errors are smaller than point size). The MRCISD and CASSCF profiles are shown for comparison.

Conclusions

VMC and LRDMC methods have been used to compute electronic and energetic properties of PSB3, using the JAGP ansatz.

VMC/JAGP and LRDMC/JAGP calculations on three different CASSCF paths reveal the fundamental role played by a balanced description of the dynamical correlation for the correct representation of the ground state energy surface in the proximity of the conical intersection. The VMC and LRDMC calculations significantly alter the CASSCF landscape, inverting the relative stability of the MEPCT and MEPDIR paths, similarly to what obtained by other correlated approaches.[11−13,15] However, the region surrounding the TSCT appears more flat when computed at the QMC levels.

On both MEPCT and MEPDIR paths, the PSB3 wave function assumes a charge transfer and single-reference character: the MEPDIR path is therefore not diradical anymore when analyzed at correlated level. This finding is immediately confirmed by the fact the CI is pushed toward values of BLA larger than that corresponding to the TSDIR (∼0.075 Å with VMC/JAGP and ∼0.06 Å with LRDMC/JAGP, while the reference calculations give a CI located at around 0.03–0.04 Å).

Furthermore, in the present study, the multiconfigurational nature of AGP has been explained in detail, following the work reported in ref (32), for the description of the diradical ππ* configuration of PSB3. Even though the AGP is defined in the seniority number Ω = 0 subsector of the Hilbert space (e.g., the set of molecular orbitals is doubly occupied), a unitary transformation of the AGP-optimized molecular orbitals allows one to always map the ππ* state of highly twisted PSB3 geometries, with two unpaired electrons, into the AGP Ω = 0 subsector.

The role of the dynamical electronic correlation has been found to be essential in order to get a reliable description of the ground state of PSB3 around the CI. Therefore, mechanistic or dynamics studies using methods which do not incorporate these effects need to be performed and interpreted critically. The recent improvements in forces calculations using QMC approaches[70] suggest that there will be soon the possibility to compute MEPCT and MEPDIR minimum energy paths based on the VMC/JAGP method, yielding to a more consistent comparison with the energies and geometries obtained from CASSCF in order to further clarify the main features of PSB3.