The shape effect and its consequences for polar surfaces and for heterogeneous catalysis.

In this paper we develop the shape effect, which is relevant for crystalline materials whose size is larger than that of the thermodynamic limit. According to this effect the electronic properties of one surface of a crystal depend upon all of its surfaces, i.e. on the overall shape. At first, qualitative mathematical arguments are presented for the existence of this effect based on the conditions for the stability of polar surfaces. Our treatment explains why such surfaces are observed even though earlier theory indicated that they should not exist. Then, models are developed from which it is found computationally that changing the shape of a polar crystal can substantially alter the magnitude of its surface charges. Apart from surface charges, it follows that the crystal shape will also significantly affect bulk properties, most notably polarization and piezoelectric responses. Additional model calculations show a strong shape effect on the activation energy for heterogeneous catalysis primarily through local surface charges rather than a non-local/long range electrostatic potential.

Perturbation Theory Treatment of Spin-Orbit Coupling. III: Coupled Perturbed Method for Solids.

A previously proposed noncanonical coupled-perturbed Kohn-Sham density functional theory (KS-DFT)/Hartree-Fock (HF) treatment for spin-orbit coupling is here generalized to infinite periodic systems. The scalar-relativistic periodic KS-DFT/HF solution, obtained with a relativistic effective core potential, is taken as the zeroth-order approximation. Explicit expressions are given for the total energy through third-order, which satisfy the 2N + 1 rule (i.e., requiring only the first-order perturbed wave function for determining the energy through third-order). Expressions for additional second-order corrections to the perturbed wave function (as well as related one-electron properties) are worked out at the uncoupled-perturbed level of theory. The approach is implemented in the Crystal program and validated with calculations of the total energy, electronic band structure, and density variables of spin-current DFT on the tungsten dichalcogenide hexagonal bilayer series (i.e., WSe2, WTe2, WPo2, WLv2), including 6p and 7p elements as a stress test. The computed properties through second- or third-order match well with those from reference two-component self-consistent field (2c-SCF) calculations. For total energies, E(3) was found to consistently improve the agreement against the 2c-SCF reference values. For electronic band structures, visible differences w.r.t. 2c-SCF remained through second-order in only the single-most difficult case of WLv2. As for density variables of spin-current DFT, the perturbed electron density, being vanishing in first-order, is the most challenging for the perturbation theory approach. The visible differences in the electron densities are, however, largest close to the core region of atoms and smaller in the valence region. Perturbed spin-current densities, on the other hand, are well reproduced in all tested cases.

CRYSTAL23: A Program for Computational Solid State Physics and Chemistry.

The Crystal program for quantum-mechanical simulations of materials has been bridging the realm of molecular quantum chemistry to the realm of solid state physics for many years, since its first public version released back in 1988. This peculiarity stems from the use of atom-centered basis functions within a linear combination of atomic orbitals (LCAO) approach and from the corresponding efficiency in the evaluation of the exact Fock exchange series. In particular, this has led to the implementation of a rich variety of hybrid density functional approximations since 1998. Nowadays, it is acknowledged by a broad community of solid state chemists and physicists that the inclusion of a fraction of Fock exchange in the exchange-correlation potential of the density functional theory is key to a better description of many properties of materials (electronic, magnetic, mechanical, spintronic, lattice-dynamical, etc.). Here, the main developments made to the program in the last five years (i.e., since the previous release, Crystal17) are presented and some of their most noteworthy applications reviewed.

Perturbation Theory Treatment of Spin-Orbit Coupling II: A Coupled Perturbed Kohn-Sham Method.

A noncanonical coupled perturbed Kohn-Sham density functional theory (KS-DFT)/Hartree-Fock (HF) treatment of spin-orbit coupling (SOC) is provided. We take the scalar-relativistic KS-DFT/HF solution, obtained with a relativistic effective core potential, as the zeroth-order approximation. Explicit expressions are given for the total energy through the 4th order, which satisfy the 2n + 1 rule. Second-order expressions are provided for orbital energies and density variables of spin-current DFT. Test calculations are carried out on the halogen homonuclear diatomic and hydride molecules, including 6p and 7p elements, as well as open-shell negative ions. The computed properties through second or third order match well with those from reference two-component self-consistent field calculations for total and orbital energies as well as spin-current densities. In only one case (At2-) did a significant deviation occur for the remaining density variables. Our coupled perturbation theory approach provides an efficient way of adding the effect of SOC to a scalar-relativistic single-reference KS-DFT/HF treatment, in particular because it does not require diagonalization in the two-component spinor basis, leading to saving factors on the number of required floating-point operations that may exceed one order of magnitude.

