{'fraction_non_alphanumeric': 0.027892561983471075, 'fraction_numerical': 0.01115702479338843, 'mean_word_length': 4.822373696872494, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Analytical infrared intensities for periodic systems with local basis sets

We report a method for the efficient evaluation of analytic infrared (IR) intensities within generalized Kohn-Sham density functional theory using Gaussian orbitals and periodic boundary conditions. A discretized form of the Berry phase is used to evaluate a periodic dipole moment and its derivatives with respect to in-phase nuclear coordinate displacements. Benchmark calculations are presented for one-dimensional chains of water molecules and poly(paraphenylenevinylene).


I. INTRODUCTION
Theoretical prediction of vibrational spectra is a valuable tool for the characterization of compounds and materials which are either not synthesized or difficult to conduct experiments on. It requires at least two types of data: positions of peaks and their intensities. In the vibrational spectroscopy of periodic systems, peak positions are associated with phonon frequencies. Calculation of phonon frequencies from analytical second derivatives have recently become available for Kohn-Sham (KS) density functional theory (DFT) and Hartree-Fock (HF) methods, with plane-waves (PW) 1,2 and Gaussian type orbitals (GTO). 3,4 On the other hand, peak intensities have received less attention, even though they are necessary for comparison of theoretical spectra with their experimental counterparts. Part of the problem is that only in the 90's was their correct theoretical expression derived. 5,6 Since then, first principle studies with KS-DFT and PWs have become common for infrared (IR) and Raman spectra of solids. 7,8,9 However, we are aware of only one work on IR intensities with GTOs. Although GTOs have many advantages, 4,10 especially for orbital dependent functionals, GTO-based techniques are mostly developed in the quantum chemistry community, and therefore, are underrepresented for periodic systems. Thus, in this work, we would like to present a simple scheme for evaluation of IR intensities for periodic systems within KS-DFT and HF with GTOs.
The first order IR intensity of a fundamental transition exciting the jth normal mode is where N A is Avogadro's number, c is the speed of light, d x−z are the Cartesian components of the dipole moment, and Q j is the normal mode coordinate. 11 Since any normal mode is a linear combination of atomic Cartesian coordinates (R i ), IR intensities are calculated from ∂d/∂R i . There are two main features arising for these derivatives under periodic boundary conditions (PBC). First, one should consider derivatives with respect to in-phase nuclear coordinate displacements where g is the index of an unit cell. Although in-phase vibrations do not represent the whole phonon spectrum of the periodic system, only these vibrations are IR active. 12 Second, under PBC, the dipole moment per unit cell cannot be straightforwardly expressed as a matrix element of the position operator. 5,6,13,14 According to the modern theory of polarization, 5,6 the dipole moment per unit cell is a geometric quantum phase or a Berry phase. 15 A nontrivial quantum phase usually appears in cases when the system under study is coupled to the rest of the universe through some external parameter ζ. 16 The Berry phase γ can then be written as where ψ(ζ) is the wavefunction and C is a closed contour. In the case of periodic dipole moment, the parameter ζ is represented by a reciprocal vector k running over the Brillouin zone. Therefore, within single-determinant methods, Eq. 3 necessitates evaluating ∂/∂k derivatives of the Slater determinant. Since the latter is obtained in the self-consistent field (SCF) procedure which usually adds to the solution arbitrary complex phase prefactors depending on the k value, it is not trivial to differentiate the SCF solution with respect to k. Two main approaches to this problem have been proposed. The first method involves a discrete representation of ∂/∂k derivatives through matrix elements at neighboring k points. This representation is invariant with respect to phase arbitrariness which could arise from the SCF procedure. 5,6 The second approach introduces phase factors which cancel the SCF arbitrary phases and allow one to evaluate ∂/∂k by employing regular finite differences techniques. 17,18,19 Here, we start from the periodic dipole expression arising in the first approach. We consider the first approach as a simpler and computationally more robust alternative, since it does not require introduction of complex phase functions and band resolution techniques. Therefore, it can be seen as a black-box technique that has only the number of k points in the Brillouin zone sampling as an input. In addition, this method has been tested with plane waves before, 20,21 and its robustness was demonstrated not only for the linear response regime but for finite electric fields as well. Yet another approach to the periodic dipole derivatives is possible in the case of PWs, where these quantities require only dipole matrix elements between occupied and virtual Bloch functions φ nk |r|φ ak . These "transition" dipole elements can be transformed into the following well-defined expression 22,23 where H SCF is the crystal SCF Hamiltonian, ǫ a (k) and ǫ n (k) are orbital energies of empty and occupied bands. However, in the case of GTOs this approach could be less beneficial since Pulay's type of terms arise. To the best of our knowledge, the only implementation of IR intensities with GTOs was reported by Jacquemin and coworkers 24 using the second approach. Thus, the focus of the current work is in adapting the discretized Berry phase approach for evaluation of IR intensities with GTOs.
This work can be also seen as complementary to our previous work on zone centered phonon frequencies. 4 One should keep in mind that IR intensities are always obtained with vibrational frequencies, and therefore, any algorithm for IR intensities can be seen as an efficient one if its execution time is negligible with respect to that for vibrational frequencies.

II. THEORY
A. Berry phase approach to periodic dipole moment We first consider the formulation of periodic dipole moment through the Berry phase approach in localized basis sets proposed by Kudin and coworkers. 25 We restrict our consideration to closed-shell one-dimensional periodic systems aligned along the z axes. This makes our description simpler without losing essential details. Throughout this work we will be using Bloch orbitals can be expanded as According to the Berry phase approach, 5,6 the z component of the dipole moment per unit cell is Here, we avoid introducing an additional z subscript, since unless otherwise stated we always will consider the longitudinal part of the dipole moment. Substituting u nk (r) in the dipole expression by Eq. (7) we obtain the following representation for the periodic dipole moment and are dipole and overlap matrices in the AO basis. Calculation of the d 1 term is more efficient to perform in the AO basis where P 0g µν is the electron density matrix. The matrix (Z 0g µν − gaS 0g µν ) can be further split into Hermitian and anti-HermitianZ parts. Both matrices Z 0g µν andZ 0g µν are real, and we refer to Hermiticity of their reciprocal space counterparts Z µν (k) andZ µν (k). Since we contract both parts with the Hermitian matrix P 0g µν , the anti-Hermitian part can be safely omitted. This leads us to the following working expression The d 2 term requires special care because CO coefficients usually contain arbitrary kdependent phase factors from a diagonalization procedure. Therefore, one cannot apply regular numerical differentiation schemes directly to CO coefficients. To cope with this difficulty we use a discretization approach 25 which starts from rewriting Eq. (11) in the form where only the occupied orbital part of C(k) is engaged. For further discussion we introduce where O is the number of occupied bands. Treating k ′ as a variable and k as a parameter we can write Keeping in mind the orthogonality relation Σ(k, k) = 1 we rewrite Eq. 20 as The trace operation commutes with the differentiation Applying the well-known matrix relation tr ln(A) = ln det(A) to the Σ(k, k ′ ) matrix we Since ln det Σ(k, k) = 0, using simple rectangular discretization for the k ′ derivatives and the Brillouin zone integration we arrive to the following discretized form 25 for d 2 : By introducing arbitrary phase factors e iθ(k j ) in front of CO coefficients, one can easily see that they cancel each other after the summation over the Brillouin zone. Please note that one cannot evaluate Eq. (18) by using a perturbation (linear response) expression to build C(k + ∆k) from C(k). As was pointed out in Ref. 6, a perturbation theory implicitly uses the so called "parallel-transport" gauge which necessitates phase equality between C(k) and C(k + ∆k), and as a consequence produces zero Berry phase.
Similar to the molecular case, the transverse components of the periodic dipole moment can be written as where q = x or y and Q 0g µν = drµ 0 (r)qν g (r).
Therefore, their differentiation is very similar to that of the d 1 part of the longitudinal dipole moment, and it will not be considered further.

B. Dipole derivatives
Differentiation of Eq. (17) with respect to nuclear coordinate displacements gives Both terms in Eq. (28) are calculated with standard techniques: the coupled perturbed SCF (CPSCF) procedure for the density derivatives, 4,26 and a regular integral derivatives evaluation technique. 27 The CPSCF procedure evaluates the same density derivatives as in the case of vibrational frequency calculations; therefore computational overhead of the d 1 part is negligible, since it only involves calculation of one-electron terms.
In order to obtain the d 2 derivatives we adapt the algebraic approach discussed in Ref. 20.
The main step can be summarized as application of the following algebraic identity which is an analog of Eqs. 22 and 23. This was proven from a different point of view in Ref. 20. The differentiation of Σ(k j , k j+1 ) gives rise to two types of terms in the d 2 derivatives which correspond to Pulay's type of forces from overlap integral derivatives and response of CO coefficients To complete the construction of the dipole moment derivatives we will express CO derivatives via the response matrix U obtained in the CPSCF procedure 3 Thus we arrive to Since we do not evaluate R i dependent phases, the application of the linear response Eq. (34) is acceptable.
The computational complexity of the presented scheme is negligible with respect to that

A. Algorithmic tests
To validate and assess the present formalism, which has been implemented in the development version of the gaussian program, 28 we have tested it on a model one-dimensional chain of water molecules (see Fig. 1 and Table I) with the geometry optimized within the HF and PBE methods. 29 In order to describe hydrogen-bond interactions properly, we employ the 6-311++G** basis set in most of our tests. 30 In Table II we compare dipole derivatives obtained in periodic calculations with those from oligomeric estimations. Given that the 6-311++G** basis set is too large to perform vibrational frequency calculations for long oligomeric chains, we have used the 6-31G basis set for this comparison. The oligomeric results are calculated by the difference scheme which is less prone to edge effects where This suggests that the periodic dipole moment converges faster than its derivatives with the number of k points.
As a matter of practical interest we illustrate the dependence of the five normal mode frequencies and their IR intensities on the number of k points involved in the discretization of the Brillouin zone with HF and DFT (see Table III). As a representative DFT method we used the non-empirical generalized gradient approximation (GGA) functional of Perdew-Burke-Ernzerhof 31 (PBE). These model calculations suggest that 100 -1000 k points are enough to obtain IR intensities converged within 1%. A large difference in sensitivity of vibrational frequencies and IR intensities to N k should not be considered as a serious issue, since the overhead from enlarging N k is negligible in our vibrational frequency evaluation algorithm. 4 We suppose that the large difference in magnitudes of HF and PBE IR intensities  Fig. 1).
Method ∂d z /∂O z ∂d z /∂O y ∂d y /∂O z ∂d y /∂O y is related to the tendency of pure DFT functionals to yield a more metallic description than the HF method.

