Anharmonic Vibrational States of Solids from DFT Calculations. Part I: Description of the Potential Energy Surface.

A computational approach is presented to compute anharmonic vibrational states of solids from quantum-mechanical DFT calculations by taking into explicit account phonon-phonon couplings via the vibrational configuration interaction (VCI) method. The Born-Oppenheimer potential energy surface (PES) is expanded in a Taylor's series in terms of harmonic normal coordinates, centered at the equilibrium nuclear configuration, is truncated to quartic order, and contains one-mode, two-mode, and three-mode interatomic force constants. The description of the anharmonic terms of the PES involves the numerical evaluation of high-order energy derivatives (cubic and quartic in our case) with respect to nuclear displacements and constitutes the most computationally demanding step in the characterization of anharmonic vibrational states of materials. Part I is devoted to the description of the PES. Four different numerical approaches are presented for the description of the potential, all based on a grid representation of the PES in the basis of the normal coordinates, that require different ingredients (energy and/or forces) to be evaluated at each point (i.e., nuclear configuration) of the grid. The numerical stability and relative computational efficiency of the various schemes for the description of the PES are discussed on two molecular systems (water and methane) and two extended solids (Ice-XI and MgH2). All the presented algorithms are implemented into a developmental version of the Crystal program.

Anharmonic Vibrational States of Solids from DFT Calculations. Part II: Implementation of the VSCF and VCI Methods.

Two methods are implemented in the Crystal program for the calculation of anharmonic vibrational states of solids: the vibrational self-consistent field (VSCF) and the vibrational configuration-interaction (VCI). While the former is a mean-field approach, where each vibrational mode interacts with the average potential of the others, the latter allows for an explicit and complete account of mode-mode correlation. Both schemes are based on the representation of the adiabatic potential energy surface (PES) discussed in Part I, where the PES is expanded in a Taylor's series so as to include up to cubic and quartic terms. The VSCF and VCI methods are formally presented and their numerical parameters discussed. In particular, the convergence of computed anharmonic vibrational states, within the VCI method, is investigated as a function of the truncation of the expansion of the nuclear wave function. The correctness and effectiveness of the implementation is discussed by comparing with available theoretical and experimental data on both molecular and periodic systems. The effect of the adopted basis set and exchange-correlation functional in the description of the PES on computed anharmonic vibrational states is also addressed.

Symmetry and random sampling of symmetry independent configurations for the simulation of disordered solids.

A symmetry-adapted algorithm producing uniformly at random the set of symmetry independent configurations (SICs) in disordered crystalline systems or solid solutions is presented here. Starting from Pólya's formula, the role of the conjugacy classes of the symmetry group in uniform random sampling is shown. SICs can be obtained for all the possible compositions or for a chosen one, and symmetry constraints can be applied. The approach yields the multiplicity of the SICs and allows us to operate configurational statistics in the reduced space of the SICs. The present low-memory demanding implementation is briefly sketched. The probability of finding a given SIC or a subset of SICs is discussed as a function of the number of draws and their precise estimate is given. The method is illustrated by application to a binary series of carbonates and to the binary spinel solid solution Mg(Al,Fe)2O4.

The electronic structure of MgO nanotubes. An ab initio quantum mechanical investigation.

The structural, vibrational and response properties of the (n,0) and (m,m) MgO nanotubes are computed by using a Gaussian type basis set, a hybrid functional (B3LYP) and the CRYSTAL09 code. Tubes in the range 6 ≤ n ≤ 140 and 3 ≤ m ≤ 70 were considered, being n = 2 × m the number of MgO units in the unit cell (so, the maximum number of atoms is 280). Tubes are built by rolling up the fully relaxed 2-D conventional cell (2 MgO units, with oxygen atoms protruding from the Mg plane alternately up and down by 0.38 Å). The relative stability of the (n,0) with respect to the (m,m) family, the relaxation energy and equilibrium geometry, the band gap, the IR vibrational frequencies and intensities, and the electronic and ionic contributions to the polarizability are reported. All these properties are shown to converge smoothly to the monolayer values. Absence of negative vibrational frequencies confirms that the tubes have a stable structure. The parallel component of the polarizability α(∥) converges very rapidly to the monolayer value, whereas α(⊥) is still changing at n = 140; however, when extrapolated to very large n values, it coincides with the monolayer value to within 1%. The electronic contribution to α is in all cases (α(∥) and α(⊥); 6 ≤ n ≤ 140) smaller than the vibrational contribution by about a factor of three, at variance with respect to more covalent tubes such as the BN ones, for which the ratio between the two contributions is reversed.