Perturbation Theory Treatment of Spin-Orbit Coupling, Part I: Double Perturbation Theory Based on a Single-Reference Initial Approximation.

We develop a perturbation theory for solving the many-body Dirac equation within a given relativistic effective-core potential approximation. Starting from a scalar-relativistic unrestricted Hartree-Fock (SR UHF) solution, we carry out a double perturbation expansion in terms of spin-orbit coupling (SOC) and the electron fluctuation potential. Computationally convenient energy expressions are derived through fourth order in SOC, second order in the electron fluctuation potential, and a total of third order in the coupling between the two. Illustrative calculations on the halogen series of neutral and singly positive diatomic molecules show that the perturbation expansion is well-converged by taking into account only the leading (nonvanishing) term at each order of the electron fluctuation potential. Our perturbation theory approach provides a computationally attractive alternative to a two-component self-consistent field treatment of SOC. In addition, it includes coupling with the fluctuation potential through third order and can be extended (in principle) to multireference calculations, when necessary for both closed- and open-shell cases, using quasi-degenerate perturbation theory.

First-Principles Calculation of the Optical Rotatory Power of Periodic Systems: Application on α-Quartz, Tartaric Acid Crystal, and Chiral (n,m)-Carbon Nanotubes.

The self-consistent coupled-perturbed (SC-CP) method in the CRYSTAL program has been adapted to obtain electromagnetic optical rotation properties of chiral periodic systems based on the calculation of the magnetic moment induced by the electric field. Toward that end, an expression for the magnetic transition moment is developed, which involves an appropriate electronic angular momentum operator. This operator is forced to be hermitian so that the chiroptical properties are real. In our formulation, the trace of the optical rotatory power matrix is gauge-origin-invariant as long as the electric dipole transition matrix elements are obtained using the velocity (rather than position) operator. On the other hand, the component along the optic axis is invariant in general for uniaxial and biaxial crystals. Under the same conditions, these properties also do not depend on the so-called missing integers that occur in the treatment of the electric dipole moment of quasi-one-dimensional periodic systems or the analogue of missing integers for the case of higher dimensionality. Tests on a model H2O2 polymer confirm the formalism and, as desired, show that the calculated properties are independent of the size and definition of the unit cell. In addition, an empirical relation to a finite oligomer gauge-including atomic orbital (GIAO) calculation is found. Applications, with comparison to experiment, are carried for α-quartz, tartaric acid crystal, and carbon nanotubes. Future developments of this initial approach to chiroptical properties in the solid state are noted.

The CRYSTAL code, 1976-2020 and beyond, a long story.

CRYSTAL is a periodic ab initio code that uses a Gaussian-type basis set to express crystalline orbitals (i.e., Bloch functions). The use of atom-centered basis functions allows treating 3D (crystals), 2D (slabs), 1D (polymers), and 0D (molecules) systems on the same grounds. In turn, all-electron calculations are inherently permitted along with pseudopotential strategies. A variety of density functionals are implemented, including global and range-separated hybrids of various natures and, as an extreme case, Hartree-Fock (HF). The cost for HF or hybrids is only about 3-5 times higher than when using the local density approximation or the generalized gradient approximation. Symmetry is fully exploited at all steps of the calculation. Many tools are available to modify the structure as given in input and simplify the construction of complicated objects, such as slabs, nanotubes, molecules, and clusters. Many tensorial properties can be evaluated by using a single input keyword: elastic, piezoelectric, photoelastic, dielectric, first and second hyperpolarizabilities, etc. The calculation of infrared and Raman spectra is available, and the intensities are computed analytically. Automated tools are available for the generation of the relevant configurations of solid solutions and/or disordered systems. Three versions of the code exist: serial, parallel, and massive-parallel. In the second one, the most relevant matrices are duplicated on each core, whereas in the third one, the Fock matrix is distributed for diagonalization. All the relevant vectors are dynamically allocated and deallocated after use, making the code very agile. CRYSTAL can be used efficiently on high performance computing machines up to thousands of cores.

