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.

Calculation of longitudinal polarizability and second hyperpolarizability of polyacetylene with the coupled perturbed Hartree-Fock/Kohn-Sham scheme: where it is shown how finite oligomer chains tend to the infinite periodic polymer.

The longitudinal polarizability, α(xx), and second hyperpolarizability, γ(xxxx), of polyacetylene are evaluated by using the coupled perturbed Hartree-Fock/Kohn-Sham (HF/KS) scheme as implemented in the periodic CRYSTAL code and a split valence type basis set. Four different density functionals, namely local density approximation (LDA) (pure local), Perdew-Becke-Ernzerhof (PBE) (gradient corrected), PBE0, and B3LYP (hybrid), and the Hartree-Fock Hamiltonian are compared. It is shown that very tight computational conditions must be used to obtain well converged results, especially for γ(xxxx), that is, very sensitive to the number of k(->) points in reciprocal space when the band gap is small (as for LDA and PBE), and to the extension of summations of the exact exchange series (HF and hybrids). The band gap in LDA is only 0.01 eV: at least 300 k(->) points are required to obtain well converged total energy and equilibrium geometry, and 1200 for well converged optical properties. Also, the exchange series convergence is related to the band gap. The PBE0 band gap is as small as 1.4 eV and the exchange summation must extend to about 130 Å from the origin cell. Total energy, band gap, equilibrium geometry, polarizability, and second hyperpolarizability of oligomers -(C(2)H(2))(m)-, with m up to 50 (202 atoms), and of the polymer have been compared. It turns out that oligomers of that length provide an extremely poor representation of the infinite chain polarizability and hyperpolarizability when the gap is smaller than 0.2 eV (that is, for LDA and PBE). Huge differences are observed on α(xx) and γ(xxxx) of the polymer when different functionals are used, that is in connection to the well-known density functional theory (DFT) overshoot, reported in the literature about short oligomers: for the infinite model the ratio between LDA (or PBE) and HF becomes even more dramatic (about 500 for α(xx) and 10(10) for γ(xxxx)). On the basis of previous systematic comparisons of results obtained with various approaches including DFT, HF, Moller-Plesset (MP2) and coupled cluster for finite chains, we can argue that, for the infinite chain, the present HF results are the most reliable.

The calculation of the static first and second susceptibilities of crystalline urea: A comparison of Hartree-Fock and density functional theory results obtained with the periodic coupled perturbed Hartree-Fock/Kohn-Sham scheme.

The static polarizability alpha and first hyperpolarizability beta tensors of crystalline urea and the corresponding first-(chi((1))) and second-(chi((2))) susceptibilities are calculated and compared to the same quantities obtained for the molecule by using the same code (a development version of CRYSTAL), basis set, and level of theory. In order to separate geometrical and solid state effects, two geometries are considered for the molecule in its planar conformation: (i) as cut out from the bulk structure and (ii) fully optimized. First, the effect of basis sets on computed properties is explored at the B3LYP level by employing basis sets of increasing complexity, from 6-31G(d,p) to 6-311G(2df,2pd) (Pople's family) and from DZP to QZVPPP (Thakkar/Ahlrichs/Dunning's family) on alpha and beta for both the molecule and the bulk. Then, five different levels of theory, namely, SVWN (local density approximation), PBE (generalized gradient approximation), PBE0 and B3LYP (hybrid), and Hartree-Fock are compared in combination with a TZPP basis set. Present results show that hybrid methods, in particular, B3LYP, are remarkably successful in predicting correctly both the first and second susceptibilities of urea bulk when combined at least with a triple-zeta quality basis set containing a double set of polarization functions. It is also shown that diffuse functions that are needed for molecular calculations are less crucial for the crystalline structure, as expected. Indeed, B3LYP/TZPP computed chi((1)) and chi((2)) tensor components (chi(aa) ((1))=1.107, chi(cc) ((1))=1.459, and chi((2))=-0.93 a.u.) are in very good agreement with experimental values. At variance with respect to previous periodic ab initio calculations, but in agreement with recent supermolecular results, the negative sign of chi((2)) is confirmed. Overall, static linear and nonlinear optical properties such as dielectric constants, refractive, and birefringence indices and second-harmonic generation coefficient of crystalline urea are very well reproduced by present calculations.