B. Poly(paraphenylenevinylene)
In order to illustrate the prediction capabilities of our approach, we chose the experimentally and theoretically well-studied one-dimensional system of poly(paraphenylenevinylene) (PPV, Fig. 2). 10,32,33,34 Such an interest is motivated by conductivity and luminescent nonlinear properties which make this system a promising candidate as a material for light-emitting diodes. 32 We would like to point out that previous theoretical work from our group on PPV was done with some limitations which are overcome in the current study. In Ref. 10, the IR intensities took into account only the changes in the non-periodic part of the dipole moment (d 1 ). This caused an underestimation of IR intensities for those vibrations that modify the longitudinal dipole moment of PPV.
As mentioned in the Introduction, IR intensities usually receive less attention and their values are presented only graphically in arbitrary units. 33,34 Therefore, in Table IV we compare IR-active calculated vibrational frequencies with their experimental counterparts. Note that in the case of IR intensities only a qualitative comparison with experimental data can be done (see Fig. 3). In our calculations we used the following functionals: the local spin density approximation 35 38 We also included the LC-ωPBE functional because we expect some overesti-mation of IR intensities by regular functionals due to the well-known problem with electric field response properties in extended systems. 39,40 According to Table IV   although the calculated vibrational frequency adequately reproduce those from experiment, the corresponding IR intensities do not always follow qualitatively correct trends. Application of the LC-ωPBE functional corrects IR intensities presumably due to the right 1/r asymptotic behavior of the LC-ωPBE exchange potential but worsens frequencies. We hope that our scheme and the results reported here will stimulate further functional development in the future. width.