Coupled Perturbation Theory Approach to Dual Basis Sets for Molecules and Solids. 1. General Theory and Application to Molecules.

We present a new coupled Hartree-Fock(HF)/Kohn-Sham DFT perturbation method that accounts for the effect of enlarging the basis set in electronic structure calculations. In contrast with previous approaches, our dual basis set treatment yields not only a correction for the total energy but also for the orbital eigenvalues and density. The zeroth order solution is obtained from the projection of the small basis set coefficients. Diagonalization of the full Fock matrix in the large basis set is avoided. In this first paper of a series, we develop the theoretical foundations of our approach for molecules, including the coupled-perturbed equations through second order and the energy expressions through fourth order-as our method complies with Wigner's 2n + 1 rule. The first-order perturbation equation turns out to be uncoupled, and odd-order terms in the energy expansion vanish. In calculations on simple molecules, our method recovers over 93% (84%) of the missing DFT(HF) energy when going from the cc-pVDZ to the aug-cc-pVDZ basis, and over about 95% in all cases if an energy extrapolation formula is used. Mulliken charges, the orbital eigenvalue spectrum, and HOMO-LUMO gaps of the large basis are well reproduced. Charge density maps show that the differences between the perturbatively corrected density and the reference nearly vanish through second-order.

Biradical rotamer states tune electron J coupling and MAS dynamic nuclear polarization enhancement.

Cross Effect (CE) Dynamic Nuclear Polarization (DNP) relies on the dipolar (D) and exchange (J) coupling interaction between two electron spins. Until recently only the electron spin D coupling was explicitly included in quantifying the DNP mechanism. Recent literature discusses the potential role of J coupling in DNP, but does not provide an account of the distribution and source of electron spin J coupling of commonly used biradicals in DNP. In this study, we quantified the distribution of electron spin J coupling in AMUPol and TOTAPol biradicals using a combination of continuous wave (CW) X-band electron paramagnetic resonance (EPR) lineshape analysis in a series of solvents and at variable temperatures in solution - a state to be vitrified for DNP. We found that both radicals show a temperature dependent distribution of J couplings, and the source of this distribution to be conformational dynamics. To qualify this conformational dependence of J coupling in both molecules we carry out Broken Symmetry DFT calculations which show that the biradical rotamer distribution can account for a large distribution of J couplings, with the magnitude of J coupling directly depending on the relative orientation of the electron spin pair. We demonstrate that the electron spin J couplings in both AMUPol and TOTAPol span a much wider distribution than suggested in the literature. We affirm the importance of electron spin J coupling for DNP with density matrix simulations of DNP in Liouville space and under magic angle spinning, showcasing that a rotamer with high J coupling and optimum relative g-tensor orientation can significantly boost the DNP performance compared to random orientations of the electron spin pair. We conclude that moderate electron spin J coupling above a threshold value can facilitate DNP enhancements.

Electronic orbital response of regular extended and infinite periodic systems to magnetic fields. I. Theoretical foundations for static case.

A theoretical treatment for the orbital response of an infinite, periodic system to a static, homogeneous, magnetic field is presented. It is assumed that the system of interest has an energy gap separating occupied and unoccupied orbitals and a zero Chern number. In contrast to earlier studies, we do not utilize a perturbation expansion, although we do assume the field is sufficiently weak that the occurrence of Landau levels can be ignored. The theory is developed by analyzing results for large, finite systems and also by comparing with the analogous treatment of an electrostatic field. The resulting many-electron Hamilton operator is forced to be hermitian, but hermiticity is not preserved, in general, for the subsequently derived single-particle operators that determine the electronic orbitals. However, we demonstrate that when focusing on the canonical solutions to the single-particle equations, hermiticity is preserved. The issue of gauge-origin dependence of approximate solutions is addressed. Our approach is compared with several previously proposed treatments, whereby limitations in some of the latter are identified.

Surface effects on converse piezoelectricity of crystals.

