Quantum Vibronic Effects on the Electronic Properties of Molecular Crystals

We present a study of molecular crystals, focused on the effect of nuclear quantum motion and anharmonicity on their electronic properties. We consider a system composed of relatively rigid molecules, a diamondoid crystal, and one composed of floppier molecules, NAI-DMAC, a thermally activated delayed fluorescence compound. We compute fundamental electronic gaps at the density functional theory (DFT) level of theory, with the Perdew–Burke–Erzenhof (PBE) and strongly constrained and approximately normed (SCAN) functionals, by coupling first-principles molecular dynamics with a nuclear quantum thermostat. We find a sizable zero-point renormalization (ZPR) of the band gaps, which is much larger in the case of diamondoids (0.6 eV) than for NAI-DMAC (0.22 eV). We show that the frozen phonon (FP) approximation, which neglects intermolecular anharmonic effects, leads to a large error (∼50%) in the calculation of the band gap ZPR. Instead, when using a stochastic method, we obtain results in good agreement with those of our quantum simulations for the diamondoid crystal. However, the agreement is worse for NAI-DMAC where intramolecular anharmonicities contribute to the ZPR. Our results highlight the importance of accurately including nuclear and anharmonic quantum effects to predict the electronic properties of molecular crystals.


S1.1 Theory
Within the Born-Oppenheimer (BO) approximation, the Schrödinger equation for the quantum nuclear motion of a system with N nuclei can be expressed as, 1 H|χ k (R)⟩ = − 1 2 where M I and R = (R 1x , R 1y , ..., R N z ) are the mass of the I-th nucleus and the Cartesian position vector of all nuclei with respect to a chosen origin, respectively. V(R), and |χ k (R)⟩ represent the 3N -dimensional adiabatic potential energy surface, and the wave function for the k-th nuclear state, respectively. Since, V(R) and |χ k (R)⟩ are indifferent to the choice of the origin, for the convenience of discussion, we set the origin at the coordinates of the equilibrium geometry where V(R) is minimum. We consider only the Γ-point vibrations within the simulation cell.
We consider an electronic observable, e.g., electronic eigenenergies or total energy, O(R) = ⟨ψ(r; R)|Ô(r; R)|ψ(r; R)⟩, that depends on one electronic state. Here r and |ψ(r; R)⟩ are the electronic coordinates and wave function, respectively. We note that as a consequence of BO approximation, the electronic wave function has only parametric dependence on the nuclear coordinates, R. When the system is at equilibrium at temperature T , the effect of electron-phonon interaction on the electronic property O(R) can be included by performing an ensemble average over all adiabatic nuclear states k, Here, denotes the probability of finding the system with the nuclear coordinates within R and R + dR. The partition function Q is defined as, where k B is the Boltzmann constant.
A molecular dynamics simulation with either a path-integral approach 2,3 or quantum thermostat approach 4-6 utilizes Eq. S2 to compute the electron-phonon renormalized electronic properties. However, being computationally expensive, such simulations are rarely adopted for computing electron-phonon renormalizations from first principles. The standard approaches in solid-state physics employ approximations to V(R), to simplify the expression in Eq. S2.
A Taylor series expansion of V(R) near the equilibrium geometry, R = 0, and subsequently, neglecting the terms beyond second order yields the potential energy with harmonic approximation (HA), with M representing a 3N × 3N diagonal matrix of nuclear masses. We note that the first order term of the Taylor expansion becomes zero by setting the condition of the minimum,

