Absorption Properties of Large Complex Molecular Systems: The DFTB/Fluctuating Charge Approach.

We report on the first formulation of a novel polarizable QM/MM approach, where the density functional tight binding (DFTB) is coupled to the fluctuating charge (FQ) force field. The resulting method (DFTB/FQ) is then extended to the linear response within the TD-DFTB framework and challenged to study absorption spectra of large condensed-phase systems.

Introduction

The theoretical modeling of large molecular systems, with application in biological and technological fields, is one of the most challenging tasks for theoretical and computational chemistry.[1,2] In fact, the description of large molecular systems requires the treatment of a large number of degrees of freedom, from both nuclear and electronic points of view.[1,2] For this reason, high-level quantum mechanics (QM) methods are usually not applicable because they are usually associated with an unfavorable scaling with the number of atoms (and electrons).[3] Different strategies, usually based on chemical intuition, can be exploited to reduce the dimensionality of the system and to make high-level QM approaches applicable.[4−9] This is, for instance, the case for local excitations, which take place in a specific part of the considered molecule.[5,8,10] However, in many cases, especially for biomolecules, such an approximation may not be chemically justified because the phenomenon is the result of changes in the whole structure.

Semiempirical QM methods have been developed to treat in a realistic way this kind of system, which can be constituted by thousands of atoms.[11−14] Such methodologies introduce a set of integral approximations and parametrizations that make the computation particularly cheap. Clearly, the accuracy of each approach strongly depends on the quality of the parametrization. Among semiempirical methods, one of the most used is the density functional tight binding (DFTB) approach.[15−17] The theoretical starting point of such method is the density functional theory (DFT) energy in the Kohn–Sham (KS) framework, expressed by means of a linear combination of atomic orbitals (LCAO) over a minimal basis set. This quantity is then approximated by means of a Taylor expansion with respect to a reference density truncated at different orders by generating a hierarchy of DFTB methods.[15] In particular, the self-consistent charge DFTB approach (SCC-DFTB), which corresponds to a second-order expression of the KS energy, has been successfully applied to the calculation of energies, geometries, and vibrational frequencies of small organic molecules; its accuracy when compared with experimental values is comparable to that of full DFT calculations performed with a double-ζ plus polarization basis set.[17] Moreover, a time-dependent DFTB (TD-DFTB) approach has been developed to calculate excitation energies in a tight-binding fashion.[18,19] However, it has been shown that the standard pure or hybrid DFT functionals are not able to accurately treat charge-transfer excitations because their extension to the corresponding long-range-corrected versions is necessary.[20,21] In this context, the time-dependent long-range-corrected DFTB approach (TD-LC-DFTB)[22] has recently been proposed and explicitly designed for the treatment of charge-transfer states in large chromophores.

The DFTB approximation allows the boundaries of the systems treatable by most ab initio approaches to be pushed, making possible the QM description of large biomolecules, such as proteins.[19,23] However, most biomolecules are typically dissolved in an external environment, as water is the most common physiological solvent.[24] As for small organic molecules, also in such cases, the external aqueous solution may strongly affect the properties of the biological system.[25] To take into account the solvent effect provided by water, the best compromise between computational cost and accuracy is to resort to the so-called focused models, in which the target system and the environment are described at different levels of theory based on the assumption that the phenomenon is carried on by the target and the environment just perturbs it.[26−28] Among the different focused models, the most accurate are the polarizable QM/molecular mechanics approaches.[29−44] In such methods, the environment molecules are classically and atomistically described by means of a polarizable force field, and mutual solute–solvent polarization is taken into account. In particular, excellent performances have been reported for the QM/fluctuating charge (QM/FQ) in the description of aqueous solutions[45] and recently for different solvents.[46] In such an approach, each solvent atom is endowed with a charge that is adjusted to the external potential generated by the solute density.[28,45,47] Such charges then polarize the QM density by entering the QM Hamiltonian in a mutual polarization fashion. In its basic formulation, the QM/FQ interaction is limited to electrostatics; however, nonelectrostatic interactions can also be considered.[48]

