TAO-DFT-Based Ab Initio Molecular Dynamics

Recently, AIMD (ab initio molecular dynamics) has been extensively employed to explore the dynamical information of electronic systems. However, it remains extremely challenging to reliably predict the properties of nanosystems with a radical nature using conventional electronic structure methods (e.g., Kohn-Sham density functional theory) due to the presence of static correlation. To address this challenge, we combine the recently formulated TAO-DFT (thermally-assisted-occupation density functional theory) with AIMD. The resulting TAO-AIMD method is employed to investigate the instantaneous/average radical nature and infrared spectra of n-acenes containing n linearly fused benzene rings (n = 2–8) at 300 K. According to the TAO-AIMD simulations, on average, the smaller n-acenes (up to n = 5) possess a nonradical nature, and the larger n-acenes (n = 6–8) possess an increasing radical nature, showing remarkable similarities to the ground-state counterparts at 0 K. Besides, the infrared spectra of n-acenes obtained with the TAO-AIMD simulations are in qualitative agreement with the existing experimental data.

1. Introduction

Molecular dynamics (MD) is a computational method for simulating dynamical processes that occur in a system consisting of atoms (for example, atoms, molecules, solids, and liquids) (Lifson and Warshel, 1968; Levitt and Lifson, 1969; Karplus and Petsko, 1990; Kresse and Hafner, 1993; Sprik et al., 1996; Silvestrelli et al., 1997; Putrino and Parrinello, 2002; Tuckerman, 2002; Chai et al., 2003; Kuo and Mundy, 2004; Jensen, 2007; Marx and Hutter, 2009; Gaigeot, 2010; Ramírez-Solís et al., 2011; Vitale et al., 2015; Hollingsworth et al., 2018). The motion of the atomic nuclei in the system is described by the classical Newtonian equations of motion, starting from pre-specified initial conditions (e.g., initial nuclear positions and velocities) and subject to boundary conditions suitable for the system studied. By performing MD simulations, one can calculate both the dynamical and equilibrium thermodynamic properties associated with a system at non-zero temperatures and can simultaneously monitor the microscopic movements of the atomic nuclei in the system. However, the specification of the forces acting on the atomic nuclei in MD simulations remains challenging, yielding two popular types of MD methods: classical MD and ab initio MD (AIMD).

In classical MD, the forces that act on the atomic nuclei are calculated using an empirical potential energy function (i.e., the potential energy of a system is expressed as a function of the nuclear coordinates) that is defined by molecular mechanics (Tuckerman, 2002; Jensen, 2007; Marx and Hutter, 2009). Because of a simple analytical function for the potential energy, classical MD is computationally very efficient, and it has thus been widely employed for many applications (Lifson and Warshel, 1968; Levitt and Lifson, 1969; Karplus and Petsko, 1990; Hollingsworth et al., 2018). For example, classical MD can reveal the details of the motion of macromolecules (e.g., proteins) (Lifson and Warshel, 1968; Levitt and Lifson, 1969), which can be very challenging using experimental approaches. However, in classical MD, empirical potential energy functions are typically parameterized for some peculiar systems and may not be transferable to others (Tuckerman, 2002). Furthermore, in classical MD, electrons are not present explicitly (i.e., their effects are approximated by empirical potential energy functions), and the electronic properties (e.g., electron density and radical nature) of systems therefore cannot be explored.

To overcome the shortcomings of classical MD, one can resort to AIMD, wherein the forces that act on the atomic nuclei are calculated on-the-fly according to the potential energy obtained with an electronic structure method along an AIMD trajectory (Tuckerman, 2002; Jensen, 2007; Marx and Hutter, 2009). Clearly, the accuracy of AIMD simulations can be limited due to the accuracy of the underlying electronic structure method. On the other hand, for AIMD simulations, it is also essential to adopt a computationally efficient electronic structure method that can provide reasonably accurate potential energy for the systems of interest, since there are typically 10⁴–10⁶ time steps per AIMD trajectory. Owing to the decent balance between performance and cost, KS-DFT (i.e., Kohn-Sham density functional theory) (Kohn and Sham, 1965) is currently the most popular electronic structure method adopted in AIMD simulations. Over the past three decades, KS-DFT-based AIMD (KS-AIMD) has been commonly used to investigate the equilibrium thermodynamic and dynamical properties associated with electronic systems at finite temperatures, and has been a powerful AIMD method for modeling various phenomena (Kresse and Hafner, 1993; Sprik et al., 1996; Silvestrelli et al., 1997; Putrino and Parrinello, 2002; Chai et al., 2003; Kuo and Mundy, 2004; Gaigeot, 2010; Ramírez-Solís et al., 2011; Vitale et al., 2015).

Nevertheless, KS-DFT using traditional exchange-correlation (XC) energy functionals could yield qualitative failures for the various properties associated with electronic systems possessing radical nature due to the lack of static correlation effects (Cohen et al., 2008, 2012). For example, consider linear n-acene (with the chemical formula C_4n+2H_2n+4), which contains n linearly fused benzene rings (for example, see Figure 1), N_e = 26n + 16 electrons, and N = 6n+6 nuclei. According to the recent findings (Hachmann et al., 2007; Chai, 2012, 2014, 2017; Rivero et al., 2013; Wu and Chai, 2015; Fosso-Tande et al., 2016), the larger n-acenes (n ≥ 6) possess increasing radical nature in their ground states. However, KS-DFT adopting commonly used XC energy functionals could perform poorly for the larger n-acenes (n ≥ 6) due to the lack of static correlation effects (Hachmann et al., 2007; Chai, 2012, 2017). It can therefore be anticipated that KS-AIMD simulations at finite temperatures may yield unreliable dynamical information for electronic systems with a radical nature [e.g., the larger n-acenes (n ≥ 6)].

FIGURE 1

Figure 1. Structure of 8-acene, containing eight linearly fused benzene rings.