∂V(R)
∂R Iα = 0. In addition, we set the reference potential energy, V(0) = 0. The elements of the dynamical matrix, which is also known as the mass-weighted Hessian matrix is given by, where α,α ′ denote the cartesian axes x,y or z and I, J denote the indices of the nuclei.
Spectral decomposition of the dynamical matrix, D = UΩ 2 U T , returns a unitary matrix, U, and a 3N × 3N diagonal matrix of normal-mode frequencies, Ω, with diagonal elements: ω 1 , ω 2 , ..., ω 3N . The unitary matrix, U, defines the cartesian to normal mode transformations, and back transformation to cartesian from normal modes, with X representing the matrix of 3N normal mode vectors.
After transformation to normal modes, the total nuclear Hamiltonian separates into 3 independent translational degrees of freedom, d r number of independent global rotational degrees of freedom, and 3N − 3 − d r number of independent vibrational degrees of freedom.
For a solid, linear isolated molecule, and a non-linear isolated molecule, the number of rotational degrees of freedom (d r ) is 0, 2, and 3, respectively. Since the translations and global rotations, which usually appear as the first 3 + d r lowest eigenvalues, do not affect the electronic properties, we focus on the vibrational Hamiltonian which becomes, The wavefunctions, |χ H ν,k ⟩, and energies, ε H ν,k of each simple harmonic oscillator is known analytically, with H k denoting the k-th order Hermite polynomial. Inserting the expression of ε H ν,k from S4 Eq. S11 into Eq. S4 yields, the partition function for the ν-th harmonic oscillator, where the Bose occupation factor is given by, Because of the separation of Hamiltonian into 3N −3−d r number of independent vibrational terms, the total vibrational wave function which is characterized by a vector of quantum numbers k = (k 3+dr+1 , k 3+dr+2 , ..., k 3N ) can be simplified by the product of the wave functions for each individual normal mode. The total partition function can be written analogously. (S14) Utilizing the expressions from Eqs.S10-S14, Eq. S2 can be rewritten as, where the harmonic probability density, W H (X, T ), reduces to a product of independent Gaussian functions, with widths related to the Bose occupation factor: We note that though Eq. S15 is a result of harmonic approximation, it does not assume any explicit dependence of electronic observable O on nuclear coordinates (R) or normal mode coordinates (X). To further simplify the expression, we expand O(X) in Taylor series, and truncate the expansion after the second order. After inserting the resulting expression into Eq.S15, we obtain the phonon renormalized electronic observable with quadratic (Q) approximation, (S19) We note that terms involving odd order of X ν or cross-coupling terms such as X ν X ν ′ , with ν ̸ = ν ′ , do not appear because harmonic density is symmetric with respect to X = 0.
For systems with strong anharmonicity, the vibrational density is no longer symmetric, and subsequently, both odd-order terms and cross-coupling terms become essential.

S1.2 Stochastic Approach
The stochastic approach employs Monte Carlo sampling of W (X, T ) and subsequently utilizes Eq. S15 to compute phonon-renormalized electronic properties. In each Monte Carlo step, a displaced coordinate in normal mode is obtained, and they showed that only a single first-principle calculation on the obtained atomic configuration, X = sσ T is sufficient to converge of ⟨O⟩ T . They also proposed that an additional first-principles calculation on the antithetic pair of the chosen atomic configuration, i.e., X = −sσ T , can improve the result, which we adopted here.

S1.3 Frozen Phonon Approach
A frozen phonon (FP) approach utilizes Eq. S19 to compute phonon-renormalized electronic properties. To compute the first and second derivatives of the observable, O ′ , O ′′ , respectively, we displaced the nuclei to X ν = +h ν and X ν = −h ν along each phonon coordinates, and subsequently, used the central difference formula.
We chose the displacements, h ν , such that the potential energy difference, δV H , is the same for all modes, assuming a parabolic dependence of potential energy on each normal mode, i.e., Throughout this work, our electronic observable (O) would be either valence and conduction band (VB) energies, E n , with n denoting the band index, or the band gap, E g .
The second derivative of n-th band energies with respect to ν-th phonon mode scaled by the phonon frequency, which appears in the Eq. S19 when O = E n , is called the electron-phonon coupling energy (EPCE) of band n with respect to phonon mode ν.
It is evident from Eq. S19, the electron-phonon renormalization of n-th band, within quadratic approximation at 0 K. EPCEs can be computed using the FP method and it describes the contributions of each mode towards the total electron-phonon renormalization of band energies.
Depending on the displacement chosen for each phonon mode, nearby eigenstates may cross each other, as we shall see during the discussion. In such cases, while calculating