Calculation of the dielectric constant epsilon and first nonlinear susceptibility chi((2)) of crystalline potassium dihydrogen phosphate by the coupled perturbed Hartree-Fock and coupled perturbed Kohn-Sham schemes as implemented in the CRYSTAL code.

The high-frequency dielectric varepsilon and the first nonlinear electric susceptibility chi((2)) tensors of crystalline potassium dihydrogen phosphate (KH(2)PO(4)) are calculated by using the coupled perturbed Hartree-Fock and Kohn-Sham methods as implemented in the CRYSTAL code. The effect of basis sets of increasing size on varepsilon and chi((2)) is explored. Five different levels of theory, namely, local-density approximation, generalized gradient approximation (PBE), hybrids (B3LYP and PBE0), and HF are compared using the experimental and theoretical structures corresponding not only to the tetragonal geometry I4d2 at room temperature but also to the orthorhombic phase Fdd2 at low temperature. Comparison between the two phases and their optical behavior is made. The calculated results for the tetragonal phase are in good agreement with the experimental data.

Calculation of the static electronic second hyperpolarizability or chi(3) tensor of three-dimensional periodic compounds with a local basis set.

The coupled perturbed Hartree-Fock (CPHF) method for evaluating static first (beta) and second (gamma) hyperpolarizability tensors of periodic systems has recently been implemented in the CRYSTAL code [Bishop et al., J. Chem. Phys. 114, 7633 (2001)]. We develop here an efficient and accurate computational protocol, along with the local basis sets needed for first and second row atoms. Application is made to several high symmetry three-dimensional systems including one (pyrope) with an 80 atom unit cell. CPHF second-order hyperpolarizabilities substantially undershoot experimental values, due to an overestimate of the band gap, but trends are satisfactorily reproduced for beta as well as gamma.

Polarization of one-dimensional periodic systems in a static electric field: sawtooth potential treatment revisited.

Various periodic piecewise linear potentials for extracting the electronic response of an infinite periodic system to a uniform electrostatic field are examined. It is shown that discontinuous potentials, such as the sawtooth, cannot be used for this purpose. Continuous triangular potentials can be successfully employed to determine both even- and odd-order (hyper)polarizabilities, as demonstrated here for the first time, although the permanent dipole moment of the corresponding long finite chain remains out of reach. Moreover, for typical highly polarizable organic systems, the size of the repeated unit has to be much larger than that of the finite system in order to obtain convergence with respect to system size. All results are illustrated both through extensive model calculations and through ab initio calculations on poly- and oligoacetylenes.

The calculation of static polarizabilities of 1-3D periodic compounds. the implementation in the crystal code.

The Coupled Perturbed Hartree-Fock (CPHF) scheme has been implemented in the CRYSTAL06 program, that uses a gaussian type basis set, for systems periodic in 1D (polymers), 2D (slabs), 3D (crystals) and, as a limiting case, 0D (molecules), which enables comparison with molecular codes. CPHF is applied to the calculation of the polarizability alpha of LiF in different aggregation states: finite and infinite chains, slabs, and cubic crystal. Correctness of the computational scheme for the various dimensionalities and its numerical efficiency are confirmed by the correct trend of alpha: alpha for a finite linear chain containing N LiF units with large N tends to the value for the infinite chain, N parallel chains give the slab value when N is sufficiently large, and N superimposed slabs tend to the bulk value. CPHF results compare well with those obtained with a saw-tooth potential approach, previously implemented in CRYSTAL. High numerical accuracy can easily be achieved at relatively low cost, with the same kind of dependence on the computational parameters as for the SCF cycle. Overall, the cost of one component of the dielectric tensor is roughly the same as for the SCF cycle, and it is dominated by the calculation of two-electron four-center integrals.

