TAO-DFT with the Polarizable Continuum Model

For the ground-state properties of gas-phase nanomolecules with multi-reference character, thermally assisted occupation (TAO) density functional theory (DFT) has recently been found to outperform the widely used Kohn–Sham DFT when traditional exchange-correlation energy functionals are employed. Aiming to explore solvation effects on the ground-state properties of nanomolecules with multi-reference character at a minimal computational cost, we combined TAO-DFT with the PCM (polarizable continuum model). In order to show its usefulness, TAO-DFT-based PCM (TAO-PCM) was used to predict the electronic properties of linear acenes in three different solvents (toluene, chlorobenzene, and water). According to TAO-PCM, in the presence of these solvents, the smaller acenes should have nonradical character, and the larger ones should have increasing polyradical character, revealing striking similarities to the past findings in the gas phase.


Introduction
Electronic structure methods are commonly used to compute the properties of gasphase molecules. Nonetheless, the properties of molecules in solutions can be very different from those of gas-phase molecules. Therefore, it is important to realize the effects of solvents on the properties of solute molecules, especially when polar solvents are involved. While it is possible to adopt explicit solvent models (wherein solvation effects are taken into consideration by explicitly including the molecular details of each solvent molecule) [1][2][3], the resulting electronic structure calculations can, however, be computationally intractable in many applications.
In order to circumvent this difficulty, several implicit solvent models have been devised to study the properties of molecules in polarizable solvents over the past few decades. Among them, the Kirkwood-Onsager model [4][5][6] is perhaps the simplest one. However, the Kirkwood-Onsager (KO) model is not appropriate for non-spherical solute molecules because of its underlying assumptions.
In TAO-PCM, the solute is treated in a quantum-mechanical manner with TAO-DFT [37], and the solvation effect is modeled implicitly with the PCM. Specifically, in the presence of the continuum solvent, the TAO-PCM solute free energy G TAO-PCM [ρ] is a functional of the solute electron density ρ(r) (atomic units were adopted throughout this work): where E TAO [ρ] (e.g., see Equation (27) of Ref. [37]) is the gas-phase solute electronic energy in TAO-DFT, V NN is the solute nuclear repulsion energy, and G PCM [ρ] (e.g., see Equation (14) of Ref. [21]) is the solute-solvent interaction free energy in the PCM. In order to obtain the solute electron density ρ(r), the TAO-PCM self-consistent equations are given by Here, v TAO (r) (e.g., see Equation (18) of Ref. [37]) is the gas-phase solute effective potential in TAO-DFT, v PCM (r) = δG PCM [ρ]/δρ(r) (e.g., see Equation (17) of Ref. [21]) is the reaction potential in the PCM, and the solute electron density ρ(r) is expressed as where f i (the occupation number associated with the i-th TAO-orbital ψ i (r)) is generated by the FD distribution (with some fictitious temperature θ): satisfying the constraints 0 ≤ f i ≤ 1 and ∑ i f i = N, with i being the energy associated with the i-th TAO-orbital ψ i (r), and µ being the chemical potential.
For the θ = 0 case, as TAO-DFT reduces to KS-DFT, TAO-PCM reduces to KS-PCM. Note also that the PCM component G PCM [ρ], as well as v PCM (r), which depends only on the solvent dielectric constant and the solute (electronic and nuclear) charge density [21], has the same expression in TAO-PCM as in KS-PCM. In addition, for an isolated solute molecule (i.e., = 1), as G PCM [ρ] = 0 and v PCM (r) = 0, TAO-PCM reduces to TAO-DFT, and KS-PCM reduces to KS-DFT.