The contribution of surface units to bulk properties are often neglected in theoretical and computational studies of crystalline systems. We demonstrate that this assumption has to be made with caution in the case of (electric field) polarization. As a generalization of an earlier work on quasi-one-dimensional systems [Springborg, et al., Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 82, 165442], it is shown that the polarization for 2D and 3D systems contains a surface contribution that can, in principle, take any value (within physical limits) and has consequences for converse piezoelectric responses. Subsequently, we determine the surface effects quantitatively for a group of ferroelectric perovskite structures. Our results indicate that such contributions can be substantial.

Numerical and exact kinetic energy operator using Eckart conditions with one or several reference geometries: Application to HONO.

For the computation of rovibrational levels and their spectroscopic intensities, the Eckart conditions are essential to achieve the optimal separation between rotation and vibration. Dymarsky and Kudin [J. Chem. Phys. 122, 124103 (2005)] proposed a procedure for a simplified calculation of the Eckart rotation matrix. In the present work, we have adapted their approach to obtain a kinetic energy operator in curvilinear coordinates using a numerical but exact procedure without resorting to finite differences. Furthermore, we have modified this approach for the study of molecular systems with several minima, for which several Eckart reference geometries are required. The HONO molecular system has been used to show the efficiency of our implementation. Using the Eckart conditions with multi-reference geometries allows for a calculation of the rotational levels as well as frequencies and intensities of the infrared spectra of both HONO isomers with a single calculation.

Computation of Second Harmonic Generation for Crystalline Urea and KDP. An ab Initio Approach through the Coupled Perturbed Hartree-Fock/Kohn-Sham Scheme.

The electronic second harmonic generation (SHG) tensor, d, of crystalline urea and potassium dihydrogen phosphate (KDP) is evaluated as a function of frequency using a Gaussian type basis set and the Coupled Perturbed Hartree-Fock (CPHF) and Kohn-Sham (CPKS) schemes as implemented in the CRYSTAL code. The results of various functionals, including LDA, GGA (PBE), and global and range-separated hybrids (B3LYP, PBE0, LC-BLYP), as well as Hartree-Fock, are compared. It is found that the calculated SHG intensity always decreases as the percentage of exact exchange increases. The hybrid functionals turn out to provide results that agree well with experiment. For urea and KDP the percentage of exact exchange determined by the inverse dielectric constant is too large. At 1064 nm the vibrational contribution for urea is found to be less than 5% of the total value. To the authors' knowledge, this is the first coupled (self-consistent) calculation of SHG for any periodic system.

Calculation of the dynamic first electronic hyperpolarizability ß(-ω(σ); ω(1), ω(2)) of periodic systems. Theory, validation, and application to multi-layer MoS2.

We describe our implementation of a fully analytical scheme, based on the 2n + 1 rule, for computing the coupled perturbed Hartree Fock and Kohn-Sham dynamic first hyperpolarizability tensor ß(-ωσ; ω1, ω2) of periodic 1D (polymer), 2D (slab), and 3D (crystal) systems in the CRYSTAL code [R. Dovesi et al., Int. J. Quantum Chem. 114, 1287 (2014)], which utilizes local Gaussian type basis sets. The dc-Pockels (dc-P) and second harmonic generation (SHG) tensors are included as special cases. It is verified that (i) symmetry requirements are satisfied; (ii) using LiF as an example, the infinite periodic polymer result agrees with extrapolated finite oligomer calculations and, likewise, for the build-up to a 2D slab and a 3D crystal; (iii) the values converge to the static case for low frequencies; and (iv) the Bishop-deKee dispersion formulas relating dc-P, SHG, and general processes are reproduced through quartic terms. Preliminary SHG calculations on multi-layer MoS2 satisfactorily reproduce experimental data.

Electronic responses of long chains to electrostatic fields: Hartree-Fock vs. density-functional theory: a model study.

The response to an electrostatic field is determined through simple model calculations, within both the restricted Hartree-Fock and density functional theory methods, for long, finite as well as infinite, periodic chains. The permanent dipole moment, µ0, the polarizability, α, and the hyperpolarizabilities ß and γ, calculated using a finite-field approach, are extensively analyzed. Our simple model allows for treatment of large systems and for separation of the properties into atomic and unit-cell contributions. That part of the response properties attributable to the terminations of the finite system change into delocalized current contributions in the corresponding infinite periodic system. Special emphasis is placed on analyzing the reasons behind the dramatic overestimation of the response properties found with density functional theory methods presently in common use.