To study the properties of electronic systems with a radical nature, one generally resorts to multi-reference (MR) electronic structure methods (Andersson et al., 1992; Hachmann et al., 2007; Gidofalvi and Mazziotti, 2008; Pelzer et al., 2011; Gryn'ova et al., 2015; Fosso-Tande et al., 2016; Battaglia et al., 2017; Mullinax et al., 2019). Although MR electronic structure methods can reliably predict the various properties associated with electronic systems with radical nature, they could be prohibitively expensive for a single-point energy + nuclear gradient calculation on large electronic systems, not to mention the respective AIMD simulations where such calculations should be performed about 10⁴–10⁶ times per AIMD trajectory. To explore the dynamical properties of nanosystems with radical nature using AIMD simulations at finite temperatures, therefore, it is essential to employ an efficient electronic structure method that can properly describe the static correlation effects during the AIMD simulations.

Recently, TAO-DFT (thermally-assisted-occupation density functional theory) (Chai, 2012) has been formulated for investigating the ground-state (GS) properties associated with nanosystems possessing a radical nature at 0 K. Note that TAO-DFT is a density functional theory adopting fractional orbital occupation numbers (i.e., generated by the Fermi-Dirac distribution using a fictitious temperature θ), which is rather different from KS-DFT. In TAO-DFT, the fictitious temperature θ is intimately correlated with the effect of the configuration mixing on the electron density of a GS system (e.g., see sections 2, 3 of Chai, 2012), and it is thus completely unrelated to the physical temperature (i.e., 0 K) of the GS system. Note that KS-DFT corresponds to TAO-DFT with θ = 0. However, for a general GS system at 0 K, the θ value in TAO-DFT can be nonzero (e.g., see sections 3.5, 4 of Chai, 2012). In TAO-DFT, an entropy contribution term, which is dependent on the fictitious temperature θ and orbital occupation numbers, can offer a proper description for static correlation even when the simplest LDA (i.e., local density approximation) XC energy functional is used (Chai, 2012). At the LDA level, even with a properly defined system-independent θ value (e.g., see section 5 of Chai, 2012), TAO-DFT, which has similar cost as KS-DFT in computation, has been shown to consistently outperform KS-DFT for electronic systems with a radical nature, while performing similarly to KS-DFT for electronic systems with a nonradical nature (Chai, 2012).

More complicated semilocal (Chai, 2014), global hybrid (Chai, 2017), and range-separated hybrid (Chai, 2017; Xuan et al., 2019) XC energy functionals could also be used in TAO-DFT. Moreover, in order to enhance its accuracy for a wide range of applications, a method that determines the θ value in TAO-DFT in a self-consistent manner has recently been developed (Lin et al., 2017). Very recently, a frequency-domain formulation of linear-response time-dependent TAO-DFT (Yeh et al., 2020) has been formulated to explore the properties of electronic excited states within TAO-DFT.

Owing to its computational efficiency and decent accuracy for exploring the properties of electronic systems at the nanoscale, TAO-DFT has recently been adopted to investigate the electronic properties (Wu and Chai, 2015; Seenithurai and Chai, 2016, 2017, 2018, 2019, 2020; Wu et al., 2016; Yeh and Chai, 2016; Yeh et al., 2018; Chung and Chai, 2019; Deng and Chai, 2019; Hanson-Heine and Hirst, 2020; Hanson-Heine et al., 2020; Huang et al., 2020; Manassir and Pakiari, 2020), hydrogen storage properties (Seenithurai and Chai, 2016, 2017, 2018), and vibrational frequencies (Hanson-Heine, 2020) of several electronic systems at the nanoscale, especially for those possessing a radical nature.

Therefore, in the present work, we propose to combine TAO-DFT (Chai, 2012) with AIMD, yielding TAO-DFT-based AIMD (TAO-AIMD). Since analytical nuclear gradients for TAO-DFT are available (Chai, 2012), TAO-AIMD is as computationally efficient as KS-AIMD. Accordingly, it is feasible to study the equilibrium thermodynamic and dynamical properties of nanosystems with a radical nature using TAO-AIMD simulations at finite temperatures. To highlight some of the present capabilities of TAO-AIMD, we perform TAO-AIMD simulations to explore the instantaneous/average radical nature and infrared (IR) spectra of n-acenes with n = 2–8 fused benzene rings at 300 K. The rest of this paper is organized as follows. The TAO-AIMD method is defined in section 2. We then describe the computational details in section 3, discuss the results in section 4, and give our conclusions in section 5.

2. TAO-AIMD

Consider a system containing N_e electrons (described by coordinates r₁, …, r_Ne) and N nuclei (described by coordinates R₁, …, R_N). Here, we resort to the adiabatic or Born-Oppenheimer (BO) approximation (Born and Oppenheimer, 1927). Because the electrons are much lighter than the nuclei, it is assumed that the nuclei move relatively slowly, and hence the electrons are able to respond to the nuclear motion almost instantaneously. In other words, the electronic motion and nuclear motion can be treated separately.

Accordingly, in the first step of the BO approximation, the kinetic energy of the nuclei is ignored. For fixed nuclear positions R₁, …, R_N, the electronic Hamiltonian is expressed as (Tuckerman, 2002; Jensen, 2007; Marx and Hutter, 2009)

$\begin{array}{l}{Ĥ}_{\text{elec}}=-\frac{{\hslash }^2}{2m_e}{\displaystyle \sum _{i=1}^{N_e}}{\nabla }_i^2-\frac{e^2}{4\pi {ϵ}_0}{\displaystyle \sum _{i=1}^{N_e}}{\displaystyle \sum _{A=1}^N}\frac{Z_A}{|{\text{r}}_i-{\text{R}}_A|}\\ +\frac{e^2}{4\pi {ϵ}_0}{\displaystyle \sum _{i=1}^{N_e}}{\displaystyle \sum _{j>i}^{N_e}}\frac{1}{|{\text{r}}_i-{\text{r}}_j|},& (1)\end{array}$

where m_e and −e are the mass and charge, respectively, of an electron, and Z_Ae is the charge of nucleus A. On the right-hand side of Equation (1), the first term is the kinetic energy of electrons, the second term is the nuclear-electron attraction energy, and the third term is the electron-electron repulsion energy. The time-independent electronic Schrödinger equation

$\begin{array}{ll}{Ĥ}_{\text{elec}}{\text{Ψ}}_k=E_k{\text{Ψ}}_k& (2)\end{array}$

is subsequently solved for the electronic energy E_k and electronic wavefunction Ψ_k. In particular, the lowest eigenvalue E₀ is the GS electronic energy, and the corresponding eigenfunction Ψ₀ is the GS electronic wavefunction. By adding the nuclear-nuclear repulsion energy to the electronic energy, one obtains the potential energy for the k-th electronic eigenstate of the system

$\begin{array}{ll}{U}_k=E_k+\frac{e^2}{4\pi {ϵ}_0}{\displaystyle \sum _{A=1}^N}{\displaystyle \sum _{B>A}^N}\frac{Z_A{Z}_B}{|{\text{R}}_A-{\text{R}}_B|}.& (3)\end{array}$

By varying the nuclear positions and solving the corresponding time-independent electronic Schrödinger equation, U_k can be expressed as a function of the nuclear positions, also known as the potential energy surface of the k-th electronic eigenstate. In the second step of the BO approximation, the kinetic energy of the nuclei is reintroduced, and in principle, the nuclear dynamics is described by the time-dependent nuclear Schrödinger equation, evolving on a potential energy surface, e.g., U_k(R₁, …, R_N).

Similar to most AIMD simulations (Tuckerman, 2002; Jensen, 2007; Marx and Hutter, 2009; Gaigeot, 2010; Ramírez-Solís et al., 2011; Vitale et al., 2015), in this work, we assume that nuclear quantum effects could be ignored, and nuclear motion occurs only on the GS potential energy surface (i.e., the potential energy surface of the electronic GS), U₀(R₁, …, R_N). Since the time-independent electronic Schrödinger equation is only applicable to very small electronic systems, a sufficiently efficient electronic structure method for the determination of electronic ground state is typically required for AIMD simulations. As mentioned previously, KS-AIMD simulations can be unreliable for electronic systems with a radical nature. On the other hand, AIMD simulations employing MR electronic structure methods are very likely to be computationally intractable for most electronic systems.

To resolve this issue with minimum computational expense, in the present work, we combine TAO-DFT (Chai, 2012) with AIMD, yielding the TAO-AIMD method. Specifically, in TAO-AIMD, the nuclei are treated as classical particles, obeying the classical nuclear Hamiltonian

$\begin{array}{l}{H}_{\text{nucl}}({\text{P}}_1(t),...,{\text{P}}_N(t),{\text{R}}_1(t),...,{\text{R}}_N(t))={\displaystyle \sum _{A=1}^N}\frac{|{\text{P}}_A(t){|}^2}{2M_A}\\ +U_0^{\text{TAO-DFT}}({\text{R}}_1(t),...,{\text{R}}_N(t)),& (4)\end{array}$

where M_A is the mass of nucleus A, and P_A(t) is the momentum of nucleus A at time t. The left-hand side of Equation (4) gives the total energy at time t [i.e., E(t)]. On the right-hand side of Equation (4), the first term is the nuclear kinetic energy at time t, and the second term is the potential energy of the electronic GS obtained with TAO-DFT for the nuclear positions at time t, i.e., R₁(t), …, R_N(t):

$\begin{array}{l}{U}_0^{\text{TAO-DFT}}({\text{R}}_1(t),...,{\text{R}}_N(t))=E_0^{\text{TAO-DFT}}({\text{R}}_1(t),...,{\text{R}}_N(t))\\ +\frac{e^2}{4\pi {ϵ}_0}{\displaystyle \sum _{A=1}^N}{\displaystyle \sum _{B>A}^N}\frac{Z_A{Z}_B}{|{\text{R}}_A(t)-{\text{R}}_B(t)|},& (5)\end{array}$

where $E_0^{\mathrm{\text{TAO-DFT}}}({\text{R}}_1(t),...,{\text{R}}_N(t))$ (e.g., see Equation 27 of Chai, 2012) is the corresponding GS electronic energy obtained with TAO-DFT. On the basis of Equation (4), the nuclei move based on Newton's equations of motion on the GS potential energy surface generated by TAO-DFT:

$\begin{array}{ll}{\stackrel{.}{\text{R}}}_A(t)=\frac{{\text{P}}_A(t)}{M_A}& (6)\end{array}$

$\begin{array}{ll}{\stackrel{.}{\text{P}}}_A(t)=-{\nabla }_A{U}_0^{\text{TAO-DFT}}({\text{R}}_1(t),...,{\text{R}}_N(t)),& (7)\end{array}$

where ${\stackrel{.}{\text{R}}}_A(t)$ is the velocity of nucleus A at time t, ${\stackrel{.}{\text{P}}}_A(t)$ is the time derivative of the momentum of nucleus A at time t, and the right-hand side of Equation (7) gives the force acting on nucleus A at time t [i.e., F_A(t)].

Equations (4)–(7) form the theoretical basis of the TAO-AIMD method. Given the initial nuclear positions R₁(0), …, R_N(0) and velocities ${\stackrel{.}{\text{R}}}_1(0),...,{\stackrel{.}{\text{R}}}_N(0)$, all the future nuclear positions R₁(t), …, R_N(t) and velocities ${\stackrel{.}{\text{R}}}_1(t),...,{\stackrel{.}{\text{R}}}_N(t)$ are determined by Equations (5)–(7), generating a TAO-AIMD trajectory (i.e., TAO-AIMD is deterministic). Note that the GS potential energy surface and the forces that act on the nuclei can be computed on-the-fly using TAO-DFT, as needed along the TAO-AIMD trajectory. According to the definitions presented here (Marx and Hutter, 2009), TAO-AIMD can also be regarded as TAO-DFT-based Born-Oppenheimer MD (BOMD), i.e., TAO-BOMD.

Note that the GS potential energy surface and the forces that act on the nuclei are computed using TAO-DFT in TAO-AIMD, while they are computed using KS-DFT in KS-AIMD. Since TAO-DFT is as computationally efficient as KS-DFT (e.g., see Chai, 2012 for details), TAO-AIMD is similar to KS-AIMD in computational expense. On the other hand, AIMD simulations employing MR electronic structure methods are very likely to be computationally infeasible for most electronic systems. Accordingly, TAO-AIMD can be a very promising method for exploring the equilibrium thermodynamic and dynamical properties of nanosystems with a radical nature at finite temperatures. In addition, existing XC energy functionals defined in KS-DFT could also be employed in TAO-DFT (Chai, 2012, 2014, 2017) and TAO-AIMD. For electronic systems with a nonradical nature, TAO-DFT performs similarly to KS-DFT, and TAO-AIMD should therefore perform similarly to KS-AIMD.

3. Computational Details

All calculations are performed using TAO-LDA (Chai, 2012), i.e., TAO-DFT employing the LDA XC and θ-dependent energy functionals, where the suggested fictitious temperature θ = 7 mhartree is employed (Chai, 2012). Almost all numerical data are obtained with Q-Chem 4.4 (Shao, 2015), using the 6-31G(d) basis set and a numerical grid consisting of 75 Euler-Maclaurin radial grid points and 302 Lebedev angular grid points, wherein isolated boundary conditions (i.e., well-suited for the study of atoms and molecules) are employed. The TAO-AIMD-based IR spectra are computed using the program TRAVIS (Brehm and Kirchner, 2011; Thomas et al., 2013, 2015; Brehm et al., 2020).

In several recent investigations (Chai, 2012; Wu and Chai, 2015; Fosso-Tande et al., 2016; Yeh and Chai, 2016; Mullinax et al., 2019), the orbital occupation numbers obtained from TAO-LDA (with θ = 7 mhartree) have been found to be similar to the natural orbital occupation numbers obtained from a very accurate MR electronic structure method that can be applied to treat relatively large active spaces, leading to a qualitatively similar tendency for the radical nature associated with various PAHs (i.e., polycyclic aromatic hydrocarbons).

On the other hand, because of the constraint of symmetry, the spin-unrestricted and spin-restricted calculations based on an exact electronic structure method should lead to the same energy values for the lowest singlet state (i.e., GS) of n-acene (Chai, 2012; Rivero et al., 2013; Gryn'ova et al., 2015). Nevertheless, KS-DFT employing conventional XC energy functionals fails to obey this constraint for the larger n-acenes (which are electronic systems possessing radical nature), leading to the unphysical symmetry-breaking effects in the corresponding spin-unrestricted calculations (Cohen et al., 2008, 2012). In our previous studies (Chai, 2012, 2014; Wu and Chai, 2015), the spin-unrestricted and spin-restricted GS (i.e., lowest singlet state) energy values of n-acene (up to n = 100), calculated by TAO-LDA (with θ = 7 mhartree), have been found to be essentially the same, leading to essentially no unphysical symmetry-breaking effects in the corresponding spin-unrestricted calculations. For computational efficiency, therefore, all TAO-LDA calculations in this work are spin-restricted calculations unless noted otherwise.

For all TAO-AIMD simulations, a time step of 20 a.u. (≈ 0.484 fs) is adopted for the integration of the equations of motion. For each TAO-AIMD simulation, the initial geometry of n-acene is chosen as the GS geometry of n-acene obtained with TAO-LDA, and the initial nuclear velocities of n-acene are randomly selected from the Maxwell-Boltzmann (MB) distribution at T = 300 K. To equilibrate n-acene at T = 300 K, the TAO-AIMD simulation is first performed in the canonical (NVT) ensemble with the aid of the Nosé-Hoover (NH) chain thermostat (Martyna et al., 1992) (with the NH chain length of three auxiliary variables and NH timescale of 200 fs in the implementation of Q-Chem) for about 10.2 ps. To avoid any interference with the dynamics (Gaigeot and Sprik, 2003), we subsequently remove the NH chain thermostat and continue the TAO-AIMD simulation in the microcanonical (NVE) ensemble for a total of 10,500 time steps (≈ 5.1 ps), at which time the total energy of n-acene is well converged. After equilibration, we continue the TAO-AIMD simulation in the NVE ensemble and collect relevant data along the equilibrated TAO-AIMD trajectory for a total of 42,000 time steps (≈ 20.3 ps), where the average temperature is 300 ± 1 K. As the initial nuclear velocities of n-acene for each TAO-AIMD simulation are randomly selected from the MB distribution at 300 K, the aforementioned processes are repeated to generate a total of four different equilibrated TAO-AIMD trajectories (≈ 20.3 ps per trajectory).

Furthermore, to demonstrate the significance of TAO-AIMD simulations for exploring the dynamical information of large molecules with radical nature, we also perform preliminary calculations on 8-acene to examine the possible symmetry-breaking effects in the spin-unrestricted TAO-AIMD/KS-AIMD simulations (see Supplementary Section 1 and Supplementary Figures 1–3).

4. Results and Discussion

4.1. Symmetrized von Neumann Entropy

In TAO-DFT, the radical nature of a GS molecule can be examined by the symmetrized von Neumann entropy (for example, see Equation 5 of Chung and Chai, 2019 for the spin-unrestricted case). For the spin-restricted case, the symmetrized von Neumann entropy of a GS molecule can be expressed as

$\begin{array}{ll}{S}_{\text{vN}}=-{\displaystyle \sum _{i=1}^{\infty }}\left\{\frac{f_i}{2}\text{ln}\left(\frac{f_i}{2}\right)+\left(1-\frac{f_i}{2}\right)\text{ln}\left(1-\frac{f_i}{2}\right)\right\},& (8)\end{array}$

where f_i (i.e., a number between 0 and 2) is the occupation number of the i-th orbital of the GS molecule, obtained with spin-restricted TAO-DFT. Note that f_i is closely related to the corresponding natural orbital occupation number (Löwdin and Shull, 1956; Chai, 2012, 2014, 2017). For a GS molecule possessing a nonradical nature, the occupation numbers associated with all orbitals are very close to either 0 or 2, yielding a vanishingly small S_vN value. Nonetheless, for a GS molecule with a significant radical nature, the active orbital occupation numbers can deviate significantly from 0 and 2 (for example, 0.2–1.8); hence, the corresponding S_vN value can greatly increase as the number of active orbitals increases and/or the active orbital occupation numbers are closer to 1 (Rivero et al., 2013; Chai, 2014, 2017; Wu and Chai, 2015; Seenithurai and Chai, 2016, 2017, 2018, 2019; Wu et al., 2016; Yeh et al., 2018; Chung and Chai, 2019; Deng and Chai, 2019; Huang et al., 2020).

On the basis of Equation (8), in a spin-restricted TAO-AIMD simulation, the symmetrized von Neumann entropy of a molecule at time t along a TAO-AIMD trajectory can be defined as

$\begin{array}{ll}{S}_{\text{vN}}(t)=-{\displaystyle \sum _{i=1}^{\infty }}\left\{\frac{f_i(t)}{2}\text{ln}\left(\frac{f_i(t)}{2}\right)+\left(1-\frac{f_i(t)}{2}\right)\text{ln}\left(1-\frac{f_i(t)}{2}\right)\right\},& (9)\end{array}$

where f_i(t) (i.e., a value between 0 and 2) is the occupation number of the i^th orbital of the molecule, obtained with spin-restricted TAO-DFT for the nuclear positions at time t, i.e., R₁(t), …, R_N(t), along the TAO-AIMD trajectory. Accordingly, the time average of S_vN(t) along the TAO-AIMD trajectory is calculated by

$\begin{array}{ll}\overline{S_{\text{vN}}}=\frac{1}{\tau }{\displaystyle {\int }_0^{\tau }}{S}_{\text{vN}}(t)dt,& (10)\end{array}$

where τ is the total time duration of the TAO-AIMD trajectory.

To investigate the radical nature of n-acene in the TAO-AIMD simulations, the $\overline{S_{\mathrm{\text{vN}}}}$ value of n-acene along each equilibrated TAO-AIMD trajectory is computed using Equation (10), and the reported $\overline{S_{\mathrm{\text{vN}}}}$ value of n-acene is an average over four different equilibrated TAO-AIMD trajectories (≈ 20.3 ps per trajectory). For comparison, we also report the S_vN value of GS n-acene (given by Equation 8), corresponding to the $\overline{S_{\mathrm{\text{vN}}}}$ value of n-acene at 0 K in this work.

As presented in Figure 2, the $\overline{S_{\mathrm{\text{vN}}}}$ value of n-acene obtained with the TAO-AIMD simulations at 300 K is very close to the S_vN value of GS n-acene Chai (2014), increasing monotonically with increasing n (also see Supplementary Table 1). This suggests that similar to the GS counterparts at 0 K (Hachmann et al., 2007; Chai, 2012, 2014, 2017), on average, the smaller n-acenes (up to n = 5) possess nonradical nature, and the larger n-acenes (n = 6–8) possess increasing radical nature in AIMD simulations at 300 K.

FIGURE 2

Figure 2. Symmetrized von Neumann entropy of n-acene, obtained with the TAO-AIMD simulations at 300 K and GS calculation, calculated by TAO-LDA.

Along each equilibrated TAO-AIMD trajectory, the instantaneous S_vN(t) value of n-acene fluctuates over time and can thus be even larger than the $\overline{S_{\mathrm{\text{vN}}}}$ value of n-acene (see Figure 3 for n = 8, and Supplementary Figures 4–9 for others). Therefore, for the larger n-acenes (n = 6–8), it is essential to perform AIMD simulations at 300 K with an efficient electronic structure method that can reliably describe strong static correlation effects, well-justifying the use of TAO-AIMD in this work.

FIGURE 3

Figure 3. Time evolution of the symmetrized von Neumann entropy of 8-acene, obtained from four different equilibrated TAO-AIMD trajectories (No.1 to No.4) at 300 K, calculated by TAO-LDA. The TAO-AIMD average and GS values are also shown for comparison.

4.2. Active Orbital Occupation Numbers

As mentioned previously, in a spin-restricted TAO-AIMD simulation, f_i(t) (i.e., a value between 0 and 2) is the occupation number of the i-th orbital of a molecule, obtained with spin-restricted TAO-DFT for the nuclear positions at time t, i.e., R₁(t), …, R_N(t), along the TAO-AIMD trajectory. Consequently, the time average of f_i(t) along the TAO-AIMD trajectory is calculated by

$\begin{array}{ll}\overline{f_i}=\frac{1}{\tau }{\displaystyle {\int }_0^{\tau }}{f}_i(t)dt,& (11)\end{array}$

where τ is the total time duration of the TAO-AIMD trajectory.

To further assess the radical nature of n-acene in the TAO-AIMD simulations, the $\overline{f_i}$ value of n-acene along each equilibrated TAO-AIMD trajectory is computed using Equation (11), and the reported $\overline{f_i}$ value of n-acene is an average over four different equilibrated TAO-AIMD trajectories (≈ 20.3 ps per trajectory). For comparison, we also present the f_i value of GS n-acene, which corresponds to the $\overline{f_i}$ value of n-acene at 0 K in this work.

For n-acene (containing N_e electrons), we define the HOMO (i.e., highest occupied molecular orbital) as the (N_e/2)th orbital, the LUMO (i.e., lowest unoccupied molecular orbital) as the (N_e/2 + 1)th orbital, and so forth (Chai, 2012, 2017; Wu and Chai, 2015; Wu et al., 2016; Yeh et al., 2018; Chung and Chai, 2019; Deng and Chai, 2019; Seenithurai and Chai, 2019; Huang et al., 2020). For brevity, HOMO and LUMO are denoted as H and L, respectively. In addition, the orbitals with an occupation number (0.2–1.8) are regarded as the active orbitals.

As presented in Figure 4, the $\overline{f_i}$ value of n-acene obtained with the TAO-AIMD simulations at 300 K is very close to the f_i value of GS n-acene (Chai, 2012) (also see Supplementary Table 2). This implies that on average, the radical nature of n-acene obtained with the TAO-AIMD simulations at 300 K is very similar to that obtained with the GS calculation at 0 K, showing consistency with the analysis of symmetrized von Neumann entropy. For smaller n values (e.g., up to n = 5), the occupation numbers associated with all orbitals are very close to either 0 or 2. Therefore, the smaller n-acenes should exhibit nonradical nature. However, as n increases, the number of active orbitals increases and/or the active orbital occupation numbers are closer to 1, apparently showing that the larger n-acenes (n = 6–8) should exhibit increasing radical character. This clearly indicates that similar to the GS counterparts at 0 K (Hachmann et al., 2007; Chai, 2012, 2014, 2017), on average, the larger n-acenes (n = 6–8) should possess increasing radical nature in AIMD simulations at 300 K.

FIGURE 4

Figure 4. Active orbital occupation numbers (HOMO−1, HOMO, LUMO, and LUMO+1) of n-acene, obtained with the TAO-AIMD simulations at 300 K and GS calculation, calculated by TAO-LDA.

In addition, the instantaneous f_i(t) value of n-acene along each equilibrated TAO-AIMD trajectory fluctuates around an average, implying that the instantaneous radical nature of n-acene can be more pronounced than the average radical nature of n-acene in the TAO-AIMD simulations (see Figure 5 for n = 8, and Supplementary Figures 10–15 for others). According to our findings, the radical nature of the larger n-acenes (n = 6–8) can persist in AIMD simulations at 300 K. For such molecules, KS-AIMD simulations can therefore be unreliable, and AIMD simulations employing MR electronic structure methods are very likely to be computationally intractable. Accordingly, this highlights the significance of TAO-AIMD simulations for exploring the dynamical information of large molecules with radical nature.

FIGURE 5

Figure 5. Time evolution of the active orbital occupation numbers (HOMO−1, HOMO, LUMO, and LUMO+1) of 8-acene, obtained from four different equilibrated TAO-AIMD trajectories (No.1 to No.4) at 300 K, calculated by TAO-LDA. For HOMO/LUMO, the TAO-AIMD average and GS values are also shown for comparison.

4.3. IR Spectra

IR spectroscopy, which involves the interaction of IR radiation with matter, has been extensively used to explore the structure of materials nondestructively in recent years (Roggo et al., 2007; Joblin et al., 2011; Beć and Huck, 2019). Recently, the IR spectra of PAHs (e.g., n-acenes) have attracted much interest since PAHs can be responsible for the unidentified infrared emission (UIR) bands from interstellar media (Allamandola et al., 1985, 1999; Hudgins and Sandford, 1998a,b; Kim et al., 2001; Joblin et al., 2011).

On the basis of the Fermi golden rule of time-dependent perturbation theory (McQuarrie, 1976), the IR spectrum of a molecule is related to the Fourier transform of the autocorrelation function (ACF) (e.g., see Equation 15 of Brehm et al., 2020) of the dipole moment (i.e., the sum of the electronic and nuclear contributions) $\overrightarrow{\mu }$ of the molecule (Thomas et al., 2013, 2015; Dutta and Chowdhury, 2019):

$\begin{array}{ll}I(\omega )\propto {\omega }^2{\displaystyle {\int }_{-\infty }^{\infty }}\langle \overrightarrow{\mu }(0)·\overrightarrow{\mu }(t)\rangle e^{-i\omega t}dt,& (12)\end{array}$

where I(ω) is the IR intensity, and ω is the vibrational frequency. Note that a quantum correction factor βℏω/(1−e^−βℏω) (Ramírez et al., 2004; Joalland et al., 2010; Thomas et al., 2013) has been taken into account in Equation (12), where β = 1/(k_BT) and k_B is the Boltzmann constant. According to the properties of the Fourier transform (Thomas et al., 2013; Lawson Daku, 2018), I(ω) can also be expressed as

$\begin{array}{ll}I(\omega )\propto {\displaystyle {\int }_{-\infty }^{\infty }}\langle \stackrel{.}{\overrightarrow{\mu }}(0)·\stackrel{.}{\overrightarrow{\mu }}(t)\rangle e^{-i\omega t}dt,& (13)\end{array}$

which is directly proportional to the Fourier transform of the ACF of the time derivative of the dipole moment $\stackrel{.}{\overrightarrow{\mu }}$ of the molecule.

Given the $\overrightarrow{\mu }(t)$ of n-acene along an equilibrated TAO-AIMD trajectory (obtained with Q-Chem 4.4, Shao, 2015), we compute the IR spectrum of n-acene using the TRAVIS program package (Brehm and Kirchner, 2011; Thomas et al., 2013, 2015; Brehm et al., 2020). To reduce the numerical noise, the IR spectrum has been smoothed using a window function applied in the time domain [specifically, each term of the ACF is multiplied by a Gaussian function $\mathrm{\text{exp}}(-\frac{\sigma t^2}{2{\tau }^2})$, where σ = 10 (as suggested in Gaigeot et al., 2005; Gaigeot, 2010; Vitale et al., 2015 for gas phase simulations) and τ is the total time duration of the equilibrated TAO-AIMD trajectory (≈ 20.3 ps)]. The reported IR spectrum of n-acene is an average over four different equilibrated TAO-AIMD trajectories (with the average temperature being 300 ± 1 K).

On the other hand, NMA (i.e., normal mode analysis) is a commonly adopted approach to compute the vibrational frequencies and intensities of GS (Harris and Bertolucci, 1978; Wilson et al., 1980; Gaigeot et al., 2007) and excited-state (ES) (Liu and Liang, 2011a,b) molecules in the harmonic approximation at 0 K. For comparison purposes, we also compute the IR spectra of GS n-acenes using NMA. To perform a GS-NMA, the computation of nuclear second derivatives of energy (i.e., the nuclear Hessian) at the GS molecular geometry is required. Since analytical nuclear Hessians for TAO-DFT are not yet available in Q-Chem, numerical nuclear Hessians are computed using finite differences of analytical nuclear gradients (with a step size of 0.001 Å, i.e., the default setting of Q-Chem) for all the GS-NMA performed in this work.

As shown in Figures 6–9 (for comparison purposes, also see Supplementary Figures 16–19 for the TAO-AIMD results obtained with different values of σ), the IR spectra of the smaller n-acenes (n = 2–5), obtained with the TAO-AIMD simulations at 300 K and GS-NMA, are in qualitative agreement with the available experimental IR spectra (Hudgins and Sandford, 1998a,b; Boersma et al., 2014; Bauschlicher et al., 2018; NIST mass spectrometry data center, 2020; Mattioda et al., unpublished). Note that a number of experimental IR bands (e.g., those in the 1,700–2,600 cm⁻¹ range for 2-acene and 3-acene and those in the 1,700–2,000 cm⁻¹ range for 4-acene and 5-acene) are completely missing from the corresponding IR spectra obtained with the commonly used GS-NMA approach, indicating that anharmonic effects on the IR spectra cannot be completely ignored. By contrast, these IR bands can be obtained with the TAO-AIMD simulations at 300 K, going beyond the harmonic approximation of GS-NMA. The overall discrepancies between theoretical and experimental results for the IR spectra can be attributed to several factors, such as the purity of samples and environmental factors (temperature, background noise level, etc.) in the experiments and the approximate nature of the theoretical methods adopted.

FIGURE 6

Figure 6. IR spectra of 2-acene, obtained with the TAO-AIMD simulations at 300 K and GS-NMA, calculated by TAO-LDA. Experimental (EXP) data (NIST mass spectrometry data center, 2020) are included for comparison. The IR spectra are normalized to have a maximum intensity of one and, for clarity, are vertically offset from each other by the same value. Subfigures show the IR spectra in the 1,000–2,600 cm⁻¹ range.

FIGURE 7

Figure 7. IR spectra of 3-acene, obtained with the TAO-AIMD simulations at 300 K and GS-NMA, calculated by TAO-LDA. Experimental (EXP) data NIST mass spectrometry data center (2020) are included for comparison. The IR spectra are normalized to have a maximum intensity of one and, for clarity, are vertically offset from each other by the same value. Subfigures show the IR spectra in the 1,000–2,600 cm⁻¹ range.

FIGURE 8

Figure 8. IR spectra of 4-acene, obtained with the TAO-AIMD simulations at 300 K and GS-NMA, calculated by TAO-LDA. Experimental (EXP) data (Hudgins and Sandford, 1998a; Boersma et al., 2014; Bauschlicher et al., 2018; Mattioda et al., unpublished) are included for comparison. The IR spectra are normalized to have a maximum intensity of one and, for clarity, are vertically offset from each other by the same value. Subfigures show the IR spectra in the 850–2,000 cm⁻¹ range.

FIGURE 9

Figure 9. IR spectra of 5-acene, obtained with the TAO-AIMD simulations at 300 K and GS-NMA, calculated by TAO-LDA. Experimental (EXP) data (Hudgins and Sandford, 1998b; Boersma et al., 2014; Bauschlicher et al., 2018; Mattioda et al., unpublished) are included for comparison. The IR spectra are normalized to have a maximum intensity of one and, for clarity, are vertically offset from each other by the same value. Subfigures show the IR spectra in the 850–2,000 cm⁻¹ range.

To our knowledge, there are no experimental IR spectra for 6-acene, 7-acene, and 8-acene. It remains highly challenging to synthesize and isolate the larger n-acenes (n = 6–8), possibly due to their radical character (Hachmann et al., 2007; Chai, 2012, 2014, 2017; Rivero et al., 2013; Wu and Chai, 2015; Fosso-Tande et al., 2016). In Figures 10–12 (for comparison purposes, also see Supplementary Figures 20–22 for the TAO-AIMD results obtained with different values of σ), we thus only compare the IR spectra of the larger n-acenes (n = 6–8), obtained with the TAO-AIMD simulations at 300 K and GS-NMA. As shown, both theoretical results are in reasonably good agreement, except for the 1,700–2,000 cm⁻¹ range, which may be attributed to anharmonic effects on the IR spectra [cf., the IR spectra of the smaller n-acenes (n = 2–5)].

FIGURE 10

Figure 10. IR spectra of 6-acene, obtained with the TAO-AIMD simulations at 300 K and GS-NMA, calculated by TAO-LDA. The IR spectra are normalized to have a maximum intensity of one and, for clarity, are vertically offset from each other by the same value. Subfigures show the IR spectra in the 1,000–2,000 cm⁻¹ range.

FIGURE 11

Figure 11. IR spectra of 7-acene, obtained with the TAO-AIMD simulations at 300 K and GS-NMA, calculated by TAO-LDA. The IR spectra are normalized to have a maximum intensity of one and, for clarity, are vertically offset from each other by the same value. Subfigures show the IR spectra in the 1,000–2,000 cm⁻¹ range.

FIGURE 12

Figure 12. IR spectra of 8-acene, obtained with the TAO-AIMD simulations at 300 K and GS-NMA, calculated by TAO-LDA. The IR spectra are normalized to have a maximum intensity of one and, for clarity, are vertically offset from each other by the same value. Subfigures show the IR spectra in the 1,000–2,000 cm⁻¹ range.

5. Conclusions

In conclusion, we have proposed TAO-AIMD (i.e., TAO-DFT-based AIMD) for the study of the equilibrium thermodynamic and dynamical properties of nanosystems with a radical nature at finite temperatures. To highlight some of the present capabilities of TAO-AIMD, we have performed TAO-AIMD simulations to investigate the instantaneous/average radical nature and IR spectra of n-acenes (n = 2–8) at 300 K. According to the TAO-AIMD simulations, on average, the smaller n-acenes (up to n = 5) possess a nonradical nature, and the larger n-acenes (n = 6–8) possess increasing radical nature, showing remarkable similarities to the GS counterparts at 0 K. Besides, the IR spectra of n-acenes obtained with the TAO-AIMD simulations are in qualitative agreement with the existing experimental data.

For GS molecules with radical nature [e.g., the larger n-acenes (n = 6–8)], on average, the radical nature of the molecules can persist in AIMD simulations at finite temperatures. For these molecules, conventional KS-AIMD simulations can therefore be unreliable, and AIMD simulations employing MR electronic structure methods are very likely to be computationally intractable for most molecules. It is thus certainly justified to perform TAO-AIMD simulations for exploring the dynamical information of these molecules.

While only TAO-AIMD (or more specifically, TAO-BOMD) is presented and discussed in this work, it is also possible to combine TAO-DFT with the popular Car-Parrinello MD (CPMD) (Car and Parrinello, 1985) method (i.e., an approximation of the BOMD method) for improved computational efficiency. In addition, it is worth mentioning that a brilliant simulation method combining the advantages of AIMD (for accuracy) and classical MD (for efficiency) has been developed, i.e., the hybrid QM/MM (quantum mechanics/molecular mechanics) method (Warshel and Levitt, 1976), where a small portion (e.g., the reactive portion) of a system is treated with QM and the remaining portion is treated with MM. The QM/MM method has been widely used for the study of very large systems where AIMD simulations are prohibitively expensive (van der Kamp and Mulholland, 2013). To further improve the efficiency of TAO-AIMD, the TAO-DFT-based QM/MM method is thus expected to be useful for the simulations of very large systems (e.g., biomolecules).

While TAO-AIMD is computationally efficient, it could be a very promising approach for the study of the equilibrium thermodynamic and dynamical properties of nanosystems with a radical or non-radical nature at finite temperatures. Nevertheless, a few limitations remain due to a number of assumptions made in TAO-AIMD. For example, the BO approximation is assumed to be valid, and the motion of the nuclei in a system is assumed to evolve only on the GS potential energy surface, wherein non-adiabatic effects are completely ignored. Besides, it is assumed that nuclear quantum effects can be neglected: the nuclear motion is assumed to be governed by the classical equations of motion, rather than the time-dependent nuclear Schrödinger equation. Moreover, it is assumed that TAO-DFT employing the approximate XC and θ-dependent energy functionals (with a system-independent fictitious temperature θ) can provide the exact GS potential energy surface. In addition, similar to KS-AIMD, the real electronic temperature is zero in TAO-AIMD, and it is thus assumed that TAO-AIMD remains applicable for systems at nonzero electronic temperatures. For a system at room temperature, TAO-AIMD should remain applicable (Gaigeot, 2010; Ramírez-Solís et al., 2011; Vitale et al., 2015). However, for a system at a considerably high electronic temperature (e.g., warm dense matter), TAO-AIMD may no longer be applicable (Rüter and Redmer, 2014; Karasiev et al., 2018; Bonitz et al., 2020), wherein AIMD based on the finite-electronic-temperature extension of TAO-DFT will be needed. To lift these limitations, we plan to investigate along some of these lines, and results may be reported elsewhere.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

J-DC conceived the project. SL performed the calculations. SL and J-DC designed the project, performed the data analysis, and wrote the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by the Ministry of Science and Technology of Taiwan (Grant No. MOST107-2628-M-002-005-MY3), National Taiwan University (Grant No. NTU-CDP-105R7818), and the NCTS of Taiwan.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank Dr. Martin Brehm for useful discussions on the evaluation of the IR spectra.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fchem.2020.589432/full#supplementary-material

References

Allamandola, L. J., Hudgins, D. M., and Sandford, S. A. (1999). Modeling the unidentified infrared emission with combinations of polycyclic aromatic hydrocarbons. Astrophys. J 511, L115–L119. doi: 10.1086/311843

PubMed Abstract | CrossRef Full Text | Google Scholar

Allamandola, L. J., Tielens, A. G. G. M., and Barker, J. R. (1985). Polycyclic aromatic hydrocarbons and the unidentified infrared emission bands: auto exhaust along the milky way. Astrophys. J. 290, L25–L28. doi: 10.1086/184435

CrossRef Full Text | Google Scholar

Andersson, K., Malmqvist, P.-Å., and Roos, B. O. (1992). Second-order perturbation theory with a complete active space self-consistent field reference function. J. Chem. Phys. 96, 1218–1226. doi: 10.1063/1.462209

CrossRef Full Text | Google Scholar

Battaglia, S., Faginas-Lago, N., Andrae, D., Evangelisti, S., and Leininger, T. (2017). Increasing radical character of large [n]cyclacenes unveiled by wave function theory. J. Phys. Chem. A 121, 3746–3756. doi: 10.1021/acs.jpca.7b00123

PubMed Abstract | CrossRef Full Text | Google Scholar

Bauschlicher, C. W. Jr., Ricca, A., Boersma, C., and Allamandola, L. J. (2018). The NASA Ames PAH IR spectroscopic database: computational version 3.00 with updated content and the introduction of multiple scaling factors. ApJS 234:32. doi: 10.3847/1538-4365/aaa019

CrossRef Full Text | Google Scholar

Beć, K. B., and Huck, C. W. (2019). Breakthrough potential in near-infrared spectroscopy: spectra simulation. A review of recent developments. Front. Chem. 7:48. doi: 10.3389/fchem.2019.00048

PubMed Abstract | CrossRef Full Text | Google Scholar

Boersma, C., Bauschlicher, C. W. Jr., Ricca, A., Mattioda, A. L., Cami, J., Peeters, E., et al. (2014). The NASA Ames PAH IR spectroscopic database version 2.00: updated content, web site, and on(off)line tools. ApJS 211:8. doi: 10.1088/0067-0049/211/1/8

CrossRef Full Text | Google Scholar

Bonitz, M., Dornheim, T., Moldabekov, Zh. A., Zhang, S., Hamann, P., et al. (2020). Ab initio simulation of warm dense matter. Phys. Plasmas 27:042710. doi: 10.1063/1.5143225

CrossRef Full Text | Google Scholar

Born, M., and Oppenheimer, J. R. (1927). Zur Quantentheorie der Molekeln. Ann. Phys. 84, 457–484. doi: 10.1002/andp.19273892002

CrossRef Full Text | Google Scholar

Brehm, M., and Kirchner, B. (2011). TRAVIS – a free analyzer and visualizer for monte carlo and molecular dynamics trajectories. J. Chem. Inf. Model. 51, 2007–2023. doi: 10.1021/ci200217w

PubMed Abstract | CrossRef Full Text | Google Scholar

Brehm, M., Thomas, M., Gehrke, S., and Kirchner, B. (2020). TRAVIS – a free analyzer for trajectories from molecular simulation. J. Chem. Phys. 152:164105. doi: 10.1063/5.0005078

PubMed Abstract | CrossRef Full Text | Google Scholar

Car, R., and Parrinello, M. (1985). Unified approach for molecular dynamics and density-functional theory. Phys. Rev. Lett. 55, 2471–2474. doi: 10.1103/PhysRevLett.55.2471

PubMed Abstract | CrossRef Full Text | Google Scholar

Chai, J.-D. (2012). Density functional theory with fractional orbital occupations. J. Chem. Phys. 136:154104. doi: 10.1063/1.3703894

PubMed Abstract | CrossRef Full Text | Google Scholar

Chai, J.-D. (2014). Thermally-assisted-occupation density functional theory with generalized-gradient approximations. J. Chem. Phys. 140:18A521. doi: 10.1063/1.4867532

PubMed Abstract | CrossRef Full Text | Google Scholar

Chai, J.-D. (2017). Role of exact exchange in thermally-assisted-occupation density functional theory: a proposal of new hybrid schemes. J. Chem. Phys. 146:044102. doi: 10.1063/1.4974163

PubMed Abstract | CrossRef Full Text | Google Scholar

Chai, J.-D., Stroud, D., Hafner, J., and Kresse, G. (2003). Dynamic structure factor of liquid and amorphous Ge from ab initio simulations. Phys. Rev. B 67:104205. doi: 10.1103/PhysRevB.67.104205

CrossRef Full Text | Google Scholar

Chung, J.-H., and Chai, J.-D. (2019). Electronic properties of Möbius cyclacenes studied by thermally-assisted-occupation density functional theory. Sci. Rep. 9:2907. doi: 10.1038/s41598-019-39524-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Cohen, A. J., Mori-Sánchez, P., and Yang, W. (2008). Insights into current limitations of density functional theory. Science 321, 792–794. doi: 10.1126/science.1158722

PubMed Abstract | CrossRef Full Text | Google Scholar

Cohen, A. J., Mori-Sánchez, P., and Yang, W. (2012). Challenges for density functional theory. Chem. Rev. 112, 289–320. doi: 10.1021/cr200107z

PubMed Abstract | CrossRef Full Text | Google Scholar

Deng, Q., and Chai, J.-D. (2019). Electronic properties of triangle-shaped graphene nanoflakes from TAO-DFT. ACS Omega 4, 14202–14210. doi: 10.1021/acsomega.9b01259

PubMed Abstract | CrossRef Full Text | Google Scholar

Dutta, B., and Chowdhury, J. (2019). Existence of dimeric hydroxylamine-O-sulfonic acid: experimental observations aided by ab initio, DFT, Car-Parrinello and Born-Oppenheimer on the fly dynamics. Chem. Phys. Lett. 732:136645. doi: 10.1016/j.cplett.2019.136645

CrossRef Full Text | Google Scholar

Fosso-Tande, J., Nguyen, T.-S., Gidofalvi, G., and DePrince III, A. E. (2016). Large-scale variational two-electron reduced-density-matrix-driven complete active space self-consistent field methods. J. Chem. Theory Comput. 12, 2260–2271. doi: 10.1021/acs.jctc.6b00190

PubMed Abstract | CrossRef Full Text | Google Scholar

Gaigeot, M.-P. (2010). Theoretical spectroscopy of floppy peptides at room temperature. A DFTMD perspective: gas and aqueous phase. Phys. Chem. Chem. Phys. 12, 3336–3359. doi: 10.1039/B924048A

PubMed Abstract | CrossRef Full Text | Google Scholar

Gaigeot, M.-P., Martinez, M., and Vuilleumier, R. (2007). Infrared spectroscopy in the gas and liquid phase from first principle molecular dynamics simulations: application to small peptides. Mol. Phys. 105, 2857–2878. doi: 10.1080/00268970701724974

CrossRef Full Text | Google Scholar

Gaigeot, M.-P., and Sprik, M. (2003). Ab initio molecular dynamics computation of the infrared spectrum of aqueous uracil. J. Phys. Chem. B 107, 10344–10358. doi: 10.1021/jp034788u

CrossRef Full Text | Google Scholar

Gaigeot, M. P., Vuilleumier, R., Sprik, M., and Borgis, D. (2005). Infrared spectroscopy of N-Methylacetamide revisited by ab Initio molecular dynamics simulations. J. Chem. Theory Comput. 1, 772–789. doi: 10.1021/ct050029z

PubMed Abstract | CrossRef Full Text | Google Scholar

Gidofalvi, G., and Mazziotti, D. A. (2008). Active-space two-electron reduced-density-matrix method: complete active-space calculations without diagonalization of the N-electron hamiltonian. J. Chem. Phys. 129:134108. doi: 10.1063/1.2983652

PubMed Abstract | CrossRef Full Text | Google Scholar

Gryn'ova, G., Coote, M. L., and Corminboeuf, C. (2015). Theory and practice of uncommon molecular electronic configurations. WIREs Comput. Mol. Sci. 5, 440–459. doi: 10.1002/wcms.1233

PubMed Abstract | CrossRef Full Text | Google Scholar

Hachmann, J., Dorando, J. J., Aviles, M., and Chan, G. K. L. (2007). The radical character of the acenes: a density matrix renormalization group study. J. Chem. Phys. 127:134309. doi: 10.1063/1.2768362

PubMed Abstract | CrossRef Full Text | Google Scholar

Hanson-Heine, M. W. D. (2020). Static correlation in vibrational frequencies studied using thermally-assisted-occupation density functional theory. Chem. Phys. Lett. 739:137012. doi: 10.1016/j.cplett.2019.137012

CrossRef Full Text | Google Scholar

Hanson-Heine, M. W. D., and Hirst, J. D. (2020). Möbius and Hückel cyclacenes with Dewar and Ladenburg defects. J. Phys. Chem. A 124, 5408–5414. doi: 10.1021/acs.jpca.0c04137

PubMed Abstract | CrossRef Full Text | Google Scholar

Hanson-Heine, M. W. D., Rogers, D. M., Woodward, S., and Hirst, J. D. (2020). Dewar benzenoids discovered in carbon nanobelts. J. Phys. Chem. Lett. 11, 3769–3772. doi: 10.1021/acs.jpclett.0c01027

PubMed Abstract | CrossRef Full Text | Google Scholar

Harris, D. C., Bertolucci, M. D. (1978). Symmetry and Spectroscopy: An Introduction to Vibrational and Electronic Spectroscopy. New York: Oxford University Press.

Google Scholar

Hollingsworth, S. A., and Dror, R. O. (2018). Molecular dynamics simulation for all. Neuron 99, 1129–1143. doi: 10.1016/j.neuron.2018.08.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, H.-J., Seenithurai, S., and Chai, J.-D. (2020). TAO-DFT study on the electronic properties of diamond-shaped graphene nanoflakes. Nanomaterials 10:1236. doi: 10.3390/nano10061236

PubMed Abstract | CrossRef Full Text | Google Scholar

Hudgins, D. M., and Sandford, S. A. (1998a). Infrared spectroscopy of matrix isolated polycyclic aromatic hydrocarbons. 1. PAHs containing two to four rings. J. Phys. Chem. A 102, 329–343. doi: 10.1021/jp9834816

CrossRef Full Text | Google Scholar

Hudgins, D. M., and Sandford, S. A. (1998b). Infrared spectroscopy of matrix isolated polycyclic aromatic hydrocarbons. 2. PAHs containing five or more rings. J. Phys. Chem. A 102, 344–352. doi: 10.1021/jp983482y

CrossRef Full Text | Google Scholar

Jensen, F. (2007). Introduction to Computational Chemistry. New York, NY: Wiley.

PubMed Abstract | Google Scholar

Joalland, B., Rapacioli, M., Simon, A., Joblin, C., Marsden, C. J., and Spiegelman, F. (2010). Molecular dynamics simulations of anharmonic infrared spectra of [SiPAH]⁺ π-complexes. J. Phys. Chem. A 114, 5846–5854. doi: 10.1021/jp911526n

PubMed Abstract | CrossRef Full Text | Google Scholar

Joblin, C., Tielens, A., and Oomens, J. (2011). Laboratory infrared spectroscopy of PAHs. EAS Publ. Ser. 46, 61–73. doi: 10.1051/eas/1146007

CrossRef Full Text | Google Scholar

Karasiev, V. V., Dufty, J. W., and Trickey, S. B. (2018). Nonempirical semilocal free-energy density functional for matter under extreme conditions. Phys. Rev. Lett. 120:076401. doi: 10.1103/PhysRevLett.120.076401

PubMed Abstract | CrossRef Full Text | Google Scholar

Karplus, M., and Petsko, G. A. (1990). Molecular dynamics simulations in biology. Nature 347, 631–639. doi: 10.1038/347631a0

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, H. S., Wagner, D. R., and Saykally, R. J. (2001). Single photon infrared emission spectroscopy of the gas phase pyrene cation: support for a polycyclic aromatic hydrocarbon origin of the unidentified infrared emission bands. Phys. Rev. Lett. 86, 5691–5694. doi: 10.1103/PhysRevLett.86.5691

PubMed Abstract | CrossRef Full Text | Google Scholar

Kohn, W., and Sham, L. J. (1965). Self-consistent equations including exchange and correlation effects. Phys. Rev. 140, A1133–A1138. doi: 10.1103/PhysRev.140.A1133

CrossRef Full Text | Google Scholar

Kresse, G., and Hafner, J. (1993). Ab initio molecular dynamics for liquid metals. Phys. Rev. B 47, 558–561. doi: 10.1103/PhysRevB.47.558

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuo, I. F. W., and Mundy, C. J. (2004). An ab Initio molecular dynamics study of the aqueous liquid-vapor interface. Science 303, 658–660. doi: 10.1126/science.1092787

PubMed Abstract | CrossRef Full Text | Google Scholar

Lawson Daku, L. M. (2018). Spin-state dependence of the structural and vibrational properties of solvated iron(ii) polypyridyl complexes from AIMD simulations: aqueous [Fe(bpy)3]Cl2, a case study. Phys. Chem. Chem. Phys. 20, 6236–6253. doi: 10.1039/C7CP07862E

PubMed Abstract | CrossRef Full Text | Google Scholar

Levitt, M., and Lifson, S. (1969). Refinement of protein conformations using a macromolecular energy minimization procedure. J. Mol. Biol. 46, 269–279. doi: 10.1016/0022-2836(69)90421-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Lifson, S., and Warshel, A. (1968). Consistent force field for calculations of conformations, vibrational spectra, and enthalpies of cycloalkane and n-alkane molecules. J. Chem. Phys. 49, 5116–5129. doi: 10.1063/1.1670007

CrossRef Full Text | Google Scholar

Lin, C.-Y., Hui, K., Chung, J.-H., and Chai, J.-D. (2017). Self-consistent determination of the fictitious temperature in thermally-assisted-occupation density functional theory. RSC Adv. 7, 50496–50507. doi: 10.1039/C7RA10241K

CrossRef Full Text | Google Scholar

Liu, J., and Liang, W. Z. (2011a). Analytical Hessian of electronic excited states in time-dependent density functional theory with Tamm-Dancoff approximation. J. Chem. Phys. 135:014113. doi: 10.1063/1.3605504

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, J., and Liang, W. Z. (2011b). Analytical approach for the excited-state Hessian in time-dependent density functional theory: formalism, implementation, and performance. J. Chem. Phys. 135:184111. doi: 10.1063/1.3659312

PubMed Abstract | CrossRef Full Text | Google Scholar

Löwdin, P.-O., and Shull, H. (1956). Natural orbitals in the quantum theory of two-electron systems. Phys. Rev. 101, 1730–1739. doi: 10.1103/PhysRev.101.1730

CrossRef Full Text | Google Scholar

Manassir, M., and Pakiari, A. H. (2020). Total non-Lewis structures: an application to predict the stability and reactivity of linear and angular polyacenes. J. Mol. Graph. Model. 99:107643. doi: 10.1016/j.jmgm.2020.107643

CrossRef Full Text | Google Scholar

Martyna, G. J., Klein, M. L., and Tuckerman, M. (1992). Nosé-Hoover chains: the canonical ensemble via continuous dynamics. J. Chem. Phys. 97, 2635–2643. doi: 10.1063/1.463940

CrossRef Full Text | Google Scholar

Marx, D., and Hutter, J. (2009). Ab Initio Molecular Dynamics: Basic Theory and Advanced Methods. Cambridge: Cambridge University Press.

Google Scholar

McQuarrie, D. A. (1976). Statistical Mechanics. New York, NY: Harper and Row Publishers.

Mullinax, J. W., Maradzike, E., Koulias, L. N., Mostafanejad, M., Epifanovsky, E., Gidofalvi, G., and DePrince III, A. E. (2019). Heterogeneous CPU + GPU algorithm for variational two-electron reduced-density matrix-driven complete active-space self-consistent field theory. J. Chem. Theory Comput. 15, 6164–6178. doi: 10.1021/acs.jctc.9b00768

PubMed Abstract | CrossRef Full Text | Google Scholar

NIST mass spectrometry data center, Wallace, W. E. director, “Infrared Spectra” in NIST chemistry webbook, NIST standard reference database number 69, eds. Linstrom, P. J., and Mallard, W. G. National Institute of Standards Technology Gaithersburg MD, 20899, https://doi.org/10.18434/T4D303, (retrieved June 17, 2020).

Pelzer, K., Greenman, L., Gidofalvi, G., and Mazziotti, D. A. (2011). Strong correlation in acene sheets from the active-space variational two-electron reduced density matrix method: effects of symmetry and size. J. Phys. Chem. A 115, 5632–5640. doi: 10.1021/jp2017192

PubMed Abstract | CrossRef Full Text | Google Scholar

Putrino, A., and Parrinello, M. (2002). Anharmonic Raman spectra in high-pressure ice from ab initio simulations. Phys. Rev. Lett. 88:176401. doi: 10.1103/PhysRevLett.88.176401

PubMed Abstract | CrossRef Full Text | Google Scholar

Ramírez, R., López-Ciudad, T., Kumar P, P., and Marx, D. (2004). Quantum corrections to classical time-correlation functions: hydrogen bonding and anharmonic floppy modes. J. Chem. Phys. 121, 3973–3983. doi: 10.1063/1.1774986

PubMed Abstract | CrossRef Full Text | Google Scholar

Ramírez-Solís, A., Jolibois, F., and Maron, L. (2011). Born-Oppenheimer DFT molecular dynamics studies of S₂O₂: non-harmonic effects on the lowest energy isomers. Chem. Phys. Lett. 510, 21–26. doi: 10.1016/j.cplett.2011.04.091

CrossRef Full Text | Google Scholar

Rivero, P., Jiménez-Hoyos, C. A., and Scuseria, G. E. (2013). Entanglement and polyradical character of polycyclic aromatic hydrocarbons predicted by projected Hartree-Fock theory. J. Phys. Chem. B 117, 12750–12758. doi: 10.1021/jp401478v

PubMed Abstract | CrossRef Full Text | Google Scholar

Roggo, Y., Chalus, P., Maurer, L., Lema-Martinez, C., Edmond, A., and Jent, N. (2007). A review of near infrared spectroscopy and chemometrics in pharmaceutical technologies. J. Pharm. Biomed. Anal. 44, 683–700. doi: 10.1016/j.jpba.2007.03.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Rüter, H. R., and Redmer, R. (2014). Ab initio simulations for the ion-ion structure factor of warm dense aluminum. Phys. Rev. Lett. 112:145007. doi: 10.1103/PhysRevLett.112.145007

PubMed Abstract | CrossRef Full Text | Google Scholar

Seenithurai, S., and Chai, J.-D. (2016). Effect of Li adsorption on the electronic and hydrogen storage properties of acenes: a dispersion-corrected TAO-DFT study. Sci. Rep. 6:33081. doi: 10.1038/srep33081

PubMed Abstract | CrossRef Full Text | Google Scholar

Seenithurai, S., and Chai, J.-D. (2017). Effect of Li termination on the electronic and hydrogen storage properties of linear carbon chains: a TAO-DFT study. Sci. Rep. 7:4966. doi: 10.1038/s41598-017-05202-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Seenithurai, S., and Chai, J.-D. (2018). Electronic and hydrogen storage properties of Li-terminated linear boron chains studied by TAO-DFT. Sci. Rep. 8:13538. doi: 10.1038/s41598-018-31947-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Seenithurai, S., and Chai, J.-D. (2019). Electronic properties of linear and cyclic boron nanoribbons from thermally-assisted-occupation density functional theory. Sci. Rep. 9:12139. doi: 10.1038/s41598-019-48560-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Seenithurai, S., and Chai, J.-D. (2020). TAO-DFT investigation of electronic properties of linear and cyclic carbon chains. Sci. Rep. 10:13133. doi: 10.1038/s41598-020-70023-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Shao, Y., Gan, Z., Epifanovsky, E., Gilbert, T. B., Wormit, M., Kussmann, J., et al. (2015). Advances in molecular quantum chemistry contained in the Q-Chem 4 program package. Mol. Phys. 113, 184–215. doi: 10.1080/00268976.2014.952696

CrossRef Full Text | Google Scholar

Silvestrelli, P. L., Bernasconi, M., and Parrinello, M. (1997). Ab initio infrared spectrum of liquid water. Chem. Phys. Lett. 277, 478–482. doi: 10.1016/S0009-2614(97)00930-5

CrossRef Full Text | Google Scholar

Sprik, M., Hutter, J., and Parrinello, M. (1996). Ab initio molecular dynamics simulation of liquid water: comparison of three gradient–corrected density functionals. J. Chem. Phys. 105, 1142–1152. doi: 10.1063/1.471957

CrossRef Full Text | Google Scholar

Thomas, M., Brehm, M., Fligg, R., Vöhringer, P., and Kirchner, B. (2013). Computing vibrational spectra from ab initio molecular dynamics. Phys. Chem. Chem. Phys. 15, 6608–6622. doi: 10.1039/C3CP44302G

PubMed Abstract | CrossRef Full Text | Google Scholar

Thomas, M., Brehm, M., and Kirchner, B. (2015). Voronoi dipole moments for the simulation of bulk phase vibrational spectra. Phys. Chem. Chem. Phys. 17, 3207–3213. doi: 10.1039/C4CP05272B

PubMed Abstract | CrossRef Full Text | Google Scholar

Tuckerman, M. E. (2002). Ab initio molecular dynamics: basic concepts, current trends and novel applications. J. Phys. 14, R1297–R1355. doi: 10.1088/0953-8984/14/50/202

CrossRef Full Text | Google Scholar

van der Kamp, M. W., and Mulholland, A. J. (2013). Combined quantum mechanics/molecular mechanics (QM/MM) methods in computational enzymology. Biochemistry 52, 2708–2728. doi: 10.1021/bi400215w

PubMed Abstract | CrossRef Full Text | Google Scholar

Vitale, V., Dziedzic, J., Dubois, S. M.-M., Fangohr, H., and Skylaris, C.-K. (2015). Anharmonic infrared spectroscopy through the Fourier transform of time correlation function formalism in ONETEP. J. Chem. Theory Comput. 11, 3321–3332. doi: 10.1021/acs.jctc.5b00391

PubMed Abstract | CrossRef Full Text | Google Scholar

Warshel, A., and Levitt, M. (1976). Theoretical studies of enzymic reactions: dielectric, electrostatic and steric stabilization of the carbonium ion in the reaction of lysozyme. J. Mol. Biol. 103, 227–249. doi: 10.1016/0022-2836(76)90311-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Wilson, E. B. Jr, Decius, J. C., and Cross, P. C. (1980). Molecular Vibrations: The Theory of Infrared and Raman Vibrational Spectra. New York, NY: Dover Publications, Inc.

Google Scholar

Wu, C.-S., and Chai, J.-D. (2015). Electronic properties of zigzag graphene nanoribbons studied by TAO-DFT. J. Chem. Theory Comput. 11, 2003–2011. doi: 10.1021/ct500999m

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, C.-S., Lee, P.-Y., and Chai, J.-D. (2016). Electronic properties of cyclacenes from TAO-DFT. Sci. Rep. 6:37249. doi: 10.1038/srep37249

PubMed Abstract | CrossRef Full Text | Google Scholar

Xuan, F., Chai, J.-D., and Su, H. (2019). Local density approximation for the short-range exchange free energy functional. ACS Omega 4, 7675–7683. doi: 10.1021/acsomega.9b00303

PubMed Abstract | CrossRef Full Text | Google Scholar

Yeh, C.-N., and Chai, J.-D. (2016). Role of Kekulé and non-Kekulé structures in the radical character of alternant polycyclic aromatic hydrocarbons: a TAO-DFT study. Sci. Rep. 6:30562. doi: 10.1038/srep30562

CrossRef Full Text | Google Scholar

Yeh, C.-N., Wu, C., Su, H., and Chai, J.-D. (2018). Electronic properties of the coronene series from thermally-assisted-occupation density functional theory. RSC Adv. 8, 34350–34358. doi: 10.1039/C8RA01336E

CrossRef Full Text | Google Scholar

Yeh, S.-H., Manjanath, A., Cheng, Y.-C., Chai, J.-D., and Hsu, C.-P. (2020). Excitation energies from thermally-assisted-occupation density functional theory: theory and computational implementation. J. Chem. Phys. 153:084120. doi: 10.1063/1.5140243

PubMed Abstract | CrossRef Full Text | Google Scholar