Computational Details
We performed all calculations using the 6-31G(d) basis set with the software package of Q-Chem 4.4 [69]. As mentioned previously, in TAO-PCM, TAO-DFT is adopted for the electronic structure calculations, and the PCM is adopted for the modeling of solvent effects. In this study, for the TAO-DFT part, we adopted TAO-LDA [37], which is TAO-DFT with the local density approximation (LDA) XC and θ-dependent energy functionals (with θ = 7 mhartree (i.e., the suggested fictitious temperature)). Specifically, we used the PW92 parametrization [70] of the LDA correlation energy functional and the LDA θ-dependent energy functional [37] obtained with the Perrot parametrization [71] of the LDA noninteracting kinetic free energy functional. While other θ values are also possible for TAO-LDA (which could lead to slightly different electronic properties) [37,40,44], TAO-LDA (θ = 7 mhartree) has recently been shown to yield a qualitatively similar MR character of various alternant polycyclic aromatic hydrocarbons (e.g., n-acenes) [46] when compared with an accurate MR electronic structure method [28]. In addition, TAO-LDA (θ = 7 mhartree) has been commonly employed to study the electronic properties of gas-phase nanomolecules with MR character [45,[48][49][50][51][52][53]. Because of its decent compromise between accuracy and efficiency, in this TAO-PCM study, we adopted TAO-LDA (θ = 7 mhartree) for the TAO-DFT part.
On the other hand, for the PCM part, we adopted C-PCM [11], which has been a popular PCM due to its simplicity in formalism and implementation as well as its decent balance between efficiency and accuracy. Note also that C-PCM is an approximation of COSMO [8] and can also be regarded as the high-dielectric limit associated with IEF-PCM [12][13][14][15]. Regarding the construction and discretization of a solute cavity surface, the solute cavity surface [3] was formed using a union of atom-centered spheres with their radii being 1.2 times the atomic vdW radii [72][73][74]. In addition, the SWIG (switching/Gaussian) approach was adopted for smooth cavity discretization [75,76], with 194 Lebedev grid points per atomic sphere. Furthermore, for the solvent dielectric constant , we adopted the Q-Chem default value, i.e., 2.38 for toluene, 5.62 for chlorobenzene, and 78.39 for water.
For comparison with the results of TAO-PCM (i.e., TAO-LDA/C-PCM), we also reveal some results obtained with the corresponding KS-PCM (i.e., KS-LDA/C-PCM), wherein KS-LDA (i.e., KS-DFT with the LDA XC energy functional, which is the same as TAO-LDA with θ = 0) was used for the electronic structure calculations, and C-PCM was used for the modeling of solvent effects.

Singlet-Triplet Energy Gap
A singlet-triplet energy gap (ST gap) can provide many insights on chemical processes. In order to obtain the ST gap of n-acene (i.e., the solute) in the continuum solvent, we fully optimized the geometries of the lowest singlet and triplet states of n-acene in the continuum solvent using spin-unrestricted TAO-PCM (i.e., TAO-LDA/C-PCM), and computed the ST gap of n-acene in the continuum solvent using where G US and G UT are the spin-unrestricted TAO-PCM free energies (given by Equation (1)) associated with the lowest singlet and lowest triplet states, respectively, of n-acene in the continuum solvent. For comparison purposes, we additionally report the results of the corresponding KS-PCM (i.e., KS-LDA/C-PCM). As presented in Figures 2 and 3 (see Tables S1 and S2 in the SI (Supplementary Information) for additional information), the ST gaps of n-acenes in three different solvents (toluene, chlorobenzene, and water) are essentially the same as those of n-acenes in the gas phase, indicating that the solvation effects play insignificant roles here. In addition, for all the cases investigated, n-acenes in different media possess singlet ground states. As molecules with MR character are commonly characterized by small ST gaps, similar to the gas-phase situations [37][38][39][43][44][45], in the three solvents, the larger n-acenes, which have smaller ST gaps, are expected to have more significant MR character than the smaller n-acenes.   On the other hand, the ST gap of TAO-PCM decreases with the size of n-acene in a monotonic manner, but the ST gap of KS-PCM exhibits an unexpected increase beyond 10-acene. Therefore, we also explored the potential causes of such discrepancies by examining the expectation values ofŜ 2 (i.e., the total spin-squared operator) associated with the lowest singlet and lowest triplet states of n-acene in the continuum solvent, calculated by spin-unrestricted KS-PCM (see Tables 1 and 2). For the KS-PCM wavefunction with an artificial mixing of several spin-states (i.e., the so-called spin contamination), the respective Ŝ 2 value will be larger than the exact Ŝ 2 value (i.e., 2 for the lowest triplet state and 0 for the lowest singlet state) [49,53,77]. Accordingly, the difference between the Ŝ 2 value of KS-PCM and the exact Ŝ 2 value can be regarded as a measure of the degree of spin contamination in the spin-unrestricted KS-PCM wavefunction. As shown, for the smaller n-acenes, the Ŝ 2 values of KS-PCM are very close to the exact Ŝ 2 values, implying that the lowest singlet and triplet states of the smaller n-acenes in different media, obtained from spin-unrestricted KS-PCM, essentially have no spin contamination. Nonetheless, for the lowest singlet states of some larger n-acenes (e.g., n = 9-12), the Ŝ 2 values of KS-PCM are considerably larger than the exact Ŝ 2 values, indicating that the lowest singlet states of some larger n-acenes in different media, obtained from spin-unrestricted KS-PCM, are heavily spin-contaminated. Consequently, the larger n-acenes can have MR character in the lowest singlet states, and the unexpected increase in the ST gap of KS-PCM beyond 10-acene could simply be an artifact related to spin contamination. Table 1. Ŝ 2 for the lowest singlet state of n-acene in three different solvents (toluene, chlorobenzene, and water) and in the gas phase, calculated using spin-unrestricted KS-PCM (i.e., KS-LDA/C-PCM). For the lowest singlet state, the exact Ŝ 2 is 0. Furthermore, we can also assess the effects of spin contamination in the lowest singlet free energies of n-acenes in different media, obtained from spin-unrestricted TAO-PCM and KS-PCM. Because of the spin-symmetry constraint [37][38][39], for the PCM with an exact electronic structure method, the lowest singlet free energy of n-acene in the continuum solvent, calculated by the spin-restricted formalism, has to be identical to that calculated by the spin-unrestricted formalism. In order to assess whether this constraint can be obeyed in TAO-PCM, we fully optimized the geometries of the lowest singlet state of n-acene in the continuum solvent using spin-restricted and spin-unrestricted TAO-PCM, and subsequently computed the E UR value (i.e., the difference between the lowest spinunrestricted and lowest spin-restricted singlet free energies) of n-acene in the continuum solvent using where G RS and G US are the TAO-PCM free energies (given by Equation (1)) associated with the lowest singlet state of n-acene in the continuum solvent, computed using the spin-restricted and spin-unrestricted formalisms, respectively. For comparison, we also report the corresponding KS-PCM results. As shown in Table 3, the E UR values of nacenes in the continuum solvent, calculated by KS-PCM, are rather large for some larger n-acenes (e.g., E UR ≥ 3.20 kcal/mol for 10-acene), yielding unphysical spin-symmetrybreaking results in the corresponding spin-unrestricted KS-PCM calculations. This implies that, in spin-unrestricted KS-PCM, the up-spin electron density can be rather different from the down-spin electron density for the lowest singlet state of larger n-acene in the continuum solvent (even when these spin densities are required to be the same due to the spin-symmetry constraint). Consequently, for the lowest singlet state of larger n-acene in the continuum solvent, an electronic property (which depends on the spin densities) obtained with spin-unrestricted KS-PCM can be rather different from that obtained with spin-restricted KS-PCM (even when these electronic properties are required to be the same due to the spin-symmetry constraint). Clearly, such an unphysical spin-symmetry-breaking feature of spin-unrestricted KS-PCM is very undesirable. By contrast, the E UR values of n-acenes in the continuum solvent, calculated by TAO-PCM, are essentially zero for all the cases investigated, leading to essentially no unphysical spin-symmetry-breaking results in the corresponding spin-unrestricted TAO-PCM calculations. Table 3. Difference E UR (in kcal/mol) between the lowest spin-unrestricted and lowest spin-restricted singlet free energies of n-acene in the gas phase and in three different solvents (toluene, chlorobenzene, and water), calculated using KS-PCM (i.e., KS-LDA/C-PCM). In summary, the ST gaps of the larger n-acenes in the continuum solvent, computed using KS-PCM, can be severely influenced by spin contamination, seriously degrading the accuracy of KS-PCM in predicting the ST gaps and possibly other electronic properties (e.g., those that depend on the spin densities). As will be shown later, the lowest singlet states (i.e., ground states) of the larger n-acenes in different media exhibit an increasing polyradical nature, wherein KS-PCM could lead to unreliable electronic properties. Consequently, we merely present the TAO-PCM results hereafter.

Vertical Electron Affinity/Ionization Potential and Fundamental Gap
For a neutral solute molecule (with N electrons) in a continuum solvent, the free energy gained when an electron is added to the neutral solute molecule (without altering the solute geometry) is defined as the vertical electron affinity the free energy that is required to remove an electron from the neutral solute molecule (without altering the solute geometry) is defined as the vertical ionization potential and their difference is defined as the fundamental gap [38,39] In this study, G N , G N+1 , and G N−1 are the spin-unrestricted TAO-PCM free energies (given by Equation (1)) associated with the neutral, anionic, and cationic states, respectively, of n-acene (i.e., the solute) in the continuum solvent. Figure 4 (also see Tables S3-S5 in the SI for additional information) presents the IP v , EA v , and E g of ground-state n-acene in three different solvents (toluene, chlorobenzene, and water) and in the gas phase, obtained with spin-unrestricted TAO-PCM (i.e., TAO-LDA/C-PCM). As shown, for each solvent studied, the IP v decreases in a monotonic manner as the size of n-acene increases. In addition, for each n-acene studied, the IP v decreases as the solvent dielectric constant increases, indicating that less free energy is needed to remove an electron from n-acene in the continuum solvent of larger . Therefore, the IP v of n-acene is the largest in the gas phase ( = 1), and is the smallest in water ( = 78.39).
A similar but opposite trend is observed from the EA v . As the size of n-acene increases, the EA v increases in a monotonic manner for each solvent studied. In addition, for each n-acene studied, the EA v increases as the solvent dielectric constant increases, indicating that more free energy is obtained when an electron is added to n-acene in the continuum solvent of larger . Accordingly, the EA v of n-acene is the smallest in the gas phase ( = 1) and is the largest in water ( = 78.39).
Because of the monotonically decreasing nature of IP v and the monotonically increasing nature of EA v with an increasing n-acene size, unsurprisingly, for each solvent studied, the E g decreases in a monotonic manner with an increasing n. Moreover, for each n-acene studied, the E g decreases as the solvent dielectric constant increases. Consequently, the E g of n-acene is the largest in the gas phase ( = 1) and is the smallest in water ( = 78.39).
In short, our study suggests that the solvation effects are rather important in the IP v , EA v , and E g values of ground-state n-acenes since their values can be easily tuned by changing solvents. This highlights the significance of TAO-PCM in the modeling of solvent effects on the GS properties of nanomolecules with MR character.

Symmetrized von Neumann Entropy
Similar to the TOONs in TAO-DFT [37][38][39]44], the TOONs in TAO-PCM can be approximately viewed as the NOONs (natural orbital occupation numbers) [78] for the ground state of a solute molecule in the continuum solvent. Hence, the strength of the MR character of ground-state n-acene in the continuum solvent can be estimated using the corresponding value of S vN (symmetrized von Neumann entropy) [38,39,79]: where { f i } are the TOONs (see Equation (4)) of ground-state n-acene in the continuum solvent, obtained from spin-unrestricted TAO-PCM (i.e., TAO-LDA/C-PCM). Note that the S vN value is close to 0 for an electronic system with single-reference character (i.e., all the TOONs are very close to either 0 or 1), and can be much larger than 0 for an electronic system with significant MR character (i.e., some of the TOONs considerably deviate from the values of 0 and 1). As shown in Figure 5 (also see Table S6 in the SI for additional information), the S vN value of ground-state n-acene in the continuum solvent increases with n, implying that the respective strength of MR character should also increase with n. In addition, the S vN values of ground-state n-acenes in three different solvents (toluene, chlorobenzene, and water) are essentially identical to those of ground-state n-acenes in the gas phase [38,39,43,45], suggesting that the solvation effects are unimportant in affecting the S vN values (i.e., a measure of the strength of MR character) of ground-state n-acenes. Symmetrized von Neumann Entropy gas phase @toluene @chlorobenzene @water Figure 5. Symmetrized von Neumann entropy for the ground state of n-acene in three different solvents (toluene, chlorobenzene, and water) and in the gas phase, obtained from spin-unrestricted TAO-PCM (i.e., TAO-LDA/C-PCM).

