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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

Ab initio analytical infrared intensities for periodic systems through a coupled perturbed Hartree-Fock/Kohn-Sham method.

A fully analytical method for calculating Born charges and, hence, infrared intensities of periodic systems, is formulated and implemented in the CRYSTAL program, which uses a local gaussian type basis set. Our efficient formalism combines integral gradients with first-order coupled perturbed Hartree-Fock/Kohn Sham electronic response to an electric field. It avoids numerical differentiation with respect to wave vectors, as in some Berry phase approaches, and with respect to atomic coordinates. No perturbation equations for the atomic displacements need to be solved. Several tests are carried out to verify numerical stability, consistency in one, two, and three dimensions, and applicability to large unit cells. Future extensions to piezoelectricity and Raman intensities are noted.

Static and dynamic coupled perturbed Hartree-Fock vibrational (hyper)polarizabilities of polyacetylene calculated by the finite field nuclear relaxation method.

The vibrational contribution to static and dynamic (hyper)polarizability tensors of polyacetylene are theoretically investigated. Calculations were carried out by the finite field nuclear relaxation (FF-NR) method for periodic systems, newly implemented in the CRYSTAL code, using the coupled perturbed Hartree-Fock scheme for the required electronic properties. The effect of the basis set is also explored, being particularly important for the non-periodic direction perpendicular to the polymer plane. Components requiring a finite (static) field in the longitudinal direction for evaluation by the FF-NR method were not evaluated. The extension to that case is currently being pursued. Whereas the effect on polarizabilities is relatively small, in most cases the vibrational hyperpolarizability tensor component is comparable to, or larger than the corresponding static electronic contribution.

The first and second static electronic hyperpolarizabilities of zigzag boron nitride nanotubes. An ab initio approach through the coupled perturbed Kohn-Sham scheme.

The coupled perturbed Kohn-Sham (CPKS) computational scheme for the evaluation of electric susceptibility tensors in periodic systems, recently implemented in the CRYSTAL code, has been extended to third-order. It is, then, used to obtain static electronic hyperpolarizabilities of zigzag BN nanotubes for the first time. This procedure, which is fully analytic in all key steps, requires a double self-consistent treatment for taking into account the first- and second-order response of the system to the applied field. The performance of different functionals is compared and the B3LYP hybrid is ultimately chosen for calculations on nanotubes having radii as large as R = 20 Å (6-200 atoms in the unit cell). Such large radii are sufficient to give the pure longitudinal component of the (hyper)polarizability tensors to within 1% of the "exact" hexagonal BN monolayer limit. Other tensor components involving the transverse direction converge more slowly. They can, however, be extrapolated to the monolayer limit to within 4% accuracy except for the pure transverse second hyperpolarizability, which has an error of 13% in that limit.

Electronic structure, bonding, spectra, and linear and nonlinear electric properties of Ti@C28.

The potential energy surface (PES) of Ti@C(28) has been revisited, and the stationary points have been carefully characterized. In particular, the C(2v) symmetry structure considered previously turns out to be a transition state lying 2.3 kcal/mol above the ground state of C(3v) symmetry at the MP2/6-31G(d) level. A large binding energy of 181.3 kcal/mol is found at the ROMP2/6-31G(d) level. Topological analysis of the generalized Ti@C(28) density reveals four bond paths between Ti and carbon atoms of the host. The character of all four contacts corresponds to a partially covalent closed shell interaction. UV-vis, IR, and Raman spectra are calculated and compared with C(28)H(4). The dipole moment and the static electronic and double harmonic vibrational (hyper)polarizabilities have been obtained. Distortion of the fullerene cage due to encapsulation leads to nonzero diagonal components of the electronic first hyperpolarizability ß, and to an increase in the diagonal components of the electronic polarizability α and second hyperpolarizability γ. However, introduction of the Ti atom causes a comparable or larger reduction in most cases due to localized bonding interactions. At the double harmonic level, the average vibrational ß is much larger than its electronic counterpart, but the opposite is true for α and for the contribution to γ that has been calculated. There is also a very large anharmonic (nuclear relaxation) contribution to ß which results from a shallow PES with four minima separated by very low barriers. Thus, the vibrational γ (and α) may, likewise, become much larger when anharmonicity is taken into account.

Communication: convergence of anharmonic infrared intensities of hydrogen fluoride in traditional and explicitly correlated coupled cluster calculations.

It is shown that the convergence of anharmonic infrared spectral intensities with respect to the basis set size is much enhanced in explicitly correlated calculations as compared to traditional configuration interaction type wave function expansion. Explicitly correlated coupled cluster (CC) calculations using Slater-type geminal correlation factor (CC-F12) yield well-converged dipole derivatives and vibrational intensities for hydrogen fluoride with basis set involving f functions on the heavy atom. Combination of CC-F12 with singles, doubles, and non-iterative triples (CCSD(T)-F12) with small corrections due to quadruple excitations, core-electron correlation, and relativistic effects yields vibrational line positions, dipole moments, and transition dipole matrix elements in good agreement with the best experimental values.

Electric field polarization in conventional density functional theory: from quasilinear to two-dimensional and three-dimensional extended systems.

The large overshoot in (hyper)polarizabilities of quasilinear (1D) chains calculated by applying density functional theory with conventional functionals is investigated for several 2D and 3D extended systems. These systems include arrays of molecular hydrogen chains, as well as 2D coronene-type structures and LiF in 1D, 2D, and 3D. Contrary to a recently proposed model it is found that the overshoot persists in all of these cases. A simple explanation is provided by an analysis of the field-induced charges for molecular hydrogen, which shows an excessive buildup at the chain ends regardless of where the chain is located within the 2D and 3D array.