In this work, we have substantially extended the applicability of the polarizable QM/FQ approach by proposing a novel polarizable QM/FQ scheme based on the DFTB approach for the QM portion, allowing for the treatment of large, complex biomolecular systems. To the best of our knowledge, this is the first time that DFTB has been coupled to a polarizable MM approach. The newly developed DFTB/FQ approach has also been extended to the linear response regime by means of the time-dependent DFTB (TD-DFTB) approximation,[18,19] and it has been tested to reproduce the excitation energies of doxorubicin (DOX), an anticancer drug, in aqueous solution and intercalated in DNA and ubiquitin (UBI) protein dissolved in aqueous solution. The Article is organized as follows: In the next section, we briefly recall the DFTB approach, and we formulate the coupling between the DFTB and FQ portions for both ground-state and excitation energies calculations. DFTB/FQ is then applied to the calculation of the excitation energies of DOX, the DOX–DNA complex, and the UBI protein in aqueous solution. Conclusions and perspectives end the Article.

Theory

In this section, the theoretical background of the DFTB/FQ approach is described. To this end, the fundamentals of DFTB and polarizable FQ approaches are briefly recapped, and the formulation of the DFTB/FQ coupling is presented. Then, the extension of the model to the linear response in a TD framework is discussed.

DFT Basis of TB Theory

In the general DFT framework, the energy functional in the KS picture reads[49]where ψ are occupied KS eigenstates, Δ is the Laplacian operator, Vext is the external potential associated with the nuclei–electron interaction, Exc is the exchange-correlation contribution, and VNN is the nuclei–nuclei repulsion term.

In the DFTB theory, the electronic density ρ is expressed as ρ = ρ0 + δρ, where ρ0 is a reference input density and δρ is a fluctuation, which is assumed to be small.[16,19,50] Within this assumption, the exchange-correlation energy contribution can be expanded in a Taylor expansion, and eq becomeswhere we have introduced the expectation value of the zeroth-order Hamiltonian (which depends only on the reference density ρ0) and the so-called repulsive energy contribution Erep. Eγ and EΓ collect the second- and third-order energy terms. Notice that in eqs –6, the usual shorthand notation such that , δρ = δρ(r), , δρ′ = δρ(r′), and so on, is used.[16,19,50]

Different DFTB methods can be defined by truncating the Taylor expansion in eq at different orders. The most basic approach consists of neglecting Eγ and EΓ in eq . This gives rise to a set of non-self-consistent KS equations because the zeroth-order Hamiltonian depends only on the reference density, ρ0. The repulsive energy Erep is approximated as a sum of repulsive, short-ranged, two-body potentials, defined in terms of a set of parameters.[51] The Hamiltonian and overlap S = ⟨ϕμ|ϕν⟩ matrix elements are calculated at a set of relevant interatomic distances and are tabulated. By this, they do not need to be computed for each DFTB calculation, and this results in substantial computational savings as compared with standard DFT. Notice that various parametrizations for Erep and the and S matrix elements have been proposed.[52]

A more sophisticated DFTB method, the SCC-DFTB, can be obtained by retaining Eγ in eq . The density fluctuation δρ is expressed as a sum of localized atomic contributions, δρ = ∑αδρα, which are subsequently approximated through the monopolar term of a multipolar expansion,[18] that is,where Fα(r) is a normalized spherical density fluctuation centered on the αth atom, whereas the net charge Δqα = qα – qα0 is computed through a Mulliken charge analysis. Within such an assumption, Eγ can be rewritten aswhere γ readsTherefore, the total Hamiltonian matrix can be written asΔqξ explicitly depends on MO coefficients through the density matrix, thus introducing a nonlinearity in the Hamiltonian. As a result, DFTB equations must be solved iteratively.

Last, the third-order term EΓ in eq may also be retained, such as in the DFTB3 approach.[53,54]

DFTB/FQ Approach