Coupled perturbed Hartree-Fock for periodic systems: the role of symmetry and related computational aspects.

A general and efficient implementation of the coupled perturbed Hartree-Fock (CPHF) scheme in the CRYSTAL06 code that applies to systems periodic in one dimension (polymers), two dimensions (slabs), three dimensions (crystals) and, as a limiting case, zero dimension (molecules) is presented. The dielectric tensor of large unit cell systems such as boehmite (gamma-AlOOH, 8 atoms/cell), calcite (CaCO3, 10 atoms/cell), and pyrope (Mg3Al2Si3O12, 80 atoms/cell) has been computed. Results are well converged with respect to the computational parameters, in particular, to the number of k points in the reciprocal space and tolerances used in the truncation of the Coulomb and exchange series, showing that the same standard computational conditions used for the self-consistent-field (SCF) step can also be used safely in a CPHF calculation. Point symmetry, being so important in determining crystal properties, also reduces dramatically the computational cost both of the preliminary SCF step and the CPHF calculation, so that the dielectric tensor for large unit cell systems such as pyrope can be computed within 2 CPU hours on a single processor PC.

Calculation of first and second static hyperpolarizabilities of one- to three-dimensional periodic compounds. Implementation in the CRYSTAL code.

A computational scheme for the evaluation of the static first (beta) and second (gamma) hyperpolarizability tensors of systems periodic in 1D (polymers), 2D (slabs), 3D (crystals), and, as a limiting case, 0D (molecules) has been implemented, within the coupled perturbed Hartree-Fock framework (CPHF), in the CRYSTAL code, which uses a Gaussian type basis set. This generalizes to 2D and 3D the work by Bishop et al. (J. Chem. Phys. 114, 7633 (2001)). CPHF is applied for beta and gamma (the polarizability tensor alpha is also reported for completeness) of LiF in different aggregation states: finite and infinite chains, slabs, and cubic crystal. Correctness of the computational scheme and its numerical efficiency are documented by the trend of beta and gamma for increasing dimensionality: for a finite linear chain containing N LiF units, the hyperpolarizability tends to the infinite chain value at large N, N parallel chains give the slab value when N is sufficiently large, and N superimposed slabs tend to the bulk value. High numerical accuracy can be achieved at relatively low cost, with a dependence on the computational parameters similar to that observed for field-free self-consistent field (SCF) calculations.

An ab initio periodic study of NiO supported at the Pd(100) surface. Part 2: The nonstoichiometric Ni3O4 phase.

The present computational study describes the structure and properties of a substoichiometric 2D monatomic in the height phase of nickel oxide, c(4 x 2)-Ni(3)O(4), which has been newly found to epitaxially grow under special deposition conditions on the (100) face of palladium. A slab model is adopted where palladium is simulated by a thin film covered on both sides by epilayers, in combination with a DFT hybrid-exchange Hamiltonian; to make convergence of the SCF procedure easier, a thermal smearing technique is used, whose consequences on the results are critically analyzed. Three adsorbed systems are considered and characterized: (i) RH, that is, the c(4 x 2)-Ni(3)O(4) phase with a rhombic distribution of Ni vacancies, as is experimentally observed; SQ, or p(2 x 2)-Ni(3)O(4), which differs from the previous one for a square, instead of a rhombic distribution of vacancies; (iii) OX, or p(2 x 2)-O, that is, a surface oxidized phase of Pd(100) which is believed to be the precursor for the formation of RH. For a better understanding of the interaction of the metal with the adlayers, the isolated substoichiometric oxides, i-RH and i-SQ, have also been studied. It is shown that RH is more stable than SQ by a few tenths of electronvolts per Ni(3)O(4) unit, which justifies its preferential formation and that the surface reaction, OX + 3NiO(ads) --> RH, is thermodynamically possible. Special attention has been devoted to characterize RH from an energetic, geometric, electronic, and magnetic viewpoint. The strong bond which is formed between surface Pd and O ions in the adlayer is responsible for some peculiar aspects of the electronic and magnetic structures of that phase.