Ab initio analytical Raman intensities for periodic systems through a coupled perturbed Hartree-Fock/Kohn-Sham method in an atomic orbital basis. I. Theory.

We present a fully analytical formulation for calculating Raman intensities of crystalline periodic systems using a local basis set. Numerical differentiation with respect to atomic coordinates and with respect to wavevectors is entirely avoided as is the determination of crystal orbital coefficient derivatives with respect to nuclear displacements. Instead, our method utilizes the orbital energy-weighted density matrix and is based on the self-consistent solution of first- and second-order Coupled Perturbed Hartree-Fock/Kohn-Sham equations for the electronic response to external electric fields at the equilibrium geometry. This method has also been implemented in the Crystal program, which uses a Gaussian type basis set.

Ab initio analytical Raman intensities for periodic systems through a coupled perturbed Hartree-Fock/Kohn-Sham method in an atomic orbital basis. II. Validation and comparison with experiments.

In this work, we validate a new, fully analytical method for calculating Raman intensities of periodic systems, developed and presented in Paper I [L. Maschio, B. Kirtman, M. Rérat, R. Orlando, and R. Dovesi, J. Chem. Phys. 139, 164101 (2013)]. Our validation of this method and its implementation in the CRYSTAL code is done through several internal checks as well as comparison with experiment. The internal checks include consistency of results when increasing the number of periodic directions (from 0D to 1D, 2D, 3D), comparison with numerical differentiation, and a test of the sum rule for derivatives of the polarizability tensor. The choice of basis set as well as the Hamiltonian is also studied. Simulated Raman spectra of α-quartz and of the UiO-66 Metal-Organic Framework are compared with the experimental data.

Raman spectrum of pyrope garnet. A quantum mechanical simulation of frequencies, intensities, and isotope shifts.

The Raman spectrum of pyrope garnet is simulated in ab initio quantum mechanical calculations, using an all-electron Gaussian-type basis set and the hybrid B3LYP functional. Frequencies calculated for the 25 Raman-active modes are in excellent agreement with the several sets of experimental data, with the mean absolute difference ranging from 4 to 8 cm(-1). Comparison of the computed and experimental spectrum shows excellent agreement for most of the intensities as well. Modes missing from experiment are shown to be characterized by low (computed) intensity. Spurious peaks in the experimental spectra are also identified. The isotopic effect has been simulated for (24)Mg → (26)Mg substitution and shows excellent agreement with shifts reported in one of the experiments. Agreement is excellent for all but one mode, which turns out to be attributed to the wrong symmetry in the experiment.

Performance of density functional theory in computing nonresonant vibrational (hyper)polarizabilities.

A set of exchange-correlation functionals, including BLYP, PBE0, B3LYP, BHandHLYP, CAM-B3LYP, LC-BLYP, and HSE, has been used to determine static and dynamic nonresonant (nuclear relaxation) vibrational (hyper)polarizabilities for a series of all-trans polymethineimine (PMI) oligomers containing up to eight monomer units. These functionals are assessed against reference values obtained using the Møller-Plesset second-order perturbation theory (MP2) and CCSD methods. For the smallest oligomer, CCSD(T) calculations confirm the choice of MP2 and CCSD as appropriate for assessing the density functionals. By and large, CAM-B3LYP is the most successful, because it is best for the nuclear relaxation contribution to the static linear polarizability, intensity-dependent refractive index second hyperpolarizability, static second hyperpolarizability, and is close to the best for the electro-optical Pockels effect first hyperpolarizability. However, none of the functionals perform satisfactorily for all the vibrational (hyper)polarizabilities studied. In fact, in the case of electric field-induced second harmonic generation all of them, as well as the Hartree-Fock approximation, yield the wrong sign. We have also found that the Pople 6-31+G(d) basis set is unreliable for computing nuclear relaxation (hyper)polarizabilities of PMI oligomers due to the spurious prediction of a nonplanar equilibrium geometry.