As stated in the Introduction, in this work, DFTB is coupled to the polarizable FQ force field, which represent each classical atom in terms of a charge q, which is allowed to “fluctuate” so as to fulfill the electronegativity equalization principle, which states that the instantaneous electronegativity χ of each atom must be the same at equilibrium. The total charge on each FQ molecule is fixed to a certain value Q by using Lagrangian multipliers λ. The FQ energy can be written as[47]where qλ is the vector of FQ charges and Lagrange multipliers, C is a vector collecting atomic electronegativities and charge constraints Q, and the M matrix takes into account the interaction kernel between FQ charges and Lagrangian blocks. In particular, the diagonal elements of the FQ–FQ block of M account for the charge self-interaction by means of the chemical hardness η.[45] The minimization of the energy functional in eq leads to a set of linear equations; their solution yields the FQ charges, that isWithin a two-layer QM/MM scheme, the total energy of the DFTB/FQ system is written aswhere EDFTB and EFQ represent the energies of the DFTB and FQ portions and EDFTB/FQ is the interaction energy between the two layers. Here, similarly to most QM/MM approaches, a purely classical interaction term is considered; that is, the DFTB and FQ portions interact through the electrostatic potential generated on the FQ charges by the total DFTB density, that is, the reference density ρ0 and the density fluctuation δρ. Within the DFTB framework, the QM/FQ interaction can be approximated by only taking into account the potential generated by δρ, similar to alternative DFTB/classical couplings.[55−57] Therefore, the corresponding approximated molecular electrostatic potential at the ith FQ charge placed at r can be written aswhere the implicit dependence of the electric potential on the density matrix through Mulliken charges is highlighted. Notice that in eq , the integration over the normalized spherical density fluctuations Fα(r) should be included. (See eq .) However, because the distance between FQ charges and QM atoms is typically larger than any intramolecular distance, we can safely assume that Fα(r) = δ(r – Rα). Therefore, density fluctuations can be described through a set of Mulliken point charges.

Moving back to the total DFTB/FQ energy functional, it can be rewritten asMinimization of the energy functional in eq with respect to both charges and Lagrangian multipliers yields the following linear systemThe right-hand side of eq collects both atomic electronegativities and the electric potential generated by the DFTB density. The latter term accounts for the mutual polarization among the DFTB and FQ portions of the system. In fact, KS equations need to be modified to include the DFTB/FQ contribution to the Hamiltonian matrix, which readswhere the kernel γFQ that takes into account the interaction between the αth Mulliken charge (or the basis function μ) and the ith FQ charge. Notice that the DFTB/FQ term to the total energy and Hamiltonian matrix is the same for both the SCC-DFTB and DFTB3 methods because the Mulliken-based expansion for the density fluctuation δρ does not change. Note finally that the formulation presented above is not limited to FQ but can easily be extended to any kind of variational polarizable MM approach.[58]

Linear Response Regime

The extension of the approach to the linear response regime allows the calculation of some spectral signals and, in particular, vertical transition energies and absorption spectra. The TD-DFTB eigenproblem can be expressed in the Casida formalism as[18]where the eigenvalues ω correspond to excitation energies and the eigenvectors X and Y correspond to single-particle excitations and de-excitation amplitudes. Similarly to DFT/FQ,[34,59−63] to take into account the FQ layer, we need to modify the DFTB response matrices A and B as followswhere the indices i, j and a, b run over the occupied and virtual molecular orbitals with energies ε. K and KFQ are the DFTB and FQ coupling matrices, respectively. K is usually simplified by exploiting the so-called γ-approximation,[18] similarly to the ground state in SCC-DFTB. (See eq .) In such an approximation, the transition density p(r) = ψ(r)ψ(r) is decomposed as a sum of atomic contributions that, after a multipolar expansion, is approximated by means of the monopole term only. (See also eq .) Therefore, K readswhere qα and qβ are Mulliken atomic transition charges. The γ̅ functional is defined as in eq ; however, the functional derivative of Exc is evaluated on ρ. For systems with small charge-transfer effects, γ slightly depends on atomic charges, so that γ̅ can be approximated with its ground-state counterpart γ.

By following refs (34), (45), and (47), the FQ contribution to the coupling matrix can be defined asAt this point, we can exploit the DFTB γ approximation, and eq becomeswhere γ̂ is defined asSimilar to the ground state, we can assume Fα(r) = δ(r – Rα). Thus we obtainwhere the interaction kernel defined in eq is considered.

Computational Details

The equations presented in the previous section were implemented in a modified version of the Amsterdam Molecular Suite (AMS), release 2020.202 program.[23,64] TD-DFTB/FQ calculations were performed on 200 configurations extracted from MD simulations already reported in the literature.[65−67]