On the use of symmetry in configurational analysis for the simulation of disordered solids.

The starting point for a quantum mechanical investigation of disordered systems usually implies calculations on a limited subset of configurations, generated by defining either the composition of interest or a set of compositions ranging from one end member to another, within an appropriate supercell of the primitive cell of the pure compound. The way in which symmetry can be used in the identification of symmetry independent configurations (SICs) is discussed here. First, Pólya's enumeration theory is adopted to determine the number of SICs, in the case of both varying and fixed composition, for colors numbering two or higher. Then, De Bruijn's generalization is presented, which allows analysis of the case where the colors are symmetry related, e.g. spin up and down in magnetic systems. In spite of their efficiency in counting SICs, neither Pólya's nor De Bruijn's theory helps in solving the difficult problem of identifying the complete list of SICs. Representative SICs are obtained by adopting an orderly generation approach, based on lexicographic ordering, which offers the advantage of avoiding the (computationally expensive) analysis and storage of all the possible configurations. When the number of colors increases, this strategy can be combined with the surjective resolution principle, which permits the efficient generation of SICs of a problem in |R| colors starting from the ones obtained for the (|R| - 1)-colors case. The whole scheme is documented by means of three examples: the abstract case of a square with C(4v) symmetry and the real cases of the garnet and olivine mineral families.

Accurate dynamical structure factors from ab initio lattice dynamics: the case of crystalline silicon.

A fully ab initio technique is discussed for the determination of dynamical X-ray structure factors (XSFs) of crystalline materials, which is based on a standard Debye-Waller (DW) harmonic lattice dynamical approach with all-electron atom-centered basis sets, periodic boundary conditions, and one-electron Hamiltonians. This technique requires an accurate description of the lattice dynamics and the electron charge distribution of the system. The main theoretical parameters involved and final accuracy of the technique are discussed with respect to the experimental determinations of the XSFs at 298 K of crystalline silicon. An overall agreement factor of 0.47% between the ab initio predicted values and the experimental determinations is found. The best theoretical determination of the anisotropic displacement parameter, of silicon is here 60.55 × 10(-4) Å(2), corresponding to a DW factor B = 0.4781 Å(2).

Nuclear motion effects on the density matrix of crystals: an ab initio Monte Carlo harmonic approach.

In the frame of the Born-Oppenheimer approximation, nuclear motions in crystals can be simulated rather accurately using a harmonic model. In turn, the electronic first-order density matrix (DM) can be expressed as the statistically weighted average over all its determinations each resulting from an instantaneous nuclear configuration. This model has been implemented in a computational scheme which adopts an ab initio one-electron (Hartree-Fock or Kohn-Sham) Hamiltonian in the CRYSTAL program. After selecting a supercell of reasonable size and solving the corresponding vibrational problem in the harmonic approximation, a Metropolis algorithm is adopted for generating a sample of nuclear configurations which reflects their probability distribution at a given temperature. For each configuration in the sample the "instantaneous" DM is calculated, and its contribution to the observables of interest is extracted. Translational and point symmetry of the crystal as reflected in its average DM are fully exploited. The influence of zero-point and thermal motion of nuclei on such important first-order observables as x-ray structure factors and Compton profiles can thus be estimated.

A new massively parallel version of CRYSTAL for large systems on high performance computing architectures.

Fully ab initio treatment of complex solid systems needs computational software which is able to efficiently take advantage of the growing power of high performance computing (HPC) architectures. Recent improvements in CRYSTAL, a periodic ab initio code that uses a Gaussian basis set, allows treatment of very large unit cells for crystalline systems on HPC architectures with high parallel efficiency in terms of running time and memory requirements. The latter is a crucial point, due to the trend toward architectures relying on a very high number of cores with associated relatively low memory availability. An exhaustive performance analysis shows that density functional calculations, based on a hybrid functional, of low-symmetry systems containing up to 100,000 atomic orbitals and 8000 atoms are feasible on the most advanced HPC architectures available to European researchers today, using thousands of processors.

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.

Ab initio modeling of TiO2 nanotubes.

TiO(2) nanotubes constructed from a lepidocrocite-like TiO(2) layer were investigated with ab initio methods employing the periodic CRYSTAL code. The dependence of strain energies, structural and electronic properties on the tube diameter was investigated in the 18-57 A range. Nanotubes constructed by a (0,n) rollup proved to be the most stable at all diameters. All three types of rollup undergo significant reconstruction at diameters <25 A. All investigated structures possess a high ( approximately 5.4 eV) band gap compared to bulk TiO(2) phases (3.96 and 4.63 eV for rutile and anatase calculated with the same functional and basis set).