Active TAO-Orbital Occupation Numbers
In order to see why the aforementioned S vN value increases with n, we reveal the active TOONs of ground-state n-acene (with N electrons) in the continuum solvent, obtained with spin-restricted TAO-PCM (i.e., TAO-LDA/C-PCM). The (N/2) th /(N/2 + 1) th TAO-orbital is denoted as the HOMO (highest occupied molecular orbital)/LUMO (lowest unoccupied molecular orbital) as convention. In addition, the TAO-orbitals with an occupation number in the range of 0.2-1.8 are considered as the active TAO-orbitals.
As plotted in Figure 6, the active TOONs of ground-state n-acene in the continuum solvent possess interesting patterns. First, the smaller n-acenes (e.g., n < 6) exhibit a nonradical nature since all the TOONs are in the vicinity of either 0 or 2. Second, as the size of nacene increases, there are more active TAO-orbitals and/or the active TOONs are closer to 1, displaying that the larger ground-state n-acenes in the continuum solvent exhibit an increasing polyradical nature. Therefore, there is a transition from a nonradical nature to polyradical nature of ground-state n-acene in the continuum solvent, causing an increase in S vN with n. In addition, the active TOONs of ground-state n-acenes in three solvents (toluene, chlorobenzene, and water) are essentially the same as those of ground-state n-acenes in the gas phase (i.e., including the previously observed curve-crossing behavior [37,39], which was recently confirmed by an accurate MR method [31]) [37,39,[43][44][45], indicating that the solvation effects are indeed insignificant in changing the MR nature of ground-state n-acenes (i.e., showing consistency with the S vN analysis).

Active TAO-Orbitals in Real Space
In real space, we reveal the active TAO-orbitals (see Figures S1-S3 in the SI), such as the HOMO and LUMO, for the ground states of some representative n-acenes (e.g., 4-acene, 6-acene, and 8-acene) in the continuum solvent, obtained with spin-restricted TAO-PCM (i.e., TAO-LDA/C-PCM). As n increases, the active TAO-orbitals exhibit an increasing trend to localize near the edges of n-acene. In addition, the active TAO-orbitals of ground-state n-acenes in three solvents (toluene, chlorobenzene, and water) and in the gas phase [45] look very much alike, suggesting that the solvation effects play unimportant roles here.

Conclusions
In conclusion, we developed TAO-PCM, combining TAO-DFT with the PCM, to study solvation effects on the GS properties of nanomolecules with MR character at a minimal computational cost. In addition, we adopted TAO-PCM (i.e., TAO-LDA/C-PCM) to explore the electronic properties (the ST gap, fundamental gap, vertical electron affinity/ionization potential, etc.) of n-acene (i.e., the solute) with n = 2-20 in three different solvents (toluene, chlorobenzene, and water) and in the gas phase. For comparison, our TAO-PCM results (e.g., the ST gap and E UR value of n-acene) were also compared with the corresponding KS-PCM (i.e., KS-LDA/C-PCM) results.
In the three solvents and in the gas phase, since the larger n-acenes have been found to reveal MR character in the lowest singlet states (i.e., ground states), KS-PCM can yield incorrect electronic properties (e.g., those that depend on the spin densities). For example, the ground states of some larger n-acenes in different media, obtained from spin-unrestricted KS-PCM, have been shown to be highly spin-contaminated, yielding unphysical spinsymmetry-breaking effects. While accurate MR methods may remedy this issue, these MR methods and related PCMs could be computationally impossible for the larger n-acenes in the gas phase and solution phase, respectively. Hence, it is well justified to use TAO-PCM to explore the electronic properties of n-acenes in different media because of its decent compromise between accuracy and efficiency.
According to TAO-PCM, in the three solvents, the smaller n-acenes (e.g., n < 6) reveal a nonradical nature, and the larger n-acenes reveal an increasing polyradical nature, similar to the previous findings in the gas-phase situations [37][38][39][43][44][45]. In addition, significant changes in some of the electronic properties (e.g., the fundamental gap and vertical electron affinity/ionization potential) of ground-state n-acene have been found due to the presence of these solvents, highlighting the importance of TAO-PCM in the modeling of solvent effects on the GS properties of nanomolecules with MR character. Among all the electronic properties studied, the solvent polarity does not have a significant impact when the net charge is not varied in the process (e.g., the singlet-triplet gap, symmetrized von Neumann entropy, active TOONs, etc.), and has a significant impact when the net charge is varied in the process (e.g., the vertical electron affinity/ionization potential).
While TAO-PCM seems to be promising for the study of solvation effects on the GS properties of nanomolecules with MR character, TAO-PCM may not work well in situations where the underlying PCM (or other implicit solvation models) could fail badly. Therefore, in the near future, we plan to explore the possibility of combining TAO-DFT with an efficient explicit solvent model to resolve this issue.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/nano13101593/s1, Figure S1: Real-space representation of active TAO-orbitals (HOMO and LUMO) for the ground state of 4-acene in (a) the gas phase and in three different solvents: (b) toluene, (c) chlorobenzene, and (d) water, calculated using spin-restricted TAO-PCM (i.e., TAO-LDA / C-PCM), at an isovalue of 0.02 e/Å 3 ; Figure S2: Real-space representation of active TAO-orbitals (HOMO and LUMO) for the ground state of 6-acene in (a) the gas phase and in three different solvents: (b) toluene, (c) chlorobenzene, and (d) water, calculated using spinrestricted TAO-PCM (i.e., TAO-LDA / C-PCM), at an isovalue of 0.02 e/Å 3 ; Figure S3: Real-space representation of active TAO-orbitals (HOMO and LUMO) for the ground state of 8-acene in (a) the gas phase and in three different solvents: (b) toluene, (c) chlorobenzene, and (d) water, calculated using spin-restricted TAO-PCM (i.e., TAO-LDA / C-PCM), at an isovalue of 0.02 e/Å 3 ; Table S1: Singlet-triplet energy gap E ST (in kcal/mol) of n-acene in the gas phase and in three different solvents (toluene, chlorobenzene, and water), calculated using spin-unrestricted TAO-PCM (i.e., TAO-LDA / C-PCM); Table S2: Singlet-triplet energy gap E ST (in kcal/mol) of n-acene in the gas phase and in three different solvents (toluene, chlorobenzene, and water), calculated using spin-unrestricted KS-PCM (i.e., KS-LDA / C-PCM); Table S3: Vertical ionization potential IP v (in eV) for the ground state of n-acene in the gas phase and in three different solvents (toluene, chlorobenzene, and water), calculated using spin-unrestricted TAO-PCM (i.e., TAO-LDA / C-PCM); Table S4: Vertical electron affinity EA v (in eV) for the ground state of n-acene in the gas phase and in three different solvents (toluene, chlorobenzene, and water), calculated using spin-unrestricted TAO-PCM (i.e., TAO-LDA / C-PCM); Table S5: Fundamental gap E g (in eV) for the ground state of n-acene in the gas phase and in three different solvents (toluene, chlorobenzene, and water), calculated using spin-unrestricted TAO-PCM (i.e., TAO-LDA / C-PCM); Table S6: Symmetrized von Neumann entropy S vN for the ground state of n-acene in the gas phase and in three different solvents (toluene, chlorobenzene, and water), calculated using spin-unrestricted TAO-PCM (i.e., TAO-LDA / C-PCM).

Data Availability Statement:
The numerical data supporting the findings of the present work are available from the authors upon appropriate request.