For DOX in the gas phase, we ran calculations on the same conformations coming from the MD after removing surrounding water molecules. Table lists an inventory of the different systems/environments addressed in this Article and some other results that will be presented in the following discussion.

Table 1

Inventory of the Number of Atoms, Water Molecules Described at the FQ Level, and Absorption Band Maxima (in Electronvolts) of All Systems under Studya

system	Natoms	NFQwaters	AbsSCC	AbsDFTB3
DOX in gas phase	69	0	2.41	2.43
DOX in water	3714	1215	2.33	2.35
DOX/water/DNA	9756	3100	2.37	2.36
UBI in gas phaseb	1231	0	5.31, 4.75	5.37, 5.05, 4.77
UBI in water	11 731	3500	5.33, 4.80	5.43, 4.86

Experimental results for the π → π* transitions of DOX and UBI in aqueous solution are 2.58 and 4.51 eV, respectively.[68,69].

From the geometry reported in ref (23).

To explore diverse DFTB Hamiltonians, we relied on the Slater–Koster-based DFTB class and performed TD-DFTB calculations for the entire set of snapshots using both the second-order self-consistent charge extension SCC-DFTB (recently also called DFTB2) and the third-order extension known as DFTB3. These two DFTB schemes are thoroughly explained elsewhere (see, e.g., ref (70)), but in short, SCC takes into account density fluctuations with improvements in the description of the polar bonds; likewise, DFTB3 describes hydrogen-bonded complexes and proton affinities, although at a little higher computational cost than SCC-DFTB calculations. SCC-DFTB and DFTB3 TD-DFTB calculations were performed by using mio-1-1[71] and 3ob-3-1[72] parameter sets, respectively.

Finally, to reduce the computational effort of TD-DFTB/FQ-based absorption spectra calculations, we tested the oscillator-strength-based truncation of the single-orbital transition space following the procedures introduced in a previous study.[23] A summary of the technical settings used in the TD-DFTB calculations can be found in Table S1 in the Supporting Information (SI). In all of the cases, absorption spectra profiles were obtained through a convolution of the TD-DFTB excitations by using Gaussian line shapes with a full width at half-maximum (fwhm) value of 0.3 eV, if not explicitly stated. A minimum number of 100 excited states were converged in each calculation. In all calculations, we exploited the FQ parameters proposed in ref (73).

Numerical Results

In this section, we apply DFTB/FQ to describe absorption spectra. We analyze the effect that different choices of the DFTB Hamiltonian and the radius of the DFTB shell have on the spectra and test the accuracy of various intensity selection thresholds for the single orbital transition basis. Also, TD-DFTB/FQ spectra are compared with those obtained by using TD-DFTB calculations in the gas phase.

As test cases, we have chosen two biologically relevant, flexible organic molecules, namely, DOX and UBI, whose structures are depicted in the left panels of Figures and 2, respectively.

Figure 1

Three environments in which UV–vis spectra of doxorubicin were computed in this work. Left: Gas phase. Middle: Snapshot of the molecular dynamics of solvated DOX. Right: Snapshot of the molecular dynamics of DOX intercalated into DNA and surrounded by water molecules.

Figure 2

Environments in which the UV–vis spectra of ubiquitin were computed in this work. Left: Gas-phase conformation by using the same geometry reported in ref (23). Right: Snapshot of the molecular dynamics of solvated UBI is shown, which is treated with the QM/FQ approach.

Doxorubicin

DOX is an anticancer drug,[74] and it is commonly studied in the context of intercalation into DNA due to the proposed mechanism of action based on the insertion of its planar aromatic chromophore portion between sequential base pairs (BPs).[75−77] The drug in water is quite investigated as well given that it travels from the purely aqueous environment to penetrate DNA helices.[65,68,78−81] Regarding the DOX/DNA/water tertiary system, binding energy studies have shown that DOX affinity is sequence-dependent.[82] Although the preferential binding of DOX to double-stranded (ds) DNA is still a subject of debate, recent works reported that among some hexameric evaluated sequences, DOX prefers to bind to the d(CGATCG) in the case of the 1:1 complexes.[67] Therefore, we only discuss the intercalation complex of DOX with that DNA model.

Because of its importance as a chemotherapy medication, there are plenty of works dealing with the spectroscopic evidence of the insertion of a DOX molecule between pairs of nitrogen-containing nucleobases and for the spectral signatures of DOX in aqueous solution. Thus theoretical[65,80,81,83] and experimentally[84−88] obtained absorption spectra can be found in the literature for both environments. The main absorption band around 480 nm has been attributed to a π → π* transition,[68,86] and some bathochromic and hypochromic effects are reported to occur upon intercalation.[86,89] Notwithstanding, it is difficult to observe those shifts because there is a vibronic component dominating the shape of the band.[65] The three environments in which we studied the absorption spectra of DOX are displayed in Figure .

From a set of snapshots (like those shown in Figure ), the array of oscillator strengths obtained at their respective peak positions yields stick spectra (see Figure ) with a natural broadening coming from the dynamical conformations of the chromophore and from the arrangements of the different molecules surrounding the system, that is, water molecules and the DNA basis. Computed stick spectra in the whole range of wavelengths are reported in Figure S1. It should be noted that the intensities of the sticks match the hypochromic effect reported to take place once the intercalation complex is formed. Furthermore, considering that the quality of the results depends on whether there is a convergence of the desired property, some test computations on the UV–vis spectra of DOX in the more complex environment were also performed with an increasing number of snapshots extracted from the MD. Figure shows the convergence of the energy and intensity of the first electronic transition with respect to the number of frames along with the associated 99% confidence intervals. The convergence behavior of the total spectra with respect to the number of frames is reported in Figure S2.

Figure 3

Stick spectra of doxorubicin in vacuo (red line), in water (blue line), and in DNA (green line) performed with different choices of the DFTB Hamiltonian.

Figure 4

Convergence test for the absorption spectra of the tertiary DOX/water/DNA system. The position of the first excitation energy (top panel, red line) and the associated intensity (bottom panel, blue line) calculated with SCC-DFTB and DFTB3 model are reported. As a measure of the convergence, the 99% confidence intervals are reported.

In addition, we evaluated the effect that solute–solvent nonelectrostatic interactions (neglected in the pure DFTB/FQ method) might have on the absorption spectra by adding the solute’s closest water molecules to the QM portion and treating them with one of the DFTB Hamiltonians, whereas we described the remaining solvent molecules by means of FQ. This was done only for DOX in aqueous solution, and the results, together with the average number of water molecules (NQM) for each radius threshold (R), are reported in Table . Clearly, the role of nonelectrostatic (mainly repulsion) effects is minimal.

Table 2

Dependence of the Maximum Absorption Energies of Solvated DOX on the Size of the QM Shella

R (Å)	NQM	ΔVEESCC (eV)	ΔVEEDFTB3 (eV)
1	0	0.00	0.00
2	10	0.01	0.01
3	47	0.01	0.01
4	85	0.02	0.01
5	126	0.02	0.01
6	180	0.03	0.02

ΔVEE is the energy difference with the maximum absorption calculated at QM/FQ level. NQM is the average number of water molecules treated in the QM portion.

DFTB Model Hamiltonians

As mentioned previously, we exploited two of the classic Slater–Koster-based DFTB Hamiltonians, SCC-DFTB and DFTB3. The resulting DFTB/FQ normalized absorption spectra obtained from the average of ∼180 structures of DOX in different environments are plotted in Figure . Table also contains the maximum absorption energies, as obtained from the averaged spectra for both DFTB schemes and for all of the DOX environments under study. Interestingly, both DFTB Hamiltonians offer a similar description of the absorption spectra and the main band attributed to the π → π* transition of the anthracycline chromophore. It should be emphasized that regardless of the DFTB model Hamiltonian and regardless of the environment, the HOMO and LUMO are, for the most part, the orbitals involved in the lowest energy transition. They are graphically depicted in Figure along with other molecular orbitals belonging mainly to the rings of the DOX structure. These results are in line with those obtained by Olszówka et al.,[65] who reported a single HOMO → LUMO transition to be responsible for the appearance of the main band in the absorption spectra.

Figure 5

Absorption spectra of doxorubicin in vacuo (red line), in water (blue line), and in a water/DNA mix (green line) performed with different choices of the DFTB Hamiltonian. Experimental excitation energies from refs (68) and (86) in water and in water/DNA mix are reported with dashed lines.

Figure 6

Most relevant MOs involved in the solvated doxorubicin absorption spectrum. For visualization purposes, virtual orbitals are depicted in different colors.

As can be seen in Figure , the main environment effect is the red shift of the main band moving from the gas phase to water and water/DNA solutions. Overall, TD-DFTB/FQ reproduces the general shape of other published electronic absorption spectra,[81] although the main band is red-shifted by ∼0.2 eV when compared with the experimental maximum absorption energy of solvated DOX. Discrepancies between calculated and experimental results have already been reported for such systems and moderately corrected by using the vertical gradient (VG) or adiabatic Hessian (AH) approaches to vibronically resolve the spectra.[65] Also, as is reported in a recent paper,[90] it would be beneficial to consider different DOX tautomers to obtain a full description of the absorption spectra.

Going from water to water/DNA, there are just slight differences in the vertical excitation energies with both Hamiltonians; however, the spectral profile does exhibit some changes, including a broader main band when DOX is intercalated into DNA and also a different spectral shape at higher energies where the nucleotides are also involved in the electronic transitions. To understand the root causes of these differences in the spectra, especially in the main peak, we have plotted in Figure the molecular orbitals that play a pivotal role in that particular excitation. By analysis and comparison of these orbitals with those displayed in Figure , it becomes clear that the aromatic rings of the DOX’s nearest nucleotides are also participating in the transitions, although the majority of them include the anthraquinone rings (the portion that intercalates between two BPs of dsDNA) and the anchor domains of DOX, where the latter are responsible for stabilizing the DOX–DNA complex via hydrogen bonds with DNA bases. It is worth noting at this point that the HOMO of the ternary DOX/water/DNA system is not localized in the DOX molecule, unlike the situation in which DOX is in the aqueous solution environment.

Figure 7

Most relevant MOs involved in the doxorubicin absorption spectrum when intercalated into DNA (represented by cyan sticks). (a–d) Occupied orbitals. (e) LUMO orbital. For visualization purposes, virtual orbitals are depicted in different colors.

Vertical Excitation Energy Dependence on the Size of the DFTB Shell

As shown in Figure , spectra obtained by varying the number of water molecules in the DFTB layer (QM/QMw/FQ) do not substantially differ each other, which is also confirmed by the data reported in Table , where ΔVEE, that is, the difference in vertical excitation energy (VEE), does not exceed 0.03 eV when compared with the QM/FQ result. It can therefore be argued that regardless of the Hamiltonian choice, the inclusion of the solvent does not play a pivotal role in the description of the bright π → π* transition of solvated DOX; however, the spectral profiles look dissimilar at shorter wavelengths, with a more pronounced contrast in the SCC case.

Figure 8

TD-DFTB/FQ absorption spectra of doxorubicin in water, varying the size of the QM portion and using two different DFTB Hamiltonians. The radius of the DFTB shell is reported in the key.

Intensity Selection Thresholds

Figure shows TD-DFTB calculated absorption spectra of DOX in the gas phase and in aqueous solution, obtained with intensity selection at different oscillator strength thresholds. It should be noticed that the reduced computational cost of the intensity-selected TD-DFTB leads to a loss in accuracy because there is a blue shift of the main band for larger thresholds. (See, for instance, f > 0.1 and f > 0.01.) Nevertheless, when a filter smaller than 0.001 is used, it is evident that the truncation of the basis in oscillator strength has a relatively small effect on the absorption spectrum. In fact, the relative intensities, number of peaks, and peak positions are kept, and the spectrum is practically unaltered compared with the nonfilter case, which is valid for both Hamiltonians. These findings indicate that a large part of the basis has a minor contribution to the spectra, as already reported for the simulation of the absorption spectra of C60 fullerene, Ir(ppy)3, and UBI.[23]

Figure 9

TD-DFTB and TD-DFTB/FQ absorption spectra of doxorubicin in vacuo (left panel) and in water (right panel) performed with different choices of the DFTB Hamiltonian and changing the intensity-selection thresholds for the single-orbital transitions basis. N.F. stands for no filter or that all single orbital transitions were considered.

Ubiquitin

UBI is a 76-amino acid polypeptide (1231 atoms) with diverse roles, mainly oriented to help in the regulation of the processes of other proteins in the body.[91−94] This small protein has been considered as a universal constituent of living cells.[95] Structurally speaking, UBI contains important chromophores like tyrosine and phenylalanine, with the former presenting higher absorbance. As a matter of fact, UBI has served as a model protein to study the sensitivity of UV–visible spectroscopy to environmental factors.[69,96]

DFTB/FQ is challenged in this section to compute the UV–vis absorption spectra of UBI in aqueous solution. The entire protein has been treated at the DFTB level, whereas water molecules are described by means of the FQ force field. Two major features are visible in the UBI experimental spectra in solution and in the gas phase, as reported by Bellina et al.:[69] (i) a broad band centered around 275–280 nm and (ii) an intense response at high energy with an onset at 250 nm. Indeed, it has been found that aromatic amino acids and proteins absorb UV light and show two main bands in UV–vis spectra, one centered on 280 nm that is the result of absorbance by the aromatic ring portion of their structure and a second one at lower wavelengths, which stems from the absorbance of peptide and carboxylic acid moieties. Because of this, it is not surprising that for UBI, a polypeptide containing tyrosine, the same bands are found. In particular, when Tyr is in aqueous solution, absorption maxima appear at ∼220 (higher absorbance) and 275 nm;[97−99] some authors have postulated that the two bands are probably arising from two well-separated π → π* transitions.[100,101]

The influence of the environment on the absorption spectra of the UBI protein has already been demonstrated.[96] This effect can also be observed in Figure , which shows the spectra in the gas phase and in aqueous solution and a comparison between the results coming from the two model Hamiltonians. Such spectra have been obtained through a Gaussian convolution of the TD-DFTB excitations using an fwhm value of 0.2 eV. As a reference, the computed stick spectra of UBI in the gas phase and in aqueous solution over the entire spectral range are reported in Figure S3.

Figure 10

Comparison between calculated absorption spectra of UBI in the gas phase and aqueous solution, as obtained with different DFTB model Hamiltonians. The experimental data of UBI in the gas phase from ref (69) are reported with green blocks representing the experimental range of the two main absorption bands.

It can be observed that both SCC and DFTB3 yield the same spectral shape in the case of solvated UBI, whereas a different behavior is observed in the gas phase. Such discrepancies may be attributed to different parametrizations exploited in the two approaches.[102] Additionally, the presence of water has a tiny but non-negligible effect on the absorption spectra, overall red-shifting the main peaks that appear in gas phase spectra. The same calculations have also been performed by employing the nonpolarizable TIP3P force field to describe water molecules. (See Figure S4.) Absorption band maxima of UBI are also summarized in Table . To assign the transitions in the full protein, we applied a filter on the strongest oscillator strengths in the region of the maximum absorbance, and we quantified the contribution of single orbital transitions, thus identifying the leading contributing molecular orbitals and the part of the protein where they were situated. As has been indicated above, the band centered around 280 nm can be assigned to a π → π* excitation, and even though orbitals from very distinct residues can participate in the transitions, the predominant ones are localized on the aromatic tyrosine and phenylalanine chromophores, with tyrosine being responsible for most of the absorbance. A more in-depth analysis indicates that frontier molecular orbitals (HOMO, HOMO–1, LUMO, and LUMO+1) in the extended tyrosine residue resemble those involved in the absorption band at 280 nm of the UBI protein. The considered orbitals are presented in Figure and are responsible for 90% of the state, having the greatest oscillator strengths in the region of the π → π* excitation.

Figure 11

Left panel: Chromophores tyrosine (framed in black) and phenylalanine, which are the dominant sources of the UBI protein absorbance. Right panel: Molecular orbitals of extended Tyr mainly involved in the transitions. (See the text.) Extended Tyr means that in the calculations, the residue was capped with N-terminal acetyl and C-terminal N-Me amide capping groups to preserve the peptide bonds inside the protein. For visualization purposes, virtual orbitals are depicted in different colors.

We move on to compare our results with experimental data on spectral shifts going from solution to the gas phase. It was previously determined that the π → π* band in the gas phase is red-shifted as compared with the absorption in solution.[69] In addition, it is known that tyrosine is located at the surface of the protein[103] and thus has a strong effect on the environment, and the addition of water is also anticipated to cause a red shift of the tyrosine absorption spectrum.[104] In this context, it should be noted that although our results are in agreement with these observations, the experimental final shift is not fully reliable due to the fact that UV gas-phase spectra were measured for the isolated deprotonated protein and UBI can change its conformation going from the gas phase to solution.[69]

Finally, we remark that it is necessary to ensure the analysis to be performed on final converged spectra. Inspection of UBI spectra in aqueous solution, as obtained by averaging out a varying number of frames (Figure ), reveals that increasing the number of frames has no impact on the final band shape; in fact, a sampling of 50 frames is sufficient to achieve convergence, with no missing features appearing in the spectra. Lastly, if single orbital transitions with an oscillator strength smaller than 0.001 are removed (results not shown here), then the spectra do not change, and as expected, the computational effort decreases.

Figure 12

TD-DFTB/FQ absorption spectra of ubiquitin in aqueous solution, as obtained by averaging out an increasing number of snapshots.

Summary, Conclusions, and Future Perspectives

We have presented a novel polarizable QM/MM approach where DFTB is coupled to the polarizable FQ force field. The model, which has been extended to TD-DFTB linear response, permits us to treat large systems in the condensed phase thanks to the favorable scaling of DFTB as compared with standard DFT and other ab initio methods. DFTB/FQ has been applied to the simulation of the electronic absorption spectra of DOX and UBI in different environments, showing that the inclusion of the FQ layer strongly affects spectral shapes and accounts for changes in both peaks’ positions and relative intensities when going from the gas phase to condensed phase.

The low computational cost of DFTB has allowed for a detailed analysis of the role of nonelectrostatic effects in the case of DOX in aqueous solution. In particular, the size of the DFTB portion has been varied by adding up a limited number of water molecules, showing that the DOX electronic response is almost unaffected by increasing the radius of the DFTB portion above 4 Å. This confirms the short-range nature of the nonelectrostatic interactions that are naturally included within the DFTB portion and that are dominated by Pauli repulsion. In addition, we show that by selecting the most intense spectral bands only, the accuracy of computed spectra is not particularly affected; however, the computational time of TD-DFTB/FQ calculations is substantially reduced. DOX and UBI spectral profiles have been obtained by taking into account solute and solvent dynamics; the phase-space sampling and the consequent configurational variability in both solute and solvent moieties have brought up natural broadening in absorption bands, which is directly obtained from the signals arising on the different snapshots extracted from MD simulations.

Finally, the results reported in this work show that the combination of DFTB and FQ permits us to model absorption spectra of large molecules embedded in complex environments at a low computational cost and in nice agreement with experimental data. This opens up the opportunity to explore more challenging spectroscopies in different environments. In fact, DFTB/FQ (similarly to other QM/polarizable MM approaches) can simulate molecular properties in any kind of complex environment, pending appropriate parametrization.

To validate the accuracy of DFTB/FQ to describe solvatochromic effects, we have also computed vertical excitation energies at the TD-DFT/FQ level (see Table S2), in line with the preliminary calculations of some of the present authors.[80,81,105] DFTB/FQ and TD-DFT/FQ values are very similar, thus demonstrating that DFTB/FQ gives a correct description, from both a qualitative and quantitative point of view, for our systems. However, a more extended benchmark analysis on several systems would be required to finally validate the accuracy of the approach, as solvatochromic effects are strongly dependent on both the system and the nature of the excitation.

Improvement in the numerical performance of DFTB/FQ can be achieved by reparametrizing the FQ force field for DFTB calculations for both aqueous and nonaqueous solutions, in line with previous studies of some of the present authors.[46,106] Also, the description of the environment can be refined by adding polarizable dipole moments on MM atoms, such as in the recently developed FQFμ approach, which appropriately simulated anisotropic solvent effects.[35,37] In addition, in this Article, we have fully relied on a purely electrostatic model. Although electrostatics often dominates solvation effects, we have recently shown that nonelectrostatic interactions, in particular, Pauli repulsion, may strongly affect computed molecular properties.[107] Such terms can be included in the DFTB/FQ approach by following recent works of our group.[48,107]