Roberto Dovesi

Fabien Pascale

Bartolomeo Civalleri

Klaus Doll

Nicholas M Harrison

Ian Bush

Philippe D'arco

Yves Noël

Michel Rérat

Philippe D ' Arco

Philippe Carbonnière

Mauro Causà

Simone Salustro

Valentina Lacivita

Bernard Kirtman

Anna Maria Ferrari

Francesco Silvio Gentile

Jacopo Baima

Mauro Ferrero

Raffaella Demichelis

Marco De

La Pierre

and beyond, a long story

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

Quantum mechanical ab initio simulation is rapidly gaining a central role in many scientific communities due to decreasing computational cost, as well as the availability of computer programs of increasing capability and ease of use. In particular, the number of computer codes devoted to periodic systems has been rapidly growing (many of them are presented in this Special Issue; see also, for example, a) Electronic mail: roberto.dovesi@unito.it the book by Evarestov 1 ). CRYSTAL was the first of such codes to be distributed publicly in 1988 through the Quantum Chemistry Program Exchange. [START_REF] Pisani | Hartree-Fock Ab Initio Treatment of Crystalline Systems[END_REF]3 Since then eight other releases have followed in 1992, 1995, 1998, 2003, 2006, 2009, 2014 and 2017 4,5 . CRYSTAL can be used to study the properties of compounds characterized by periodicity in one dimension (quasilinear and helical polymers, nanotubes), two dimensions (monolayers, slabs) and three dimensions (crystals, solid solutions, substitutionally disordered systems). As a limiting case, molecules can also be studied (it is very useful, for example, to investigate molecular crystals and the corresponding isolated molecule by using the same computational setup: basis set, functional, thresholds, computer code). Despite the many improvements and generalizations that have been introduced since its first release, the basic aspects of CRYSTAL have remained the same. Thus, this program computes the electronic structure of periodic systems within the HF and density functional theory (DFT) single particle models using Bloch functions (BF; see section IV). As regards DFT, the local-density and generalized-gradient approximations (LDA and GGA), already available in the early versions of the program, have been complemented with exchange-correlation functionals corresponding to the thirdand fourth-rung of Jacob's ladder [START_REF] Perdew | Jacob's ladder of density functional approximations for the exchange-correlation energy[END_REF] . To the latter rung global and range-separated hybrids of the previous rungs are also available. DFT methods can be augmented with Grimme's D2 [START_REF] Grimme | Semiempirical GGA-type density functional constructed with a long-range dispersion correction[END_REF][START_REF] Civalleri | B3LYP augmented with an empirical dispersion term (B3LYP-D*) as applied to molecular crystals[END_REF] and D3 9 dispersion corrections [START_REF] Grimme | Dispersion-corrected mean-field electronic structure methods[END_REF] . Also, "3c" composite methods (i.e. HF-3c, PBEh-3c, HSE-3c, B97-3c) [START_REF] Caldeweyher | Simplified DFT methods for consistent structures and energies of large systems[END_REF] have been made recently available in the code along with a recent revision for solids [START_REF] Doná | Extending and assessing composite electronic structure methods to the solid state[END_REF] . At variance with respect to other periodic codes, in particular (but not only) based on the use of a plane wave basis set, calculations with hybrid functionals (as well as HF) are only three to five times more expensive than GGA. In CRYSTAL the crystalline orbitals (CO) are expanded as linear combinations of atom-centered Gaussian-type functions (GTFs; also denoted as AOs in the following). Powerful screening techniques are employed to exploit real space locality, which is another distinguishing characteristic of the code. Restricted (closed shell) and unrestricted (spin-polarized) calculations can be performed with all-electron or with valence-only basis sets using effective core pseudo-potentials. Another unique feature is the extensive exploitation of symmetry to achieve computational efficiency. Besides the 230 space groups, 80 two-sided plane groups, 99 rod groups, and 32 crystallographic point groups, there is provision for molecular point group symmetry (e.g. icosahedral) as well as helical symmetry. An extensive discussion of the use of symmetry, in particular for low-dimensional systems, and with particular reference to its implementation in the CRYSTAL code, can be found in the excellent recent book by Evarestov [START_REF] Evarestov | Theoretical Modeling of Inorganic Nanostructures. Symmetry and ab initio Calculations of Nanolayers, Nanotubes and Nanowires[END_REF] . Automatic tools allow users to obtain lower dimensionality systems from higher dimensionality structures by specification of a few geometrical parameters. Thus, slabs (2D periodic), nanorods (1D periodic) and nanocrystals (0D) are easily generated from 3D crystalline structures as well as nanotubes (1D) and fullerenes (0D) from 2D sheets or multilayered slabs. Full use of symmetry is involved in all steps of a CRYSTAL calculation, leading to drastically reduced computation time and required memory, along with improved task farming in parallel computations. [START_REF] Dovesi | On the role of symmetry in the ab initio Hartree-Fock linearcombination-of-atomic-orbitals treatment of periodic systems[END_REF][START_REF] Zicovich-Wilson | On the use of symmetry-adapted crystalline orbitals in SCF-LCAO periodic calculations. I. The construction of the symmetrized orbitals[END_REF][START_REF] Zicovich-Wilson | On the use of symmetry-adapted crystalline orbitals in SCF-LCAO periodic calculations. II. Implementation of the self-consistent-field scheme and examples[END_REF][START_REF] Zicovich-Wilson | On the use of symmetry in SCF calculations. The case of fullerenes and nanotubes[END_REF][START_REF] Orlando | On the full exploitation of symmetry in periodic (as well as molecular) self-consistent-field ab initio calculations[END_REF] Symmetry is also used to select the independent elements of tensor properties to be computed, and plays a crucial role in the selection of dominant configurations in determining the thermodynamic properties of solid solutions, or, more generally, disordered systems. [START_REF] D'arco | Symmetry and random sampling of symmetry independent configurations for the simulation of disordered solids[END_REF][START_REF] Lacivita | On the use of the symmetry-adapted Monte Carlo for an effective sampling of large configuration spaces. The test cases of calcite structured carbonates and melilites[END_REF] A wide variety of crystal properties can be computed automatically, through a simple, single line, keyword. Those that can be classified as total energy derivatives are summarized in table I. The list of properties obtained from wave function analysis is very long and can be found in the CRYSTAL17 Users' Manual [START_REF] Dovesi | CRYSTAL17 User's Manual[END_REF] . Some of them will be presented in section VIII. Many of the new features introduced in the code in the last 8-9 years are described in two papers presenting the 2014 (reference 21) and 2017 (reference 4) releases, to which we refer for an extended description. Technical details about how to run calculations and extract the corresponding information (input/output structure) can be found in the CRYSTAL17 Users' Manual [START_REF] Dovesi | CRYSTAL17 User's Manual[END_REF] , which also contains information on formal aspects. The many tutorials available on the CRYSTAL website 22 are a more direct and agile introduction to the investigation of specific properties. The choice of the Gaussian type basis set is of paramount importance for solids, and requires some attention (see section III). An extended list of basis sets is available on the CRYSTAL website. 22 CRYSTAL may be run either sequentially on a single processor or in parallel. Parallel processing, in turn, can be done either through a replicated data procedure (PCRYSTAL), wherein all the most relevant quantities (i.e. matrices and vectors) are copied on each processor, or according to a Massive Parallel Processing (MPP) strategy in which the large matrices are partitioned and distributed amongst the processors.

While this summarizes at a glance the main features of the code, the rest of the paper is organized into eight further sections, where the most important points are discussed with more details. Section II describes the structure of the code including inner checks of the accuracy and consistency. In section III the choice of basis set is discussed, and the reasons why the use of a standard molecular basis set for crystalline solids can be unbalanced, inefficient and expensive are illustrated. Some background material on HF and on the LDA, GGA and hybrid functionals of DFT, as implemented in CRYSTAL, is presented in Section IV. Then Section V describes the many properties (often tensors of rank 2, 3 and 4) that can be obtained automatically (a single keyword) and, in most cases, at low cost. Two sub-sections are devoted to properties obtained through the CPHF/KS (Coupled Perturbed Hartree-Fock/Kohn Sham) computational scheme and those associated with vibrational degrees of freedom. Section VI deals with symmetry, which plays a crucial role in the efficiency and accuracy of the code. In this area CRYSTAL out-performs most other codes, whether periodic or molecular. An important issue is the limit separating systems that can be tackled with ab initio quantum mechanical methods and

β tuv = -∂ 3 E ∂ ε t ∂ ε u ∂ ε v
Energy E and electric field ; A Raman intensity 3

I a tuv = -∂ 3 E ∂ ε t ε u ∂ r a,v
Energy E, electric field and atomic position vector r a ; A elasticity 4

C tuvw = ∂ 2 E ∂ η tu ∂ η vw
Energy E and rank-2 strain tensor η; A/N photoelasticity 4

p tuvw = ∂ ∆ε -1 tu ∂ η vw
Rank-2 dielectric tensor and rank-2 strain tensor η; A/N second hyperpolarizability 4 γ tuvw = -

∂ 4 E ∂ ε t ∂ ε u ∂ ε v ∂ ε w
Energy E and electric field ; A TABLE I. Physical properties that can be computed with CRYSTAL. For each property, its formula and tensor rank are given along with the general method of computation, which may be either analytical (A) or involve a mixed analytical/numerical (A/N) scheme. Here, t, u, v, w = x, y, z represent Cartesian indices.

those that require approximate methods (force fields, semiempirical Hamiltonians, and so forth). This varies with many parameters (level of ab initio theory, size of basis set, property, etc.). It is a moving target; as the performance of the hardware improves rapidly, so should the software. In section VII we try to define where the border lies with reference to our code. Subsequently, in section VIII, a few examples of the many properties that can be obtained with our code from the wavefunction are reported. The idea is to show that physical understanding of a problem can be achieved and facilitated by the several available interpretation tools. Finally, a few conclusions are drawn in section IX.

II. THE STRUCTURE OF THE CODE AND ITS ACCURACY

The CRYSTAL project was formulated in 1976. The basic computational structure of the code, as already mentioned, is very similar to the one illustrated in three papers published in 1980 (see reference 23 and 24) and 1983 (see reference 25). In 1988 this material was presented more systematically, in connection with the first public release of the code, 3 in a book by Pisani, Dovesi and Roetti 2 (see also reference 26). Fundamental to the early development was the contribution of V.R.Saunders, who joined the project in 1980 and implemented a package for computing two-and one-electron integrals up to d orbitals, afterwards extended to f, and included an elegant procedure for handling the infinite Coulomb series (see below). Subsequent extensions of the code have required efforts from tens of collaborators. As regards the basic structure, briefly: -All relevant integrals (Coulomb, exchange, and monoelec-tronic) are evaluated in direct space. 23,24 -Efficient screening techniques allow for truncation of interactions to a small subset of neighbors of a given atom. 23 -Coulomb infinite series are treated using a combination of ingredients: multipolar expansion of atomic and shell charge distributions in spherical harmonics, 25 together with an Ewald summation of the infinite series that exploits the properties of Gaussian functions, and the grouping of terms to reduce integrals from four to three centers when the distance between interacting electrons is such that the distributions do not penetrate one another. An efficient evaluation of the Coulomb infinite series in a Gaussian basis set for 3D systems was presented in 1992 27 (subsequent proposals do not add much in our opinion). The treatment of 2D and 1D systems was implemented at roughly the same time, although only the latter case was published. 28 -The exchange series is truncated by a screening procedure based on the behavior of the density matrix. 29 This implementation avoids problems that affect other implementations, even if recent developments and the use of screened Coulomb range-separated have mitigated this limitation [30][31][32][33][34][35] .

-Symmetry is used in many steps of the calculation. It has been implemented over many years for one-and two-electron integrals and their storage, [START_REF] Dovesi | On the role of symmetry in the ab initio Hartree-Fock linearcombination-of-atomic-orbitals treatment of periodic systems[END_REF] as well as for construction of Symmetry Adapted Crystalline Orbitals (SACO) that lead to factorization of the Fock matrix [START_REF] Zicovich-Wilson | On the use of symmetry-adapted crystalline orbitals in SCF-LCAO periodic calculations. I. The construction of the symmetrized orbitals[END_REF][START_REF] Zicovich-Wilson | On the use of symmetry-adapted crystalline orbitals in SCF-LCAO periodic calculations. II. Implementation of the self-consistent-field scheme and examples[END_REF] with a consequent major reduction in memory allocation. [START_REF] Orlando | On the full exploitation of symmetry in periodic (as well as molecular) self-consistent-field ab initio calculations[END_REF] -The DFT exchange-correlation contribution to the Fock matrix is evaluated by numerical integration over the unit cell volume. Following a standard numerical protocol 36 radial and angular points for the integration grid are obtained from Gauss-Legendre radial quadrature and Lebedev twodimensional angular distributions respectively.

-Once the Fock (or Kohn-Sham) matrix is built, it is Fourier transformed to reciprocal space and diagonalized. 23 The resulting density matrix is transformed back to real space in anticipation of the next SCF iteration.

-The Unrestricted Hartree-Fock (UHF) option was implemented, by Aprà, Saunders, Harrison in 1992-93, and applied to many transition metal compounds [37][38][39] .

-Analytical gradients with respect to atomic coordinates and cell parameters, needed for geometry optimization, 40 were implemented by Doll et al. in 2000 and following years [41][42][43][44] .

-Vibrational frequencies [START_REF] Pascale | The calculation of the vibrational frequencies of crystalline compounds and its implementation in the CRYSTAL code[END_REF][START_REF] Zicovich-Wilson | Calculation of the vibration frequencies of α-quartz: The effect of Hamiltonian and basis set[END_REF][START_REF] Pascale | Vibration frequencies of Mg 3 Al 2 Si 3 O 12 pyrope. An ab initio study with the CRYSTAL code[END_REF] and second order elastic constants [START_REF] Perger | Ab-initio calculation of elastic constants of crystalline systems with the CRYSTAL code[END_REF][START_REF] Erba | Elastic properties of six silicate garnet end members from accurate ab initio simulations[END_REF] are computed semi-analytically (numerical derivatives of the analytical gradients); -The most recent general documentation of the theory behind the code and of its implementation dates back to 1996 (24 years!), when many contributions from various authors of the code were collected by Cesare Pisani 39 .

A. Accuracy: inner and outer checks

Generally speaking, it is reasonable to assume that the numerical accuracy of periodic codes is the same, in the best cases, or lower, than that of molecular codes, due to the infinite summations and series truncations required for periodic calculations. The accuracy can obviously vary from property to property, but it is difficult to estimate since comparison with the experiment can be biased by the error bar of the latter and by the many approximations implicit in the adopted computational model such as basis set and functional. Cross comparison with other codes is more difficult than in the molecular context, where most of the programs use the same basic ingredients. For example, the large majority of plane-wave (PW) codes use pseudopotentials (PP) that are often unique to that code. One good way to check numerical accuracy (and the formalism, as bugs can always be there!) is to perform inner checks. Three types of inner check will be used below to evaluate the accuracy of CRYSTAL. In subsection II A 1, we show what happens when the size and shape of the unit cell is modified (size and shape consistency), or when the symmetry requirement is reduced/eliminated (symmetry consistency), for a given physical system. In subsection II A 2 we compare a particular property, the infrared (IR) intensity of a large unit cell compound, obtained with three different algorithms and then also with three different parts of the code. A third check consists of a comparison between systems of different dimensionality and verifying that, in the limit of the same physical situation, the numerical values coincide. When starting from the 0D case, the results obtained with molecular codes can be used for the comparison. This scheme has been applied successfully to check the computed polarizability, α 51 , first hyperpolarizability, β , and second hyperpolazizability γ 52 obtained with the CPHF scheme (see section V B). LiF in different aggregation states -finite and infinite chains, slabs, and cubic crystal -was considered. For a finite linear chain containing N LiF units, the (hyper)polarizabilities tend 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. In subsection II A 3 we provide two other examples of this strategy, in which a nanotube is compared with the corresponding monolayer. In that case the compared properties are the vibrational wavenumbers, and IR intensities. When the radius of the tube becomes very large, its properties must tend to the properties of the monolayer. These are also additional examples of size, shape and symmetry consistency, as we compare very large nanotubes, containing 1200 (BN) or 2124 (chrysotile) atoms, and characterized by a given symmetry and shape, to monolayers containing just a few atoms in the unit cell (2 for BN, 18 for chrysotile) with a different shape and symmetry. In all of these cases the agreement turns out to be excellent.

Size, shape, and symmetry consistency checks

In principle, the unit cell of a periodic system can be chosen in many different, but equivalent, ways so that all intensive properties are expected to be the same. This might appear to be a not very useful academic exercise, but that is not the case. We provide here two examples for which the size consistency check is critical. The most common technique for studying defects in solids is the supercell scheme. Consider, for example, a substitutional defect formed by inserting a hetero-atom in a supercell (a multiple of the primitive cell). Thus, the system remains periodic and the effect of the defect concentration may be explored by increasing the size of the supercell. A necessary (but not sufficient) condition for attributing the results obtained with the two supercells to the effect of defect concentration is that both supercells provide exactly the same results (per formula unit) for the perfect host system. A second example is provided by the study of solid solutions. Consider, for example, the Mg-Ca-O system. Again, a supercell approach is usually adopted. Suppose we are interested in a 25% Mg cation composition. The minimum supercell describing this composition is, obviously, MgCa 3 O 4 . But this provides a very poor statistical sample of the many possible configurations. In a supercell 8 times larger, Mg 8 Ca 24 O 32 , many more configurations are possible, although the cost is obviously higher. One possible strategy is to start with small supercells (cheap) and, then, go to increasingly larger supercells (progressively more expensive), checking at each step for convergence of the mean value of the properties of interest. Again, a necessary condition is that the reference perfect system, investigated with different supercells, provides the same results. In table II a set of properties for bulk MgO is evaluated with supercells of different size, shape, and symmetry. The properties are the total energy, the three elastic constants, the electronic contribution to the dielectric constant, along with the fundamental IR mode and its intensity. The first five lines refer to cubic supercells containing 2, 8, 64, 216 and 512 atoms. 10 8 C(1,1,1) 8 48 B3LYP 0.0000 370.3 111.2 189.2 1.71412 6.7553 436.9 4249.1 14 2 7 8 C(1,1,2, The SC column indicates the size and shape of the cell, where P stands for primitive, and C for the conventional cell expanded along the three lattice vectors by i,j,k. C(1,1,1) contains 8 atoms. IS is the shrinking factor, chosen to assure a convergence of the total energy to at least 10 -8 E h . G is the number of symmetry operators. H is the Hamiltonian or functional, i.e. either Hartree-Fock (HF) or hybrid (B3LYP). [START_REF] Becke | Density-functional thermochemistry. III. The role of exact exchange[END_REF] E c is the total energy difference (in micro-Hartree per formula unit) with respect to the conventional cell, whose energy is -1101.72507583 E h . The interatomic distances are the same in all cases (Mg-O = 2.10208181 Å). C 11 , C 12 and C 44 are the elastic constants in GPa. ε ∞ and ε 0 are the optical (high-frequency) and static dielectric constants. ν is the fundamental vibrational frequency in cm -1 and I its IR intensity in km/mol (normalised to the 8 atoms case). ND is the number of SCF plus gradient calculations for constructing the Hessian matrix; n cycle is the number of SCF cycles required to reach convergence. Other computational conditions: grid for the exchange correlation integration XXLGRID (see reference 5); convergence in the energy for the SCF cycle in geometry optimization 10 It is seen that all properties remain essentially the same for the 5 cases. The maximum difference in energy, for example, is always smaller than 0.25 µE h . It can be reduced by one order of magnitude by increasing the parameters controlling the number of direct and reciprocal lattice vectors in the Ewald sums, with negligible cost. It should be noticed that when the convergence threshold on energy TE=10 -11 is used, the vibrational frequency ν of the smallest cell is 6.7 cm -1 higher than that of all other cases. Using TE=10 -12 , however, produces a value in line with the other (super)cells. This higher sensitivity to SCF convergence is due to the fact that, for the cell with only 2 atoms, the full-Hessian matrix is built with a single point (see the ND column in the table), whereas for the other cases at least 2 displacements are used. Lines 6 to 11 in the table show the effect on the various properties of a more severe set of tolerances for the truncation of the Coulomb and exchange series (for more details see the CRYSTAL manual 5 ). Roughly speaking, tolerances act as thresholds for overlaps among Gaussian functions. Thus, T=10 means that overlaps smaller than 10 -10 are disregarded.

There is a cost related to these more severe computational conditions; for example, increasing T from 7 to 10 the cost of the SCF cycle increases for the various MgO supercells by 40%-60%. As regards the results, all properties remain essentially unaltered.

For the shape consistency check, three cases were considered.

Starting from the cubic conventional cell with 8 atoms, the c lattice constant was elongated by a factor 5, 10 and 15, for a total of 40, 80, 120 atoms in the supercell. In this case as well the difference for all reported properties is negligible.

The symmetry consistency check is also quite important. In CRYSTAL, indeed, symmetry is used in many steps of the calculation. In a cubic system, for example, with 48 point symmetry operators, the computation time for the two-electron integrals is reduced by a factor ranging from 20 to 48, according to the case, by the use of symmetry. When comparing results for different symmetry one must be careful to avoid bias due to the different number of symmetry operators. To use again the defect example, it happens frequently that the highest symmetry configuration is not the energy minimum. The broken symmetry case is then investigated. The energy difference, that sometimes is very small, should not depend on the different number of point symmetry operators characterizing the two cases. In table II the 8 and 64 atom cubic supercells, with 48 and 1 (C 1 ) symmetry operators respectively, were investigated. The largest difference along all properties occurs for the IR intensity, which varies by less than 0.3% (4244 vs 4255 km/mol) for the B3LYP functional.

The numerical noise can also depend on the adopted Hamiltonian. DFT results might be affected by the choice of the numerical grid adopted for the integration of the exchangecorrelation functional; Hartree-Fock (HF), by the truncation of the exchange series. In the table the HF results for 3 supercells, S 2 , S 8 , and S 64 , are reported. The differences are even smaller than the ones between the corresponding three B3LYP cases (for example E c for the supercell with 64 atoms is 1.35 • 10 -8 E h for B3LYP, and only 4 • 10 -10 E h for HF ), indicating that part of the residual noise in the latter case is attributable to the numerical integration grid. As a final comment, it should be noticed that, in most cases, standard computational conditions (T=7) were used for the above calculations. Working with more severe tolerances (T=8 or 9) drastically reduces the very few minor anomalies.

Using alternative algorithms

A completely different check on the accuracy of the code consists in evaluating the same property through alternative algorithms. Since each alternative involves, for the most part, quite different ingredients (quantities in direct or reciprocal space, obtained from Fourier or anti-Fourier transforms, analytical or numerical derivatives, an SCF process or not), finding the same value represents a check (although not complete) of the accuracy of the algorithms and of their implementation. An example is the IR intensity of the vibrational modes, which can be evaluated in three different ways [START_REF] Dovesi | Calculation of the infrared intensity of crystalline systems. A comparison of three strategies based on Berry phase, Wannier function, and coupledperturbed Kohn-Sham methods[END_REF] , that probe three quite different paths through the CRYSTAL code.

• One algorithm exploits the Berry phase (BP) theory due to King-Smith and Vanderbilt, [START_REF] King-Smith | Theory of polarization of crystalline solids[END_REF] and Resta [START_REF] Resta | Quantum-mechanical position operator in extended systems[END_REF] to evaluate the dipole moment difference between the equilibrium geometry and each one of the distorted geometries used for numerical differentiation. [START_REF] Dall'olio | Spontaneous polarization as a Berry phase of the Hartree-Fock wavefunction: The case of KNbO 3[END_REF][START_REF] Noël | Polarization properties of ZnO and BeO : An ab-initio study through the Berry phase and Wannier functions approaches[END_REF] • A second scheme utilizes localized Wannier Functions (WFs) to obtain the dipole moment difference from the sum of the reference WF centroids. [START_REF] Zicovich-Wilson | A general method to obtain well-localized Wannier functions for composite energy bands in linear combination of atomic orbital periodic calculations[END_REF][START_REF] Zicovich-Wilson | Beyond Standard Quantum Chemistry: Applications from Gas to Condensed Phases[END_REF][START_REF] Zicovich-Wilson | Ab initio simulation of the IR spectra of pyrope, grossular, and andradite[END_REF] The geometries and numerical differentiation procedure are the same as in the BP method.

• The third method employs the Coupled Perturbed Hartree-Fock (CPHF) scheme of Hurst et al. [START_REF] Hurst | Ab initio analytic polarizability, first and second hyperpolarizabilities of large conjugated organic molecules: Applications to polyenes C 4 H 6 to C 22 H 24[END_REF] , as adapted for electric fields in periodic systems and for Coupled Perturbed Kohn-Sham DFT. [START_REF] Ferrero | The calculation of static polarizabilities of 1-3D periodic compounds. The implementation in the CRYSTAL code[END_REF][62][START_REF] Otto | Calculation of the polarizability and hyperpolarizabilities of periodic quasi-one-dimensional systems[END_REF][START_REF] Kirtman | Extension of the Genkin and Mednis treatment for dynamic polarizabilities and hyperpolarizabilities of infinite periodic systems. I. Coupled perturbed Hartree-Fock theory[END_REF][65] This is done in combination with the standard analytical treatment of geometric energy gradients, based on the orbital energy weighted density matrix, to yield the desired second derivatives of the energy with respect to the field and atomic displacements. 66,67 An application of the three algorithms to pyrope (four units of Mg 3 Al 2 (SiO 4 ) 3 , 80 atoms/cell, space group Ia3d, 48 symmetry operators) is reported in table III. Here the B3LYP functional [START_REF] Becke | Density-functional thermochemistry. III. The role of exact exchange[END_REF]68 was used together with the same basis set as in references 60 and 47: an 8-511G* contraction for Mg, 8-611G* for Al, 8-6311G* for Si and 8-411G* for O (1488 contracted AOs in the unit cell). The table shows that the results obtained with the three schemes are very similar. The largest difference is 57 km mol. -1 , for mode 15 at 865 cm -1 , which is however a difference of only 0.4 %. The percentage difference increases up to 4.5% for mode 7, but in that case the intensity is relatively low, 43 or 45 km mol. -1 , and the absolute difference is only 2 km mol. -1 . It should be stressed that these results were obtained under standard computational conditions. By increasing the accuracy requirements a bit further, the differences can be reduced by an order of magnitude. As regards the cost of the three schemes, we observe that: -The cost of the calculation of the IR intensity is always a small fraction of the cost required by the construction of the Hessian matrix.

-The cost of the BP and CP schemes is about the same, whereas the WF scheme is more expensive. -These figures reflect, more than the ideal efficiency of the schemes, their story from the first implementation.

-The Berry Phase approach has been the first to be implemented (1997). It is the simplest of the three: it does not rely on iterative procedures, involves only little linear algebra, and does not require additional Fock matrix builds. It has been ported to a massively parallel implementation already in 2014 and applied to the study of large systems, up to an order of magnitude larger than pyrope.

-The parallelization of the WF and CP schemes is not complete yet and as efficient as the BP one.

-As a final comment, we observe that the BP and WF schemes, with a calculation for each configuration explored when building the Hessian matrix (displacements along x,y and z for each atom) can be used in conjunction with the fragment strategy (see Section V A 3), that is mandatory when the IR spectrum of very large systems (containing, say, more than 1000-2000 atoms) is required.

When changing the dimensionality of the system

Consider the properties of a set of (n, 0) BN nanotubes. As n increases the curvature of the tube decreases and, in the limit of extremely large n the properties of the tube must be very close to those of the monolayer. Let us recall that the unit cell of the monolayer contains only two atoms and 3N -3 = 3 fundamental modes of vibration. They are the A mode at 836.7 cm -1 , with IR intensity, I = 13.07 km mol. -1 , and the doubly degenerate E mode, at 1376 cm -1 with I = 2340 km mol. -1 . In the nanotube there are 12n -4 fundamental modes of vibration since there are 4n atoms in the unit cell and there is a degree of freedom corresponding to the rotation around its axis. For very large n the effect of the curvature is negligible. Then, in the case of a (12, 0) nanotube, for example, the 48 atoms in the unit cell correspond to the unit cell of the slab, plus 23 copies of this cell; that is, to a supercell in which the two lattice vectors of the 2D cell have been multiplied by 2 and 12. The set of eigenvalues of the tube will correspond to the eigenvalues of the monolayer at Γ plus a set of 23 other k-points of the first Brillouin zone. In this limit, -only the modes corresponding to those of the monolayer (A and E ) have IR and Raman intensities different from zero; the modes at any k = Γ are inactive for symmetry reasons; -the frequency of three modes is going to zero, and its derivative as a function of n can be related to the elastic constants of the monolayer in directions orthogonal to the tube axis (by convention oriented along the x direction).

For the calculations, the B3LYP functional and a split valence basis set has been used. Let us consider, first, the two modes with the highest wavenumbers, that tend towards the highest frequency mode of the monolayer at 1376 cm -1 . Figure 1 shows the behavior of their wavenumbers as a function of n.

The A-symmetry case (atomic displacement along the periodic direction of the tube) converges rapidly to the slab limit.

The difference in wavenumber is only 1 cm -1 for n = 30. The intensity, shown in the second panel, is also very close to that of the monolayer; for n = 30 it is 2382 km mol.

-1 compared to 2340 for the slab.

The situation is completely different for the E-symmetry modes (in-plane modes along y and z, the two directions orthogonal to the tube axis on each half of the circumference of the nanotube). The frequency at the infinite n limit for these modes must also approach the monolayer value of 1376 cm -1 . However, for the (30, 0) nanotube the frequency is 1474 cm -1 and the intensity is 1191 km mol. -1 , as compared to the target value of 2340 km mol. -1 . When n = 150, the wavenumber is reduced to 1404 cm -1 (within 28 cm -1 of the target) and the intensity is 2010 km mol. -1 . For the (300,0) nanotube the frequency and intensity differ by only 14 cm -1 and 160 km mol. -1 , respectively, from the slab values.

The same behavior with respect to n is obtained for the Aand E-symmetry modes around 830 cm -1 (not shown) corresponding to the A out-of-plane, and much less intense, mode of the slab. The connection between monolayer and nanotubes is more delicate and less evident for the nanotube modes that tend to zero. In that case it is useful to represent the tube wavenumbers as a function of the inverse of n (see figure 2). It can be shown that the slope of these curves is related to the elastic components of the monolayer (i.e. C 11 and C 66 ). 69,70 For example, for the A 1 -mode

ω A 1 = C 11 (M B + M N ) 2π a 1 n (1)
where a is the lattice parameter of the monolayer, and M B and M N are the atomic masses of B and N . The calculated value of C 11 (= C 22 ) of the slab is 3.755 E h leading to ω A 1 = 2656/n cm -1 which is in good agreement with the fitted value of 2644 ± 15 cm -1 of figure 2 for the A 1 mode. The connection between C 66 = 1 2 (C 11 -C 12 ) = 1.480 E h and the E 1 soft mode of the nanotube leads to ω E 1 = 1667/n cm -1 , to be compared to 1659 ± 15 cm -1 from the fitting. Finally, the relationship for the E 2 soft mode is similar to that for the A 1 mode, but with half the total mass in Eq.(1) so that ω E 2 =

√

2ω A 1 = 3756/n cm -1 , which is to be compared with 3745 ± 15 cm -1 from the fitting. The BN case might be considered too simple since the monolayer is planar, the unit cell contains only two atoms, and there are many symmetry constraints involved in determining the 1D-2D correspondence. Thus, a second case was considered, namely the lizardite 2D monolayer (18 atoms in the unit cell) that corresponds to the infinite n limit of a chrysotile nanotube with an ideal formula unit corresponding to Mg 3 Si 2 O 5 (OH) 4 . The crystal structure of lizardite consists of so-called TO layers, where a brucite-type octahedral sheet (O) is attached to a silica tetrahedral sheet (T). The latter consists of di-trigonal rings formed by SiO 4 tetrahedra units through sharing one oxygen atom (an apical oxygen, O a ). What makes their structures so peculiar and difficult to determine is the bending flexibility of these slabs, which can give rise to a range of layered structures where the layers can be flat (like in lizardite), or wavy (like in antigorite), or curled in the form of concentric spirally wrapped nanotubes (chrysotile). 71,72 In figure 3 a chrysotile tube, two lizardite layers seen along two orthogonal views, and an antigorite bilayer are shown. A family of nanotubes from (9,-9) to (59,-59) was considered; the unit cell of the latter contains 2124 atoms. The B3LYP functional and a split valence basis set (see Ref. 73) were used. For the comparison of vibrational wavenumbers, we selected the two most intense Raman peaks of lizardite, at 702.1 and 3800.9 cm -1 , and examined how two frequencies of the nanotubes tend to this limit when the radius increases. The Raman spectrum of the nanotubes presents only two modes with important intensity for large radii, as expected. These modes are indicated as A 1 -large and A 1 -small in figure 4, which reports the variation of the wavenumber of the two peaks as a function of n. When fitted with the function ν = a + b/n + c/n 2 + d/n 3 , the limiting values shown in the figure are 700.67 ± 0.21 and 3798.17 ± 0.11 and differ from the slab values by only 1.5 and 3 cm -1 respectively. These results show that a) the accuracy of the frequency calculation is very similar for unit cells containing either a few (18) or more than two thousand atoms, and (b) the accuracy is the same for both 1D and 2D systems. It is interesting to note that for the (59,-59) nanotube, the high frequency mode is much further from the monolayer (11 cm -1 ) than the low frequency mode (0.2 cm -1 ) (note the different scales at the right and left hand sides of the figure). The larger difference is a measure of the residual interaction in the nanotube due to the small but not null curvature, that generates a red shift for the hydrogen stretching motion.

The chrysotile example will also be used in Section VI, showing the effect of symmetry, that is particularly high for nanotubes, and in Section VII, devoted to the code performance, because the geometry of very large tubes 73 containing up to 5000 atoms has been optimized 74 for investigating the relative stability of lizardite and chrysotile tubes, and because the vibrational frequencies and the IR spectrum have been computed for tubes containing up to 2124 atoms thanks to symmetry.

III. THE CHOICE OF THE BASIS SET.

A peculiar feature of CRYSTAL that makes it unique amongst other codes, is the use of atom centered Gaussiantype basis functions. The choice of the Gaussian basis set is not only crucial to determine the accuracy of the results, but also a delicate matter for tightly packed crystalline systems. In the naive strategy, frequently adopted, molecular basis sets are used as such. While for small-to-medium basis sets this might be not a problem (for example when used for molecular crystals) extended molecular basis sets often include diffuse Gaussians with exponents as small as, say, 0.05-0.02 Bohr -2 needed to describe the tails of the wavefunction. Several diffuse Gaussians are necessary to recover the correct long-range exponential behavior in the external region. In contrast with a molecule, however, a tightly packed 3D system like diamond or MgO does not have any empty space external to the system and, therefore, such functions are not needed either ex- ternally or internally 24,75,76 . Moreover, the Bloch functions constructed from diffuse Gaussians tend to be very similar. This leads to an overlap matrix that is nearly singular which, in turn, creates serious numerical problems with regard to orthogonalization of the basis functions.

The major advantage of a local basis set is that matrices remain small and the number of integrals limited. In the following we document how it, then, becomes possible to calculate the properties of very large unit cell systems at low cost. This feature is lost progressively when the exponents of the most diffuse Gaussian-type functions decrease. It is why specific basis sets have been implemented for periodic systems, both by ourselves (see the CRYSTAL website) and other groups. [77][78][79] For those circumstances where diffuse functions are desirable one may be able to add them subsequently by means of perturbation theory, as suggested in a very recent paper by Maschio and Kirtman. 80

IV. HARTREE-FOCK AND DFT

As already mentioned in the introduction, CRYSTAL was formulated for performing Hartree-Fock (HF) calculations of solids, as opposed to the various codes that were appearing in these years, all of them based on LDA, with many other mixed ingredients. All the papers published by the group for about 15 years, including covalent (graphite), 23,24 ionic (MgO and the relative stability of its phases), 81 semi-ionic systems like quartz, 82 molecular crystals, 83 defects treated with the supercell scheme, 84 transition metal compounds, 37,85,86 as well as many properties, were based on HF. There are many reasons for this choice.

-The molecular, or quantum chemistry, background of Pisani, Dovesi, Roetti and Saunders.

-At that time (1970-80) most of the attention of the solid state community (physicists) was focused on the band structure and density of states, not on total energy and related properties (formation energy, equilibrium geometry, bulk modulus), that were on the contrary the key outcomes of the calculation for a chemist. The accurate calculation of total energy was still a challenge for most of the codes.

-At the beginning of the seventies, Euwema et al. [START_REF] Surratt | Hartree-Fock lattice constant and bulk modulus of diamond[END_REF][START_REF] Wilhite | Charge-conserving integral approximations for ab initio quantum chemistry[END_REF] implemented a code for 3D solids (applications to diamond, BN and others), using many of the ingredients (HF, Gaussian basis) that Pisani, Dovesi and Roetti were planning to use, showing that this kind of calculations was possible, although difficult.

-Probably the main reason for this choice (HF) was the belief that in any case HF could represent a valid solution by itself, and a good starting point for post HF treatments (see the Colle-Salvetti correction below). This opinion was quite unconventional in the solid state community. One of the authors (RD) was invited in 1984 to a workshop, the attendance of which was LDA-PW-PP based. When he presented results for bulk silicon, most of the audience was surprised to know that HF could produce a reasonable total energy, equilibrium geometry, band structure, charge density maps and Compton profiles, because it was a widespread opinion that the pathological behavior of HF at the Fermi level of the electron gas would just represent an example of the inadequacy of HF in describing any system and any property. So the final comment of the chairman was: why to do expensive calculations that are completely wrong, instead of using LDA that produces excellent results at much lower cost? We will come back to this point in section IV B, where we will try to demonstrate that HF can be very useful in some cases, in spite of the opinion of the chairman. We will also show that, as implemented in CRYSTAL, HF is not much expensive, at least for ionic and covalent compounds, for which it performs best. CRYSTAL authors were well aware of the limitations of the pure HF treatment for many properties of crystalline solids. Note, however, that these were not larger than the ones affecting in these years the results produced with other codes, due to the adopted functional (LDA), to problems with the description of core electrons, to basis set limitations, to continuity problems of the wavefunction at the frontiers between zones in which the cell volume was partitioned in some of the schemes in use at that time. To try to overcome the inaccuracies of HF, an a posteriori correction to the HF energy according to a scheme formulated in 1975 by Colle and Salvetti 89 was adopted. The implementation dates back to 1987 [START_REF] Causà | Correlation correction to the Hartree-Fock total energy of solids[END_REF] with a further extension in 1988. [START_REF] Causà | Correlation correction to the Hartree-Fock total energy of solids. II[END_REF][START_REF] Causà | Erratum: Correlation correction to the Hartree-Fock total energy of solids[END_REF] DFT schemes were inserted in the self consistent structure of the code as an alternative to HF a few years later, mainly by Causá, Zupan, Towler, in the LDA and GGA version. [START_REF] Causà | Density functional LCAO calculation of periodic systems. A posteriori correction of the Hartree-Fock energy of covalent and ionic crystals[END_REF][START_REF] Zupan | Density functional LCAO calculations for solids: Comparison among Hartree-Fock, DFT local density approximation, and DFT generalized gradient approximation structural properties[END_REF][START_REF] Towler | Density functional theory in periodic systems using local Gaussian basis sets[END_REF][START_REF] Zupan | Density-gradient analysis for density functional theory: Application to atoms[END_REF] The XC (exchange correlation) matrix elements were evaluated by best fit, representing the XC potential as a set of Gaussian functions. This part of the code became very efficient when V. R. Saunders implemented a package in which the fitting set was extended up to I functions (l=6), and the integrals were evaluated by exploiting recursion relations (unpublished). In 2000-2001, the integration technique for the XC matrix elements was changed: the fitting strategy was substituted by the numerical integration technique based on Becke's scheme, 36 quite effective and accurate. The latter is based on an atomic partition method, in which the radial and angular points are obtained from a Gauss-Legendre quadrature and a Lebedev two-dimensional distribution, respectively. The choice of a suitable grid is crucial for both numerical accuracy and cost. A pruned grid (the number of angular points reduces with the distance from the nuclear position) with 75 radial and 974 angular points is usually employed. The reason for changing the strategy was due to the difficulty in reaching high accuracy with the fitting scheme, in particular with all electron basis sets (linear dependence problems when a too large set of fitting functions was used), and a reasonable scaling with large unit cells (say more than 1000 atoms). For a few years then it has been possible both to perform the SCF with various functionals of the GGA family, and the a posteriori correction for the correlation energy. Among the various attempts, an SCF HF scheme supplemented by a correlation only functional was used for a while.

The limitations of LDA and GGA became evident in the eighties, and many proposals mixing HF and DFT exchange were formulated. In their paper: Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density, Lee, Yang and Parr (LYP, 1988) 68 provided the background for the formulation, in 1993, of the B3LYP hybrid functional, [START_REF] Becke | Density-functional thermochemistry. III. The role of exact exchange[END_REF] one of the most cited and used, at least in the molecular community. Many of the combinations mixing the various ingredients (exchange and correlation, Hartree-Fock and LDA or GGA) where rapidly incorporated in the code. Later, mGGA functionals were also implemented along with range-separated hybrids (CRYSTAL14). Therefore, the last version of the code incorporates exchange-correlation functionals up to the forth-rung of the Jacob's ladder. Definitely, hybrid functionals were available in CRYSTAL since mid-nineties [START_REF] Towler | Density functional theory in periodic systems using local Gaussian basis sets[END_REF] (i.e. a few years after the seminal work by Becke [START_REF] Becke | Density-functional thermochemistry. III. The role of exact exchange[END_REF] ) and used for many investigations in the next years, whereas in the PW-PP community the migration from LDA to GGA was slow. But the real challenge was the implementation of hybrid functionals, whose cost was orders of magnitude higher than for pure DFT, when used with a PW basis. A partial solution to this problem arose when range separated hybrids appeared on the market, and hybrids became formally available in many public PW-PP codes, although the number of systems investigated at this level remains a small fraction of the total, suggesting that the cost is still high with respect to pure DFT, or possibly that technical problems persist.

A. On the choice of the XC functional

There is now a widespread consensus on the opinion that mGGA and GGA are better than LDA, and hybrids are better than mGGA and GGA, although the evidence is sometimes not strong, and the number of exceptions and counterexamples is high. For instance, some of the present authors in collaboration with A. Savin have proposed that presently the selection of a functional is largely a personal choice and the relative performance of different approximations can depend on different statistical measures [START_REF] Civalleri | On choosing the best density functional approximation[END_REF][START_REF] Pernot | Prediction uncertainty of density functional approximations for properties of crystals with cubic symmetry[END_REF][START_REF] Civalleri | On the use of benchmarks for multiple properties[END_REF] . In spite of this, during the last decades, the inclusion of HF exchange in hybrid functionals has been proved to be crucial for many properties from molecules to solids, so that in many cases they are the methods of choice. As mentioned, CRYSTAL offers a computationally efficient implementation of different hybrid functionals, so that their performance can be easily tested. Therefore, when investigating a new class of problems, the attitude of many of the present authors is to compare various functionals, most of them being hybrids, for instance: B3LYP, PBE0 and HSE06, but also PBE and LDA and finally HF, knowing that HF on one side, and LDA and PBE on the other side, define the extreme values, and hybrids are in between and in general much closer to what should be considered the correct result. While there could be some differences between the hybrid functionals, mostly due to the way (and amount) HF exchange is incorporated, many of those tested achieve a useful accuracy. On the basis of previous applications and general experience, some of the present authors use B3LYP as a first choice. More than 25 years of accumulated experience brought to a very rich statistical file, that has allowed us to know limits and merits of the functional in the different situations. This is why in the present work most of the reported results have been obtained with the B3LYP functional.

B. Hartree-Fock and the superexchange interactions in Transition Metal compounds

When in 1992-93 the UHF (as well as the Restricted Open Shell Hartree Fock, ROHF) option was implemented, it was applied to a set of transition metal (TM) oxides, (VO, MnO, NiO) 37,38 , and, in a few years, to many other TM ionic compounds: Fe 2 O 3 hematite 100 , KNiF . Surprisingly, the description of the spin density and the energetics of the various spin states was reasonably accurate. The many spin states, that in TM molecular complexes are often separated by tiny energies, on which electron correlation plays a dramatic role, are on the contrary properly described in solid state compounds. Two factors simplify the scenario in solid state: a) the ionic character of many of these compounds, that increases the energy differences among spin states, b) the high symmetry, that imposes orthogonalities reducing the number of possible states and again separates the levels. In these years the main tool for describing these compounds, when approached with quantum mechanical periodic codes, was through the LDA functional plus the so called scissor operator 108,109 . Many of these systems, say for example NiO, that are large gap insulators, when treated at the LDA level turn out to be metallic, and the spin localization on Ni is spread over many atoms. The scissor operator opens the gap, maintains artificially the insulating character of the system, and forces somehow the localization of the uncoupled electrons. The amount of imposed separation between valence and conduction bands is usually a function of the target results. A second, slightly more sophisticated version is LDA+U, where U stands for Hubbard 110,111 . The U additional term plays the same role as the scissor operator; in principle the numerical values used for it could be obtained from atomic data 112 . In practice, in many cases, it was used as a parameter for reproducing the experimental results. It should be underlined that the structural properties of these compounds, like the lattice parameter, the bulk modulus, the elastic tensor and even the formation energy, are not very sensitive to the degree of localization of the unpaired electrons, or to the quality of the description of the short range exchange interaction. Other properties, on the contrary, like the superexchange coupling constants, are dramatically dependent on the quality of this description. The superexchange coupling constant, usually indicated with J (for historical reasons measured in Kelvin), is the quantity that is tabulated by experimentalists, and that can be considered a measure of the energy difference between the ferromagnetic (FM) and antiferromagnetic (AFM) solution.

The usual way the calculated data are compared with experiment is through an Ising model, in which the interaction between magnetic centers is considered additive. As this interaction is falling down exponentially with the distance between the magnetic centers, in many cases it is sufficient to consider the interaction between first magnetic neighbors.

∆E = 2zS 2 | J| ( 2
)
where z is the number of magnetic first neighbors with opposite spin. Consider for example the KMF 3 compounds, where M can be Mn, Fe, Co, Ni and Cu. The ideal cubic structure (see Figure 5) is perturbed by the Jahn-Teller distortion in the Fe, Co and Cu cases, whereas it is preserved in KMnF 3 and KNiF 3 , thanks to the d 5 and d 8 electronic configuration of the TM, that allows to fill completely the t 2g and e g α levels in Mn, and also the t 2g β levels in Ni. The structure is a simple cubic array of MF 6 octahedra. The K + ions fill the empty space (dodecahedra) between the octahedra (see Figure 5, left). The M-F distance is a/2 (a is the lattice parameter).

There are six M second nearest neighbors of the central M atom. The M-M distance is a. At low temperature, the systems are AFM. The AFM structure consists of ( 111) sheets of M atoms of common spin, the spin alternating between the sheets (see the central panel). As a consequence, each M atom is surrounded by 6 second nearest M atoms of opposite spin. For obtaining ∆E in Eq. 2, the unit cell is doubled, and the HF-SCF calculation is performed with (say for the Mn compound) the 5 electrons up on both atoms (FM solution), and with the spin up on one atom and down in the second in a second run (AFM solution). The energy difference, divided by two (we are interested in the energy difference per magnetic center), is shown in Table IV B, first column. Thus, -The PBE0, B3LYP and PBE vallues are about factors 1.5-5 higher for J and 5-13 for J y (see below), where the ratio is similar for Mn and Ni. It should be underlined that this qualitative agreement provided by HF remains valid in many situations: -when a second star of magnetic neighbors is considered, as shown in Ref. 86 for the KMF 3 set.

| J| = ∆E 2zS 2 (3
-when the set of six first magnetic neighbors splits in two sets, as in KFeF 3 and KCoF 3 -when the situation becomes even more complicated, as in KCuF 3 , where both FM and AFM interactions appear 85 .

It should be underlined that this kind of comparison between experiment and simulation is delicate, as it can be biased by the models used in going from the raw experimental data to the tabulated J and from the calculated total energies to the calculated J. For example, the calculations are performed at the experimental geometry, in order to eliminate part of the error due to the different (experimental vs calculated) lattice parameter (remember that J fall down exponentially). The main limitation on the side of the calculation is however the following. The UHF solution is an eigenfunction of the S z operator, but not of S 2 . Closed shell solutions, and the highest S z UHF solutions are also eigenvalues of S 2 . In the other cases, the S z solutions are a mixing of different S 2 states, so that one must be careful when comparing experimental and calculated data. In DFT the situation is less clear, but in general to the UKS (Unrestricted Kohn Sham) solutions the same meaning is attributed as to the UHF wavefunctions.

Recently, in CRYSTAL it has been implemented an option for evaluating the degree of contamination of the spin state. The S 2 H and S 2 L columns of Table IV B report the mean value of the S 2 operator in the high (H, FM solution) spin state (S z =2 and S z =5 for KNiF 3 , and KMnF 3 , remember that the unit cell contains two atoms). As expected, the high spin mean value is close to the theoretical value (6 for the Ni compound, 30 for the Mn one). For the low spin state, with S z =0, the S 2 L value should be zero, which is not, due to contamination. The S 2

L value is very close to S z (2 and 5 for the Ni and Mn compound, respectively). Yamaguchi et al. 113 proposed a formulation for obtaining J, indicated as Jy and alternative to Eq. 3, that somehow takes into account the spin contamination in the low spin state: the new data are shown in the last column of table IV B. It turns out that now the agreement with experiment is not only qualitative, but also quantitative. Using this same formula, the PBE, B3LYP and PBE0 data are even farther from experiment than when using Eq. 3. HF then seems to produce, ab initio, that is without any parametric arrangement as when using the scissor operator or LDA+U, quite reasonable superexchange coupling constants. A more systematic investigation must be performed to confirm the excellent results reported in Table IV B. A comment can be added here concerning a CRYSTAL feature that we can call time consistency. Calculations for KMnF 3 and KNiF 3 were performed in 1997 107 , obtaining for ∆E 0.293 and 0.569 mE h , respectively, to be compared with 0.294 and 0.569 mE h , obtained when preparing this manuscript and appearing in the first column of Table IV B. So the same calculation repeated after 23 years produces the same total energy, with just a µE h difference. The stability of the results for so many years is probably not shared by many other periodic codes.

| Jy | = ∆E S 2 H -S 2 L (4) System Method ∆E(mE h ) S 2 H S 2 L J(K) Jy (K) KMnF 3 HF 0.

V. ENERGY DERIVATIVES AND TENSORS

In the investigation of crystalline solids, a large fraction of the literature is still focused on relatively simple and low-cost properties, like the band gap and the density of states, for which a single SCF calculation is required. Nevertheless, the interest and computational effort is now moving towards more complex (and more interesting, in the authors' opinion) observables, such as the physical tensors (elastic, dielectric, piezoelectric, photoelastic, hyperpolarizability, and many others) and the vibrational properties (frequencies, Infrared and Raman intensities and spectra, thermodynamics). All these properties are related to derivatives of the total energy with respect to four variables, namely the coordinates of the atoms in the unit cell (r a in table I), the lattice parameters (a i ), one (or more) external electric fields and one or more magnetic fields (not implemented in CRYSTAL yet). The first two lines of table I refer to the gradient, both with respect to the Cartesian atomic coordinates 41,42 and the cell parameters. 43,44 They are the basic ingredients for the geometry optimization 40 and transition state search 116 . The former is very flexible and accurate (structures containing more than thousand atoms in the unit cell have recently been optimized without any particular problem, as for the cases with a few tenths of atoms). Several algorithms rely on the high numerical accuracy of the geometry optimizer: search for equilibrium structures and transition states, volume-or pressure-constrained minimizations 117 for the determination of the equation of state of bulk crystals, semi-analytical determination of the vibrational frequencies, nuclear relaxation of strained lattices for the computation of elastic, piezoelectric, photoelastic tensors, etc. In table I examples of computed tensors of rank 2 or 3 or 4 are shown. To this latter set belongs the elastic tensor (related properties are the seismic wave velocities), [START_REF] Perger | Ab-initio calculation of elastic constants of crystalline systems with the CRYSTAL code[END_REF][START_REF] Erba | Elastic properties of six silicate garnet end members from accurate ab initio simulations[END_REF]69 the photoelastic Pockels' tensor, 118 the second order hyperpolarizability. 119 The direct and converse piezoelectricity, 120 as well as the first hyperpolarizability, [START_REF] Ferrero | Calculation of first and second static hyperpolarizabilities of one-to three-dimensional periodic compounds. Implementation in the CRYSTAL code[END_REF]121 are third order tensors. The dielectric or polarizability is a second order tensor, and is evaluated, as its higher order terms, analytically via the Coupled Perturbed Hartree-Fock/Kohn-Sham (CPHF/KS) method. [START_REF] Ferrero | The calculation of static polarizabilities of 1-3D periodic compounds. The implementation in the CRYSTAL code[END_REF][START_REF] Ferrero | Calculation of first and second static hyperpolarizabilities of one-to three-dimensional periodic compounds. Implementation in the CRYSTAL code[END_REF]65,119,121,122 The last column indicates that some of these third or fourth order tensors are computed semianalytically: this means that a second derivative is computed as a numerical derivative of an analytical derivative. It must be underlined that in all cases: a) a simple keyword is sufficient for the calculation of the full tensor. For example, for the elastic, piezoelectric and photoelastic tensors the input lines reduce to a single keyword: ELASTCON, or PIEZOCON, or PHOTOELA.

b) All the components of the tensor (they are, for the elastic tensor, 3 and 21 for a cubic or triclinic system, respectively) are obtained with a single calculation. This is not the case of the experiment. Many different measurements are required; each one of them produces a linear combination of components of the tensor; a system of linear equations must be solved for determining the individual elements of the tensor. As usually small and large numbers appear, the accurate experimental determination of the small components is extremely difficult; this difficulty is increased by the fact that in many cases the sign of the components is undetermined. Two sets of properties listed in table I require some additional comment, for their generality and importance. The first one includes the Hessian matrix, that produces the vibrational frequencies, and the IR and Raman intensities; they permit a complete characterization of the vibrational related properties, including the thermodynamics. The second set includes all properties that can be obtained by the SC-CP (self-consistent coupled perturbed scheme, a very general, analytical, versatile, tool).

A. Vibrational frequencies, infrared and Raman spectra

Vibrational frequencies and related properties are available in CRYSTAL since 15 years (CRYSTAL06). Frequencies at the Γ point are obtained within the harmonic approximation by diagonalising the mass-weighted Hessian matrix, W , whose elements are defined as , [START_REF] Pascale | The calculation of the vibrational frequencies of crystalline compounds and its implementation in the CRYSTAL code[END_REF][START_REF] Zicovich-Wilson | Calculation of the vibration frequencies of α-quartz: The effect of Hamiltonian and basis set[END_REF]123,124 :

W Γ αi,β j = H 0 αi,β j M α M β with H 0 αi,β j = ∂ 2 E ∂ u 0 αi ∂ u 0 β j , (5)
where M α and M β are the masses of the atoms associated with the i and j coordinates of atoms α and β . Energy first derivatives with respect to the atomic positions, g α, j = ∂ E ∂ u α, j , are calculated analytically for all the u α, j coordinates (E is the total energy, u α, j is the displacement coordinate with respect to equilibrium). Second derivatives at u = 0 are calculated numerically using a single displacement along each coordinate (N=2, the central point and a point on the positive direction of the coordinate):

∂ g α j ∂ u β i ≈ g α j (0, . . . , u β i , . . . ) u β i (6)
or averaging two displacements (N=3):

∂ g α j ∂ u β i ≈ g α j (0, . . . , u β i , . . . ) -g α j (0, . . . , -u β i , . . . ) 2u β i (7)
If the number of atoms is n, then (3n-3)*1 or (3n-3)*2 SCF+gradient (G) calculations are performed, for N=2 or N=3, respectively. The default value for N is 2. The default value for the step is 0.003 Å, much smaller that in other codes. It is possible to use such a small step because the accuracy of both energy and gradient is very high, as high must be the accuracy in the definition of the equilibrium position. A value of 0.001 Å can be safely used in many cases. However, when very soft modes are present in the spectrum, 0.003 Å for the step, and N=3 are to be preferred. The accuracy in this case is higher than with N=2, because the two displaced points along each coordinate are symmetric with respect to the equilibrium position.

Infrared (IR) intensities can be evaluated through three alternative schemes, as discussed in the previous section.

The CPHF scheme is used for generating, analytically, the Raman intensities. IR and Raman spectra are then easily generated. [124][125][126] Frequencies at k points different from Γ are similarly obtained by making use of supercells and exploiting the translational symmetry of the Hessian matrix to include phonon dispersion effects.

Imaginary frequencies and the SCANMODE keyword

A special case arises when one (or more) of the eigenvalues of the dynamical matrix is negative (imaginary frequency). This happens when the optimized geometry is not a real minimum; the possible reason is a constrained optimization obtained imposing a symmetry higher than the real one, or a unit cell too small. An option is available, SCANMODE, that follows the eigenvector corresponding to the negative eigenvalue, looking for a minimum along this normal coordinate that, if not corresponding to the A 1 symmetry, reduces the point symmetry to the one of the eigenvector. A subsequent optimization usually permits to find the real minimum; the diagonalization of the dynamical matrix produces a full set of positive frequencies.

Isotopic shift and graphical animation of the modes

Identifying the exact nature of modes from observations alone, is extremely difficult, and as a result, unreliable. There are two particular tools, which are implemented in the CRYS-TAL program, that greatly facilitate this task. The first is the ability to change, arbitrarily and selectively, any of the masses within the cell, or fragment thereof, without any additional computational cost, as indicated in equation 5. Within the Born-Oppenheimer approximation, the Hessian matrix, H αi,β j , depends only on the electronic structure of the system, so that a single calculation of this can be used with multiple choices of masses, to give 'isotopically' shifted vibrational frequencies. Since the choice of mass (or masses) is arbitrary, there is no restriction to the mass of a 'real/physical' isotope, so that highly exaggerated changes of the mass can be used solely for the purpose of identifying and characterising a particular vibration. The second is the graphical animation of the modes, whether isotopically shifted or not, which provides a clear description of both the spatial nature and amplitude of the vibrations. A tool for such visualization is available on the CRYSTAL web-site, www.crystal.unito.it. Many other options are available in CRYSTAL related to vibrations, including the calculation of phonon bands and density of states, and of anisotropic displacement parameters (ADP), 123,127 , inelastic neutron scattering spectrum and thermodynamic properties within the quasi-harmonic approximation 128,129 .

The fragment strategy

Suppose we are interested in a system containing ten thousand atoms, that generate three times more vibrational modes, spanning a wavenumber window of, say, 4000 cm -1 of fundamental frequencies including the H stretching. This means that there is a continuous band in many parts of the explored wavenumber range. Such a continuum is essentially useless from the point of view of the characterization of the system (obviously this is not the case for the thermodynamic properties). Only the isolated IR or Raman peaks are of interest in these cases. So, instead of computing the complete Hessian (extremely expensive, but not far from reach with CRYSTAL), one can focus the attention and the computational effort on a fragment of the full system, by building a dynamical matrix involving only a subset of atoms of the full system, and diagonalizing it. This strategy, based on the hypothesis of the local character of these isolated modes is at hand with a simple keyword, FRAG-MENT, and is quite effective. 130 It has been applied recently to the investigation of a set of defects in diamond [131][START_REF] Gentile | Vibrational spectroscopy of hydrogens in diamond: a quantum mechanical treatment[END_REF][133] .

B. The coupled-perturbed Hartree-Fock and Kohn-Sham method

Static electric field

For a molecule in a static electric field (indicated here as E , whereas in Table I and in the rest of the manuscript the symbol ε is used), the permanent dipole moment (µ) and (hyper)polarizabilities (α, β , γ, • • • ) correspond to its energy (E) derivatives of order 1, 2, 3, 4, • • • with respect to the Cartesian components of the field (indicated by the subscripts t, u, v, w in the following equation):

E(E ) = E(0) -∑ t µ t E t - 1 2! ∑ tu α tu E t E u + - 1 3! ∑ tuv β tuv E t E u E v - 1 4! ∑ tuvw γ tuvw E t E u E v E w + . . .
with E(0) the field-free electronic energy. One way of obtaining these derivatives is by a numerical fit to finite field energy values. In the case of periodic systems, however, the finite electric field potential breaks the periodicity. In order to circumvent that difficulty it is necessary to introduce a so-called saw-tooth potential, 134,135 which requires the use of supercells large enough to generate a constant field with a potential that changes sign at the center of each half supercell.

The finite field (FF) method 136 was the first treatment of (hyper)polarizabilities to be implemented in the CRYSTAL code.

An alternative to the FF method is the analytical Coupled-Perturbed-Hartree-Fock/Kohn-Sham (CPHF/KS) scheme [START_REF] Hurst | Ab initio analytic polarizability, first and second hyperpolarizabilities of large conjugated organic molecules: Applications to polyenes C 4 H 6 to C 22 H 24[END_REF] which, in the adaptation for periodic systems, does not have need for supercells. [START_REF] Kirtman | Extension of the Genkin and Mednis treatment for dynamic polarizabilities and hyperpolarizabilities of infinite periodic systems. I. Coupled perturbed Hartree-Fock theory[END_REF] The latter is a self-consistent coupled perturbed approach that takes into account the relaxation of the crystalline orbitals under the effect of an external electric field. In principle, this treatment leads to the same results as the FF method, but is more accurate, more easily extended to (time/frequency)-dependent fields, and is more economical.

The adaptation to periodic systems, while avoiding supercells, is not trivial since the molecular form of the perturbation operator (E • r) is not bound and breaks translational invariance.

Hence, an alternative formulation has to be adopted 62,[START_REF] Otto | Calculation of the polarizability and hyperpolarizabilities of periodic quasi-one-dimensional systems[END_REF]137,138 wherein r is replaced by

Ω k ≡ ıe ık•r ∇ k e -ık•r ≡ r + ı∇ k .
The new perturbation operator, Ω k , is block-diagonal in the same reciprocal k-space as the unperturbed Fock matrix. Hence, the factorization due to translational symmetry is preserved. The CPHF/KS method was first implemented in CRYS-TAL09 for static fields, allowing for the calculation of the electronic contribution to the static polarizability tensor α (or, equivalently, the dielectric tensor ε ∞ ) of closed-and open-shell periodic systems. [START_REF] Ferrero | The calculation of static polarizabilities of 1-3D periodic compounds. The implementation in the CRYSTAL code[END_REF]65 This scheme was further extended to second-order in the perturbed wavefunction in CRYSTAL14, thus allowing for the calculation of static nonlinear properties (namely, first hyperpolarizabilities and second hyperpolarizabilities with the 2n + 1 rule). [START_REF] Ferrero | Calculation of first and second static hyperpolarizabilities of one-to three-dimensional periodic compounds. Implementation in the CRYSTAL code[END_REF]119 The same general approach has also been used for the evaluation of vibrational (hyper)polarizabilities as well as observables obtained as mixed derivatives of the energy with respect to the field and nuclear displacements such as infrared and Raman intensities 66,67,139,140 and the piezoelectric tensor. 141

Electric Field Frequency Dependence

In the two most recent versions of CRYSTAL [START_REF] Dovesi | Quantum-mechanical condensed matter simulations with CRYS-TAL[END_REF][START_REF] Dovesi | CRYSTAL17 User's Manual[END_REF]21,142 the frequency of the field ω has been introduced into the CPHF/KS method. 143,[START_REF] Ferrari | Ab initio calculation of the ultraviolet-visible (UV-vis) absorption spectrum, electron-loss function, and reflectivity of solids[END_REF] The calculated dependence of the electronic linear polarizability on ω (or wavelength λ ) can be directly compared, for example, with experimental data on the photoelastic tensor in the high frequency limit (see reference 118). However, the wavelength of the field must be much larger than the size of the unit cell (λ 1Å) to retain the block-diagonal structure of the Ω k matrix. 145 If the field is frequency-dependent, the determination of the dynamic polarizability for a closed-shell system is obtained as follows:

α uv (-ω; +ω) = -∂ 2 E T OT ∂ E u[-ω] ∂ E v[ω] = -2 n k ℜ P u[-ω],v[+ω] BZ ∑ k Tr n Ω (E u[-ω] ) (k)U (E v[+ω] ) (k) , (8)
where ℜ indicates the real part and Tr the trace, n k is the number of k-points sampling the Brillouin zone (BZ), and P is a permutation operator. 146,147 Ω (E u ) (k) is the matrix representation of the operator u + ı∇ k u in the AO basis; t, u, v are Cartesian directions; and n is the diagonal occupation matrix with eigenvalues equal either to 2 (for occupied orbitals in a closed-shell case) or 0 (for virtual orbitals). The U(k) matrices, that determine the first-order perturbed orbitals, are obtained as: [START_REF] Ferrari | Ab initio calculation of the ultraviolet-visible (UV-vis) absorption spectrum, electron-loss function, and reflectivity of solids[END_REF]148,149 U

(E t[±ω] ) ia (k) = lim η→0 + G (E t[±ω] ) ia (k) ε (0) a (k) -ε (0) i (k) ± ω + ıη , (9)
where G

(E t[±ω]
) ia (k) is given by the derivative of the AO Fock matrix element with respect to the frequency-dependent field along the t-direction, which is subsequently projected onto the unperturbed crystalline orbital basis set,

G (E t [±ω]) ia (k) = ∑ µ,ν C * µi (k)F (E t [±ω]) µ,ν (k)C νa (k) . (10)
In equation ( 9), ε

i(a) (k) is the unperturbed eigenvalue of the occupied i (virtual a) crystalline orbital for each k-point (k) of the reciprocal space. Since F (E t [±ω]) (k) and, thereby, G (E t[±ω] ) (k) depends upon U(k) the solution of equation ( 9) for the U (E t [±ω]) (k) matrices is obtained by fixing the value of ω and then solving iteratively. Knowledge of the U (E t [±ω]) matrices also allows one to study the variation of the high-frequency dielectric matrix ∞ (ω) through the relation:

ε ∞ uv (ω) = δ u,v + 4πα uv (ω) V ( 11
)
where α is the polarizability of a unit cell and V its volume. Far from resonance (ω ε

(0) a (k) -ε (0)
i (k)) the parameter η in equation 9 can be set to zero. On the contrary, near resonance η is determined by the lifetimes of the excited states. Conventionally, the value of this parameter is set equal to the energy resolution of the experimental equipment, which is typically 0.1 eV in UV-visible measurements. The dielectric matrix, then becomes a complex function of frequency and, in terms of the refractive index n = √ ∞ , one can determine the reflectivity R(ω) = | 1-n(ω) 1+n(ω) | 2 as well as the electron loss function ELF(ω)= -ℑ(1/ (ω)).

A formula for the dynamic first hyperpolarizability can be obtained by starting with the 2n + 1 rule expression for the static limit (see equation ( 58) of reference 52) and subsequently taking advantage of the time-dependent HF formulation for molecules developed by Karna and Dupuis. 149 This leads to a general expression for the first hyperpolarizability of closed-shell periodic systems in the presence of frequencydependent fields, which has been implemented in CRYSTAL and validated. 150,151 Equation ( 8) above holds for the HF model. The extension to KS-DFT has been carried out in reference 119 for the static case and in reference 150 for dynamic fields. The dynamic hyperpolarizability is a third rank tensor. It depends on three Cartesian field directions t, u, v and on two incoming frequencies ω 1 and ω 2 (the out-going frequency is ω σ = ω 1 + ω 2 ). Among all possible choices of ω 1 and ω 2 we take as an example the Second Harmonic Generation (SHG). In V β, with V being the unit cell volume (the other non-vanishing component, d zxy , has a very similar dispersion behavior). The curves obtained for the various Hamiltonians remain everywhere parallel, except when approaching the first resonance. The same figure shows the well-known tendency of LDA and GGA functionals to largely overestimate high-order electric susceptibilities, which is exaggerated near the first resonance. When the percentage of HF-exchange decreases, d xyz increases at any frequency. This behavior is not entirely unexpected since it correlates with the predicted band gap, that decreases from HF (14.0 eV) to PBE0 (7.4 eV) and B3LYP (6.9 eV) to PBE (5.2 eV) and LDA (4.8 eV).

The second order perturbation frequency-dependent U(k) matrix, required for the calculation of the dynamic second hyperpolarizability, has not been coded yet in CRYSTAL. However the static finite field (FF) method 136 , combined with CPHF/KS, provides, by best fit of the high-frequency dependent polarizability and first hyperpolarizability, all the non linear electric χ (3) -susceptibilities involved in non linear optics (NLO) processes, such as the electric field induced secondharmonic generation (EFISH) or dc-Kerr effect 151 , with the exception of the third-harmonic generation (THG) and the intensity dependent refractive index (IDRI).

VI. USE OF SYMMETRY IN CRYSTAL

The role of symmetry in the simulation of crystalline compounds is crucial, much more than in the molecular context. This statement, that is obvious when referred to translational symmetry, applies as well to point symmetry. Translational symmetry is the mandatory tool for reducing the infinite system (infinite number of atoms, of basis functions, and then matrices of infinite size) to an infinite set of finite problems (one for each point of the First Brillouin Zone, FBZ, in reciprocal space), whose size is the one of the unit cell (see the first two panels in figure 6). Continuity properties in reciprocal space permit then to sample the FBZ at a relatively small number of points. The exploitation of point symmetry is not mandatory as the translational symmetry is; it is however very useful from many points of view: i) it can improve performance dramatically; ii) it can reduce enormously the required memory storage; FIG. 6. Block-factorization of the Fock matrix for periodic systems. F g : in the basis of AOs, non-packed form (borders are blurry to indicate that such matrix is infinite in principle); F k : in the basis of AO Bloch functions; F k : in the basis of SACO Bloch functions. Also in the last two cases the matrix is infinite (as there is an infinite number of k points in the Brillouin zone), but it is block diagonal.

iii) it helps in comparing simulation with experiment (for example symmetry labels can be assigned to IR and Raman peaks both experimentally and theoretically); iv) last but not least, it greatly simplifies the set of data to be given as input. The data needed to define a crystalline structure is greatly reduced if its space group is known, because the unit cell can be generated automatically from the knowledge of the asymmetric unit only. The same applies to slabs, rods, polymers and molecules. Two significant examples are nanotubes and fullerenes, as obtained by geometrical construction from graphene (only the fractional coordinates of one atom and the length of one lattice parameter must be defined) or, more generally, from a monolayer that is then rolled up. The same kind of simplification applies to many structure manipulations that modify either translational or point symmetry, or both. For instance, a 2D graphene monolayer can be automatically cut from a 3D graphite crystal simply by using the keyword SLABCUT. Obviously i-iv apply when the crystalline system has more symmetry operators than just the identity. Note, however, that: a) the number of crystalline compounds characterized by some symmetry is very high. Whereas large molecules, with few exceptions (say fullerenes), in most of the cases have null symmetry, large unit cell systems with high symmetry are very common. b) even when symmetry is apparently lost, as is the case of defects (translational symmetry broken) or of solid solutions, it turns out that it somehow reappears; in the latter case, for example, symmetry plays a fundamental role in optimizing the sampling of large configurational spaces.

A. Symmetry and efficiency

The key steps where symmetry (translational and point) is used are given in the following for the case where the BFs are constructed from a local basis:

1. Diagonalization of the Fock matrix is restricted to the subset of k points in the irreducible BZ. The eigenval-ues in a star of k points (a set of points that are symmetry related) are the same, and the eigenvectors can be generated by applying symmetry operators.

2. The time required for the calculation of one-and twoelectron integrals is reduced by a factor of up to the number of symmetry operators in the point group. This feature has been implemented in the first release of the code. [START_REF] Dovesi | On the role of symmetry in the ab initio Hartree-Fock linearcombination-of-atomic-orbitals treatment of periodic systems[END_REF] 3. The Fock matrix diagonalization time can be reduced dramatically if performed in the basis of the Symmetry Adapted Crystalline Orbitals (SACO)/Symmetry Adapted Molecular Orbitals (SAMO), shown in the last panel of figure 6. SACO/SAMO are generated automatically in CRYSTAL from the selected basis set of Atomic Orbitals (AO) in the unit cell, with no need for additional information about irreducible representations or characters. This part of the code was implemented about 20 years ago. [START_REF] Zicovich-Wilson | On the use of symmetry-adapted crystalline orbitals in SCF-LCAO periodic calculations. I. The construction of the symmetrized orbitals[END_REF][START_REF] Zicovich-Wilson | On the use of symmetry-adapted crystalline orbitals in SCF-LCAO periodic calculations. II. Implementation of the self-consistent-field scheme and examples[END_REF] The savings factor in computation time is proportional to the second/third power of the ratio between the number of AOs in the basis set (N AO ) and the size of the largest block in the Fock matrix, when represented in the SACO/SAMO basis.

4. Construction of the density matrix scales with the third power of the basis set size (N 3 AO ), as each of the N 2 AO matrix elements is obtained by summing over all occupied crystalline orbitals (a fraction of N AO ). Summations are much shorter if the matrix is constructed in the SACO/SAMO basis first, and then transformed back to the AO basis.

Reduction of computing time thanks to the use of symmetry is not the only issue when handling very large unit cell cases. Memory requirements can also become a bottleneck, if not properly managed at every step of a calculation. Storage of the Fock F, overlap S and density P matrices as full square matrices in the AO-Bloch function basis represents the main bottleneck and needs to be avoided. To make clear the following discussion, we underline that the F, P and S matrices can be represented in the AO basis (as in a molecular context), in the Bloch function basis obtained from the latter by Fourier transform (BF-AO), and in the Bloch function basis obtained from SACOs again by Fourier transform from a direct lattice representation starting from SACOs (BF-SACO).

-As regards the first representation, since the first formulation of the code a vector representation has been used in which, for example, the matrix elements between core orbitals (say Ni 1s ) and other orbitals in the system are very few: only matrix elements are evaluated and stored between this AO and the other AOs on the same atom, but all matrix elements between Ni 1s and its neighbors are disregarded, on the basis of the screening performed with the tolerances defined in the input deck. So these matrices, to be indicated as F g , S g and P g where g labels the (infinite) lattice vectors, are very compact, and increase only linearly with the size of the system. For example, the ratio between the dimension of the Fock matrix F g for chrysotile nanotubes with n=59 and n=29 is 2.02 (see table V), exactly the same as the ratio between the number of atoms. Obviously the number of matrix elements between Ni 2s or Ni 3s and their neighbors increases rapidly. This strategy permits to perform all electron calculations at a cost and memory occupation only moderately larger (say by a factor two to three) than a PP calculation, also for systems containing first row transition metal atoms.

-Since both one-and two-electron integrals are evaluated in the AO basis, the original strategy when SACOs have been implemented [START_REF] Zicovich-Wilson | On the use of symmetry-adapted crystalline orbitals in SCF-LCAO periodic calculations. I. The construction of the symmetrized orbitals[END_REF][START_REF] Zicovich-Wilson | On the use of symmetry-adapted crystalline orbitals in SCF-LCAO periodic calculations. II. Implementation of the self-consistent-field scheme and examples[END_REF] consisted in building first the Fock matrix in the BF-AO basis (large square matrices), and then to transform it to the BF-SACO basis, and back once diagonalized, at each SCF cycle. In this way the advantages of the speed up were not accompanied by the reduction of the memory occupation.

-In the most recent formulation, in 2013-2014 [START_REF] Zicovich-Wilson | On the use of symmetry in SCF calculations. The case of fullerenes and nanotubes[END_REF][START_REF] Orlando | On the full exploitation of symmetry in periodic (as well as molecular) self-consistent-field ab initio calculations[END_REF] F, S and P are directly transformed from the AO to the SACO/SAMO basis (many small matrices instead of a single big one), and then transformed to BF-SACOs (see figure 6, right panel), avoiding then the bottleneck represented by the intermediate step in which large BF-AO matrices appear. For many years the maximum number of symmetry operators in all vector and matrix allocations in CRYSTAL was 48, corresponding to the most symmetric cubic groups. Nevertheless, a more general concept of symmetry was considered more recently, that allowed to extend this limit as exemplified in the following.

-The number of point symmetry operators of icosahedral fullerenes is 120; the group presents irreducible representations (IRREP) with dimensionality from 1 to 5. In the case of largest degeneracy, only one of the five matrices, each one corresponding to one row of the IRREP, must be diagonalized, the eigenvectors of the others being obtained by rotation.

-In section VI C it is shown that for investigating solid solutions supercells must be used in which the number of operators is as large as the product of the number of symmetry operators of the primitive cell of the end members and of the number of translation vectors that are used to generate the supercell adopted for describing the configurations. In the example of section VI C the symmetry operators are 144.

-In section VII the convergence of the thermodynamic properties of MgO is investigated by including the phonon dispersion through the SCELPHONO option. Supercells of increasing size are used; for the largest one, the conventional cell of MgO is expanded by a factor 11 along the three lattice parameters, and contains 10648 atoms, with 95832 AOs. Frequencies are evaluated for these large cells. Symmetry permits to obtain the frequencies just performing in all cases only 3 SCF+G calculations (the equilibrium point, and two points corresponding to a single displacement of Mg and O).

-In section II A 3 the exploitation of the rototranslational symmetry of nanotubes permits to investigate a BN tube as large as (300,0), containing 1200 atoms and characterized by the same number (1200) of symmetry operators.

B. Symmetry in anisotropic properties

A number of physical properties can be given a tensorial representation. For example, linear elasticity is described by a fourth order tensor. Such a tensor consists of 21 independent components for a triclinic crystal, whereas they reduce to 3 in the case of a cubic crystal. Thus, symmetry is key to study tensorial properties: which components must be computed? Which elements are null by symmetry? Which are symmetry-related? Moreover, experimental data are generally reported with some standard orientation of the cell parameters with respect to the Cartesian frame. For comparison between calculated and experimental data to be consistent, the orientation issue needs to be clearly stated. CRYSTAL performs such a symmetry analysis of the tensor using the TENSOR keyword. That is particularly useful in low-symmetric cases and in the case of three-fold rotation axes. Such analysis is automatically performed before the calculation. A single keyword (ELASTCON, for example) is sufficient for generating the full tensor of interest.

The magic symmetry of nanotubes

In the present section an example of the effect of symmetry is provided, with numerical data, and reference to the (n,-n) chrysotile nanotubes, already introduced in section II A 3. Nanotubes are a special case for symmetry, because the number of symmetry operators increases with n, as shown in table V.

It is worth remembering that the input for generating a nanotube from the monolayer is extremely simple, and consists of only two lines: NANOTUBE n1 n2 In the present case, n2 = -n1. This is enough for generating both the symmetry operators and the Cartesian coordinates of all the atoms. More information about formal aspects, examples and inputs/outputs can be found at the CRYSTAL Tutorials web page (see "How to model nanotubes") along with animated graphical illustrations. Let us consider now the data shown in table V, where results are reported for n = 9, 19, 29, 39, 49 and 59. Due to the exploitation of the helical symmetry, 73,[152][153][START_REF] Ferrari | Ab initio periodic study of the conformational behavior of glycine helical homopeptides[END_REF] the computational cost for the full self-consistent field (SCF) and gradient calculation is expected to increase very slowly when passing from the lizardite monolayer (18 atoms and 236 AOs in the unit cell) to the largest tube. Calculations have been performed, for the geometry optimization step of the calculation, with the replicated data version of the code (PCRYSTAL), that is compatible with the use of SACO's, crucial for the drastic reduction of the computational time in diagonalization. The number of atoms, N at , the size of the basis set N AO , the number of symmetry operators (roto-translations and planes), N sym , are reported in columns 2-4. N sym increases from 36 to 236 in such a way that the ratio N at /N sym remains constant. This means that the number of irreducible atoms is the same for all tubes. The size of the irreducible part of the Fock F g matrix, F IRR , reported in the table, is nearly constant and slightly decreasing when the size of the tube increases. For (9,-9) there are about 1753 more matrix elements than for (19,-19), because, due to the small radius R of the tube, there are matrix elements between atoms on the opposite walls of the tube whose value is larger than the thresholds (tolerances of the code). When R increases further, the difference from tube to tube decreases (by 269, 170, 38 and 13 matrix elements). The (very small) difference for the largest tubes is not due to interactions across the tube, but rather to the different curvature of the tubes. The (nearly) constant size of the Fock matrix implies that the number of bielectronic (Coulomb and exchange) and monoelectronic integrals to be computed is constant for all the (n,n) tubes. The N max column gives the maximum size of the matrices to be diagonalized: it is 268 in ALL cases. This means that instead of diagonalizing a matrix whose dimension is 27848×27848, in the (59,-59) case, we diagonalize many matrices whose size is 268×268 or less. We can have a rough estimate (lower limit) of the saving factor from the ratio (27848/268) 2 = 11000, four orders of magnitude! The ratio Ra between F ALL , the number of matrix elements of the complete Fock matrix, and F IRR (the irreducible Fock matrix), that generates F ALL by rotation, is a measure of the real exploitation of symmetry in the various cases. In the ideal case, it should coincide with N sym ; it turns out that actually the two numbers are very close: for example for the largest tube, they are 220 vs 236. In the last column the time for one SCF cycle is reported. The cost of the 2124 atom case is only 7 times larger than the cost of the monolayer, in spite of the fact that the ratio between the atoms in the unit cell is as large as 118. The same ratio between (59,-59) and (9,-9) is only 2.5 for the time and 6.6 for the number of atoms. Partial costs for the various steps of the SCF are given in the three P rows of table VI. We notice that: -the time required by the calculation of the bielectronic integrals, t biel , is essentially constant. It just increases by 10% in going from the smallest to the largest tube.

-The cost of the diagonalization step t diag is a small fraction of the total time per cycle, t cycle . t diag increases from 12% (n=19) to 21% (n=59). As t biel is nearly constant, and t diag is increasing linearly with n, this is not surprising. -t biel + t diag are as much as 67% of t cycle at the lower extreme, but only 47% for (59,-59). Other parts of the SCF are growing rapidly in cost.

-The cost of the monoelectronic part, t mono , remains always a small fraction (10%) of the total. -t P , the reconstruction of the density matrix, that being essentially a matrix multiplication should scale with the third power of the size of the basis, increases by a factor 9 from the two extremes. For the largest tube, it is the most expensive step. -t XC , the time required by the reconstruction of the XC (exchange-correlation) contribution to the Fock matrix, increases by a factor 2, and uses less than 10% of the total time 19, 29, 39, 49 and 59. N at is the number of atoms per cell, N AO is the size of the basis set (number of AOs), N sym the number of symmetry operators, N max is the size of the largest matrix to be diagonalized, F IRR is the size of the irreducible part of the Fock matrix, F ALL is the size of the complete Fock matrix, Ra is the ratio F ALL /F IRR , t is the time for a single SCF cycle; a SCF calculation requires about 25 cycles, and the optimization about 30 steps (small oscillations from case to case), so that the complete optimization requires T=t×25×30. Calculations for the present and next two tables performed on a computer using Intel(R) Xeon(R) CPU E5-2683 v4 @ 2.10GHz cores. Each node contains 2 sockets of 16 physical cores (so 32 physical cores per node) and 128 GB of RAM.

for the largest tube. These data show that the exploitation of symmetry for some steps (reconstruction of the density matrix) can be improved.

-overall, the cost of the calculation increases by only a factor 2.3 when the number of atoms triplicates. Note also that all calculations have been obtained with a very limited amount of CPU time. As a last comment, we notice that the exploitation of symmetry requires: -the classification of atoms, shells, atomic orbitals, construction of the tables that describe symmetry relationships between atoms, -construction of the character table of the group, that in the largest case contains as many as 236 symmetry operators, -the rotation of matrix elements when all the symmetry operators are applied to the irreducible Fock matrix for generating the complete Fock matrix, that will then be Fourier transformed to F(k), the Fock matrix in reciprocal space.

-a similar set of rotations that is applied to the eigenvectors obtained in the SACO basis for generating the density matrix P IRR in direct space, that is then combined with the bielectronic integrals for the next SCF cycle. The cost of all that is not larger than 1-2% of the total.

C. Symmetry and Solid Solutions

An interesting area in which symmetry is very useful are solid solutions. [START_REF] D'arco | Symmetry and random sampling of symmetry independent configurations for the simulation of disordered solids[END_REF][START_REF] Lacivita | On the use of the symmetry-adapted Monte Carlo for an effective sampling of large configuration spaces. The test cases of calcite structured carbonates and melilites[END_REF]155,156 Suppose we are interested in Ca x Mg 1-x CO 3 , and consider for simplicity a particular composition, say 75% Ca and 25% Mg. The space group of the two end members, MgCO 3 and CaCO 3 , is R 3C, with 12 point symmetry operators. A statistically reasonable and computationally affordable choice for the supercell is a 2×2×1 extension of the hexagonal cell, that is 3 times larger than the rhombohedral primitive cell containing two formula units (10 atoms). The supercell contains then 120 atoms (10×3×4), and 144 symmetry operators (12 in the original point group times (3×4) translation operators imported in the supercell). The number of ways of distributing the 6 Ca atoms and the 18 Mg atoms in the 24 lattice positions is 24!/(6!18!) = 134596. Many of these configurations are symmetry related, and then equivalent. The quantum mechanical calculation can be limited to the subset of non symmetry equivalent configurations. Then symmetry is used in CRYSTAL for the selection of the SIC, Symmetry Independent Configurations. It should be noticed that, due to the large numbers of operations (and then computing time) required for the identification and sorting of these equivalences, this is not an easy task. [START_REF] D'arco | Symmetry and random sampling of symmetry independent configurations for the simulation of disordered solids[END_REF][START_REF] Lacivita | On the use of the symmetry-adapted Monte Carlo for an effective sampling of large configuration spaces. The test cases of calcite structured carbonates and melilites[END_REF]155 The number of SIC is, in this case, just 1033, more than 2 orders of magnitude less than the total number of configurations. The SIC are characterized by a residual symmetry (number of symmetry operators that leave invariant the SIC); the number of operators is comprised between 1 (no symmetry) and 144 (maximum symmetry). The number of symmetric SIC (indicated as SSIC) is a small fraction of the SIC: 178 out of 1033; of these 178, the largest part, 143, has multiplicity 72 (only 2 operators leave these configurations invariant). Now, for obtaining at a given temperature the statistical average of a property of the solid solution, one should in principle explore the 1033 SIC; this means to perform a complete quantum mechanical geometry optimization, possibly followed by the calculation of a specific property (for example the vibrational frequencies) for each of them. However it is shown [START_REF] D'arco | Symmetry and random sampling of symmetry independent configurations for the simulation of disordered solids[END_REF][START_REF] Lacivita | On the use of the symmetry-adapted Monte Carlo for an effective sampling of large configuration spaces. The test cases of calcite structured carbonates and melilites[END_REF] that the calculation can be limited to the subset SSIC, that represents a good statistical sample of the full set of SIC, as the two panels in figure 7 suggest. The figure shows that the SSIC span the same interval of the SIC. It also shows that the extremes in the interval are occupied by the SSIC (in particular the low energy side, the most important for the statistical averages). A deeper analysis con-firms that the distribution of the SSIC in each energy interval is proportional to the SIC one, as evident from figure 8, where the distribution of SIC and SSIC per energy interval is shown (histogram) as well as the corresponding integrated quantity (normalized to one). The figure shows that the blue and red curves are very close, and then that the statistical average can be performed including just the SSIC instead of the SIC. It is a common practice in the literature to consider, for cost reasons (as well as, possibly, for the lack of tools able to generate all the configurations in an ordered unique list), a small subset of configurations, obtained, for example, by using some random generator for deciding the occupancy of each lattice site. The above discussion shows that this strategy brings systematic errors, because the low symmetry configurations (the black ones in the figure) are on one hand the ones with the highest probability to be extracted (their number is dominating, 855 vs 178), and, on the other hand, they are the less relevant ones from the point of view of the statistical average as they do not appear in the low energy region in figure 7. Vice versa, the SSIC, the very important ones, have low probability to be extracted. One note should be added: inspection of the low energy region in figure 7 shows that actually a small number of black configurations appears at energies that might have some role for the statistical averages. It can be shown [START_REF] Lacivita | On the use of the symmetry-adapted Monte Carlo for an effective sampling of large configuration spaces. The test cases of calcite structured carbonates and melilites[END_REF] that they are very close to one high symmetry configuration, and can be obtained by simple permutation of two positions. This complement, easily implemented for the black configurations at relatively low energy, further improves the computational scheme, and again requires the application of symmetry. In summary, a very effective strategy is the identification of the SIC and SSIC, and the quantum mechanical calculation of the latter only. For building their list, and for determining their statistical weights, symmetry is crucial. the time required for the calculation of the bielectronic and monoelectronic integrals, of the exchange-correlation contribution to the Fock matrix, for the reconstruction of the density matrix and for the diagonalization of the Fock matrix, respectively. The first two lines refer to the lizardite monolayer, whose unit cell contains 18 atoms. Calculations performed on 64 cores, as described in the caption of table V.

VII. CRYSTAL PERFORMANCE A. Moving the border between quantum mechanical and approximate methods

The frontier separating the systems that can be treated at the quantum mechanical level, or with semi-empirical methods, or classical force fields is moving at high speed, thanks to the rapid progress of hardware and to the (less rapid) improvements of software. Quantum mechanical codes, mostly based on the various flavours of DFT, are becoming more general and more efficient in many respects, although large differences exist among them (one example is the degree and efficiency of parallelization). So one might be tempted to define the largest system that can be investigated at the quantum mechanical level. The question is, however, not sufficiently precise, because many additional variables must be defined besides the size of the system, to make it meaningful. The adopted basis set, the level of the theory (for example LDA, or GGA, or hybrid DFT, or even higher level) as well as the property to be investigated must be specified, because their cost and use of memory can be quite different from case to case. Also the more or less dense nature of the system (the cost per atom of the unit cell of diamond is much higher than the one of a zeolite, when a local basis set is used), the nature of the atoms (a system containing transition metal atoms is 10) 1728 128 65. 4 3.3 5.8 3.5 1.4 51.1 2223(34) 854( 9) 1728 256 36.2 2.4 3.4 3.2 0.9 25.9 1231(34) 457( 9) 1728 512 22.1 1.9 2.2 3.1 0.6 13.8 751(34) 257(9) 2744 128 244.6 5.1 13.0 5.7 3.9 216.2 9784(40) 6575(21) 2744 256 121.7 3.3 7.2 4.4 2.2 104.0 4868(40) 3287(21) 4096 128 778.6 7.8 26.5 8.4 11.2 723.1 31144(40) 19284(21) 4096 256 387.5 4.9 13.9 6.4 5.9 355.6 15500(40) 9466(21) 8000 256 2626.1 9.4 44.7 12.6 20.7 2536.1 89287(34) 57250( 20) 8000 512 1420.0 5.9 25. 4 10.5 11.2 1365.3 48280(34) (48 and 8 symmetry operators, respectively). The number of cycles, that appears in parentheses, is smaller for the displaced point, because in this case the density matrix of the central point is used as an initial guess (GUESSP option in the code). This option permits to reduce the number of cycles by a factor 1.5 to 4. Standard conditions are used for the calculation (T=7, XXLGRID for the exchange correlation integration, a shrinking factor equal to 1 (Γ point only). The SCF convergence threshold is 10 -11 E h , certainly more than necessary for the SCF. Such a large value was used to check the stability of the SCF process. In the case of the largest supercell, if 10 -6 and 10 -8 E h is used, which is certainly sufficient for the evaluation of all ground state properties, the number of cycles reduces to 16 and 24, respectively. VIII. MgO room temperature thermodynamic data computed with supercells of increasing dimension. N at is the number of atoms in the cell, N AO is the corresponding number of atomic orbitals (AO), IS the shrinking factor defining the Monkhorst-Pack grid of points in which the reciprocal space is sampled. S and C v (in J mol. -1 K -1 ) and H 0 and E T (in kJ mol. -1 ) are the entropy, the heat capacity, the enthalpy and the vibrational contribution to the enthalpy. ∆E is the total energy difference per formula unit (in micro-Hartree) with respect to the 8 atoms case (conventional cell), whose energy is -1101.72507579 E h . RMS(ν) (in cm -1 ) is the Root Mean Square of the frequency difference between the 21 modes of the conventional cell and the corresponding modes of the larger supercells. ∆ν max is the maximum difference in the set. Standard computational conditions have been used. more expensive than one containing only first or second row atoms), and the level of required accuracy, that can be controlled by input parameters, can have a dramatic effect on the computational cost. With respect to the many variables listed above, we remind that all the examples of large systems to be presented in the following have been investigated with the global hybrid B3LYP. A global hybrid is considered extremely expensive when used with many periodic codes; this is not the case with CRYSTAL.

N at N AO IS ∆E ∆ν max RMS(ν) S C v H 0 E T 2 18

P and MPP CRYSTAL, and their relationship with symmetry

As anticipated, there are two parallel versions of CRYS-TAL. A replicated data procedure is used by PCRYSTAL, wherein all the most relevant quantities are copied on each node. In MPPCRYSTAL the large matrices are partitioned and distributed amongst the cores; this version of the program, first released in 2010, is advantageous in general for systems with a large unit cell and low symmetry. One of the most important differences is that in PCRYSTAL each matrix to be diagonalized is attributed to a single core, whereas in MPPCRYSTAL it is distributed over as many cores as required. This would suggest to use the latter in any case when the system is large enough. A reasonable strategy consists in attributing no less than 30 -50 basis functions to each core in the diagonalization step. So, for example, with a basis set of 1000 AOs/cell, 20 cores can efficiently diagonalize the Fock matrix. When, however, the point symmetry is non-null, SACOs (see previous section) are very effective in reducing the cost of the diagonalization step. For example, in the case of the already mentioned (10,10) icosahedral fullerene, [START_REF] Zicovich-Wilson | On the use of symmetry in SCF calculations. The case of fullerenes and nanotubes[END_REF]157 the large matrix to be diagonalized (86000 AOs) splits in many small matrices, the largest of which is just 3600x3600 large. This huge reduction suggests obviously the use of SACOs in all the cases in which at least one point symmetry operator exists, apart from the identity. At the moment, however, SACOs are implemented only in the P version of the code, not in MPP. So, in spite of the strong factorization due to SACOs, when the matrices in the SACO basis become larger than, say, 5000, the more economical strategy consists in shifting from P to MPP. The implementation of an MPP strategy, compatible with the use of SACOs, would permit a further step forward in the study of very large systems.

Comparing P and MPP, with and without symmetry: an example

Data in Table VI refer to chrysotile nanotubes (subsection VI B 1), and can help in understanding the interrelation between SACOs, P and MPP. Let us consider the three lines referring to the (39,-39) tube whose unit cell contains 1404 atoms and 18408 AOs. The system is very symmetric (156 operators); so, when running in P mode for, say, geometry optimization, using 64 cores, execution is incredibly fast: 7.2s (CPU time) per SCF cycle (t cycle in the table), only 2.8s for the evaluation of the bielectronic integrals (t biel ) and 1.0s for the diagonalization step, t diag . Let us consider now exactly the same case run with MPP (one of the effects being that the construction of SACOs is switched off, the same effect is obtained by using the NOSACO directive). t cycle increases from 7.2 to 301.9s, mostly due to t diag that jumps from 1.0s to 295.2 s, all the other parts of the calculation remaining essentially unaltered. So the huge advantage of SACOs is in the diagonalization step. Consider finally the same case without symmetry (see the N sym column). In this case t cycle increases to 827.8s, about 115 times more than in the case with symmetry, a factor not far from the number of symmetry operators. Comparing, however, the various contributions to this number, it turns out that the diagonalization time of MPP is the same (plus or minus a couple of seconds) for the case with 156 and 1 symmetry operators (no effect of symmetry in the diagonalization step, as SACOs are not active), whereas the difference is large for all other steps, as symmetry in these parts of the code is active in MPP as is in P.

So for example the time for the calculation of the bielectronic integrals increases by a factor even larger than the number of symmetry operators (from 2.8s to 404.7s). In the last three lines of the table the MPP results for the (59,-59) system without symmetry are reported, when 128, 256 and 512 cores are used. The scaling is quite satisfactory in general, and perfect for t biel and t XC . The diagonalization scaling between 256 and 512 cores is smaller than expected, but in this case the number of AO-BF/core is 54, near the lower border mentioned above for the best use of the MPP diagonalization. Two more comments on P and MPP CRYSTAL: -sequential execution is also possible, and in fact it is a special case of PCRYSTAL, when a single core, instead of n, is attributed to the job.

-a single keyword permits to switch from the P to the MPP version of the code. A multi-task option also exists, that is very useful when many medium size independent tasks must be performed, as it is the case of the many SCF+G calculations required by a frequency calculation of a relatively large system, or when investigating solid solutions, for which many independent configurations must be studied.

Thermodynamic properties of bulk MgO

Among the evidences that CRYSTAL can tackle large systems, we present in tables VII and VIII a set of data referring to MgO supercells containing 2 to 10,648 atoms. The adopted split valence basis set contains 9 AOs for both Mg and O, so that in the largest supercell there are 95832 AOs. For each supercell the vibrational modes are obtained, by using the SCELPHONO option, that uses also the translational symmetry within the supercell in the selection of the irreducible atoms to which the small x, y and z displacements in the numerical calculation of the Hessian matrix are applied. As a result, for all the supercells listed in Table VIII, the Hessian matrix is constructed by computing 3 SCF+G points only, the central one and two points in which one of the two atoms is displaced by 0.003 Å. The point symmetry in the three cases is 48, 8 and 2, respectively. For each supercell the eigenvalues are obtained in an increasing number of k points. For example, in the case of the 8,000 atoms supercell, we obtain the eigenvalues at 4,000 k points. The thermodynamic properties of the crystalline system are obtained by summing up all the eigenvalues of the infinite set (here 4,000) of k points of the BZ.

In the following we list some comments to tables VII and VIII:

1. they provide the evidence that a system with more than 10000 atoms and nearly 100000 AOs can be investigated with CRYSTAL.

2. they provide the evidence that the results are very accurate, as table VIII documents (see section II A).

3. All calculations are performed with MPPCRYSTAL, with use, then, of symmetry in the calculation of the integrals, but not in the diagonalization step (SACOs non-compatible with MPP).

4. In all cases the time per cycle is dominated by the diagonalization step. It costs about 60% for the cells with 1000 atoms, and 96% for the cell with 8000 atoms. This is because the other potentially expensive step, the calculation of the bielectronic integrals, is strongly reduced by the use of symmetry: there are 48 operators in the central point (equilibrium configuration).

5. The same case has been investigated with a variable number of cores. The results for all computed properties are exactly the same. For the construction of the table up to 512 cores have been used. t cycle scales quite well with the number of cores, as a consequence of the good scaling of the diagonalization step.

The scaling of the other steps of the calculation is not easily verified, as the corresponding CPU times are extremely small. It happens to see in the literature calculations referring to very large systems, producing properties whose accuracy can hardly be checked. To eliminate any doubt on the quality of the results obtained in the present case, we add one more evidence to the one provided in section II A about the numerical accuracy of CRYSTAL. Columns 4-6 of table VIII report the total energy (per MgO unit), and two statistical indices obtained from the 21 frequencies that are common to the conventional cell, containing 8 atoms, and to all larger supercells. The two indices are the maximum difference with respect to the wavenumbers of the conventional cell, ∆ν max , and the Root Mean Square RMS(ν) among the sets. It turns out that both ∆ν max and RMS(ν) are always smaller than 0.3 cm -1 , confirming that the accuracy is not only very high, but also constant with the size of the unit cell. We can then consider the thermodynamic properties, resulting from the summation extended to all k points. It turns out that the data for all quantities converge to a well defined value (this means that the numerical error is affecting all these properties on some decimal figures beyond the ones reported in the table). Convergence is rapid: for N at equal to 1000, the difference with respect to the largest calculation is 0.5% for S, 0.1% for C v and even smaller for the other two properties. The graphical representation (not shown) of the four quantities as a function of the number of k points used in the integration is extremely regular, confirming once more the above comments on the negligible effect of the numerical noise.

B. A few examples of investigation on large systems with CRYSTAL

In the following we illustrate a few cases of large systems investigated in the last 7-8 years. 1) The first documentation of the MPP CRYSTAL structure and performance dates back to 2012. In the paper at Ref. 158, whose title is: A New Massively Parallel Version of CRYSTAL for Large Systems on High Performance Computing Architectures, the feasibility of SCF calculations of large supercells containing up to 14 units of MCM-41 (8106 atoms, 108584 AOs) on 512 cores was documented. The mesoporous silica MCM-41 structure consists of a 41x41x12 Å unit cell containing 579 atoms. The supercells were built expanding the MCM-41 unit cell along the c axis. The basis set was a 6-31G * one. The B3LYP hybrid functional was used.

2) In 2016 the structure of the simplest protein, crambin, was investigated 159 in its bulk, hydrated form, exploiting in part neutron scattering information for the position of some of the H atoms, and filling the rest with a set of water molecules. The crambin molecule has 642 atoms; there are 2 molecules in the unit cell. The large voids in the structure were filled with 0, 84 and 172 water molecules, resulting in 1284, 1536, 1800 atoms per cell, and, for the largest case, 16482 AOs. The B3LYP functional was used. At variance with respect to the previous case, here the geometry was fully optimized. This step was a real challenge for the code, due to the presence of the many water molecules, with a relatively large freedom to move, the null symmetry and the very large size of the cell. In all cases the optimization process was successful. Calculations were performed on 1920 cores. It should be underlined that even larger protein systems are at hand, with a sufficiently large number of cores.

3) Between 2012 and 2014 a full family of icosahedral fullerenes was investigated, [START_REF] Zicovich-Wilson | On the use of symmetry in SCF calculations. The case of fullerenes and nanotubes[END_REF]157 up to (10,10), containing 6000 atoms and 86000 AOs. The geometry was optimized. In this case, thanks to symmetry, a single node with 8 cores was enough for the full optimization with the B3LYP functional. This is a strong evidence that symmetry can, in some cases, be as effective as thousands of cores. 4) Recently, the interaction of CO 2 within the pore of the giant-metal organic framework MIL-100 has been investigated. MIL-100 shows an extraordinarily complex structure with about 2800 atoms in the primitive unit cell and contains open metal sites (e.g. Al, Sc, Cr, Fe) that can interact with CO 2 . Dispersion corrected hybrid functionals have been successfully adopted to fully optimise the geometry and predict the binding of CO 2 with the metal sites 160 . In the four previous examples, the code was used for just SCF+G calculations, for geometry optimization, and then for the investigation of many ground state properties. In the next example the analysis extends to the vibrational calculations, that are much more demanding. 5) In 2016 the thermodynamic properties of two garnet systems, pyrope Mg 3 Al 2 Si 3 O 12 and grossular, Ca 3 Al 2 Si 3 O 12 , were investigated 161 computing the phonon dispersion data for six supercells of increasing size, from 80 to 2160 atoms. A size consistency check, similar to the one documented by table VIII for MgO, was performed and was very successful. This was the first example of very large unit cell calculation extended to the vibrational properties. In the same years, large unit cell calculations for defective systems were performed. 162,163 The vibrational spectrum and infrared (IR) intensity of supercells as large as 1000 atoms were computed. 164 Such a large supercell permits to explore geometries in which the lateral interaction among defects is progressively reduced.

C. And then what about scaling?

We think that the previous sections provided sufficient documentation concerning large systems and scaling. Thanks to the efficiency of the code, and to the presence of symmetry in many cases, the need of an extremely large number of cores is rare. A rich documentation can be found in the presentation paper of the 2014 release of CRYSTAL, and in a recent paper. 165 Here we report the response to a challenge concerning the possibility of running codes up to 32.000 cores. In this case two MCM-41 supercells were considered, X16 (9264 atoms/cell) and X24 (13896 atoms/cell and 186144 AOs/cell, resulting from a 6-31G * * basis set). The PBE functional was used, for memory reasons. Calculations were run on the SuperMUC (LRZ, Germany) HPC IBM iData-Plex machine powered by 16 Intel cores per node running at 2.7 GHz, with 2 GB/core. Figure 9 documents: a) that actually the code can run with this large number of atoms and basis functions, on this large number of cores. b) that the scaling (see the numbers on the curves; perfect scaling is the diagonal of the figure, and corresponds to 1) is excellent, larger than 0.8 for X24 up to 8000 cores, and up to 4000 cores for X16. c) that also at 16000 cores the scaling is as large as 0.7 for X24. d) that the larger the system, the better the scaling. e) that also at the limit of 32000 cores, it is not too bad for X24, 0.45.

VIII. PROPERTIES FROM THE WAVEFUNCTION

A very large number of tools has been developed in CRYSTAL for the analysis of the wavefunction, in direct and reciprocal space. The subroutines for the calculation of many of them have recently been parallelized for accessing large systems (see Ref. 166). It is not the aim of the present document to go through the full list; rather, we will give a few examples referring to a specific scientific area, that of defects, that offers the opportunity for illustrating a few of these tools.

A. Defects in diamond

Natural and synthetic diamonds attract significant interest for many possible applications, such as high temperature diodes, transistors, thermistors and detectors, 167,168 thanks to their high melting point, thermal conductivity, hardness and wide band gap. The presence of selected intrinsic and extrinsic defects, which can be incorporated in diamond, [169][170][171][172] can change dramatically these properties. Nitrogen is, by far, the most common impurity, and nitrogencontaining diamonds are classified as type Ib, when a single N atom is substituted for a carbon atom (N s or C defect), or Ia when the N atoms aggregate. The most common aggregates are vicinal N s pairs (A centers) or a combination of four nitrogen atoms surrounding a carbon vacancy V (B centers) 173,174 .

According to the dynamics of the diffusion process, and the geological or preparation history, intermediate situations exist in which there is just one, or two, or three N atoms around the vacancy, before the final stage (the B defect) is reached. The A, B and C centers are then the most common N defects in diamond. Many variants are also possible, with hydrogen atoms saturating part or all the C atoms around the vacancy. 131,175 The A 163 , B 164 and C 176 centers in diamond have been investigated with CRYSTAL, together with many other intrinsic and extrinsic defects.

Spin states

Whereas the ground state of perfect diamond has a closed shell structure (each atom forming strong covalent bonds with the four neighbors), many of the defects present an open shell ground state. Consider for example the carbon vacancy V. Four covalent bonds are broken, and four uncoupled electrons remain unsaturated on the four first neighbors of V. In principle, one could imagine that the defect region reorganizes in such a way to allow the C atoms to couple forming somehow double bonds, as it is the case, for example, when the (111) diamond surface is allowed to relax. The structure of diamond is however too rigid for that, and atoms cannot The ideal situation consists in having a few numbers summarizing the changes with respect to the perfect diamond structure. To this aim, the distances in the defect region (supplemented by some angles, if necessary), the atomic charges and spin momenta, and some indication on the kind of bond linking the atoms, is sufficient to define the chemical situation around the defect. It is clear that from the wavefunction the charge and spin densities can be generated. These are however 3D functions, whose representation (2D maps or 3D plots) might not satisfy the need of a concise description of the local situation. In CRYSTAL various ways of integrating the charge (and spin) distribution functions are available: -The topological analysis of the electron charge density, performed according to Bader's prescriptions using Gatti's TO-POND package, 178 that is integrated with the CRYSTAL program; -The Hirshfeld population analysis.

The reader can obtain more information and references for the above two tools in the CRYSTAL manual 5 ; for the latter, see also Ref. 179 -The Mulliken population analysis. It is the simplest, easy to read and quickly computed tool for the analysis of the charge and spin densities. It provides net charges, bond populations, for each type of atom surrounding the vacancy for the VN 0 2 , VN - 2 and VN 0 1 defects in diamond. Data have been obtained for the doublet (d) and quadruplet (q) spin states for VN - 2 and VN 0 1 , and for the triplet (t) and singlet (s) for VN 0 2 . Energy differences between high and low spin states ∆E (in eV) are reported. In all cases the low spin state is the ground state. Multiplicity (number of equivalent atoms) is given in parentheses. Note that in the VN 0 1 case there are two sets of carbon atoms, that are shown in two contiguous lines. Distances R (in Å) and bond populations b refer to the N atom(s) and its (their) first neighbors, and to the carbon atoms around the vacancy, and to their first neighbors. When more than two values exist, ranges are given in square brackets. atomic magnetic moments. The main criticism to the Mulliken analysis concerns the way the charge (spin) density is partitioned among atoms (halfand-half). So the results of this analysis can change when the basis set is changing. However, when comparing similar situations treated with the same basis set, trends and differences are meaningful, and this analysis turns out to be very useful. In figure 10 the distances, the net charges and the bond populations around the VN 0s 2 defect (the neutral VN 2 defect in its ground singlet state) are reported. In table IX the data of the excited triplet state, as well as of a couple of other defects in diamond, namely the negatively charged variant of VN 0 2 , indicated as VN - 2 , and VN 0 1 are reported for comparison. For the negatively charged defect, the neutrality of the unit cell, required by the infinite nature of the system, is obtained by adding a uniform charge background whose integral is equal to a unitary positive charge. A few comments: 1) in perfect diamond, each carbon atom is obviously neutral (net charge Q C =0 |e|). So the net charges shown in the figure and table measure the charge displaced as a consequence of the presence of the vacancy V and of the two nitrogen atoms.

System Q C Q N µ C µ N R CC R NC b CC b NC ∆E VN 0 2 s +0.
2) The net charge on nitrogen is about the same in the six cases, varying from -0.38 to -0.41 |e|. Nitrogen is accumulating nearly half an electron.

3) The net charge on the first three carbon neighbors of N is +0.13 |e|; this means that the NC 3 cluster is essentially neutral. This is not surprising, as the electronic structure of diamond is extremely rigid, and it is very difficult (energetically expensive) to pull electrons from or add them to a C-C bond. 4) Also the net charges of the C V atoms (the unsaturated carbon atoms around the vacancy) have limited variations from case to case: from +0.06 to +0.12 |e|. 5) The net charges of second and third neighbors of the vacancy fall rapidly to zero, and the perturbation turns out to be very localized. 6) The integral of the spin density extended to the unit cell is equal to 0. (singlet), 1. (doublet), 2. (triplet), 3. (quadruplet). 7) A large fraction of this spin density is concentrated on the two (VN 0 2 and VN - 2 ) or three (VN 0 1 ) C V atoms. For example, it is 0.88 |e| in the quadruplet state of VN 0 1 ; so 0.88x3 = 2.64 |e| out of 3 electrons are concentrated on the three C V . In the case of the charged defect VN - 2 the spin density is more diffuse, and the sum must be extended to more neighbors for approaching the formal spin momentum of the unit cell. 8) The magnetic moment µ on C V is about the same for the three high spin cases: 0.89. 0.88 and 0.88 |e| for VN 0 2 , VN - 2 and VN 0 1 , respectively. 9) In the low spin cases µ decreases to 0.81 |e| (VN 0 2 ) and 0.68, 0.80, 0.80 |e| (VN 0 1 ). In the charged defect, in which the Pauli repulsion is higher in the vacancy region, µ on C V is as small as 0.47 |e|.

10) The magnetic moment on the nitrogen atoms is always very small: close to zero in the low spin cases, and in between 0.05 and 0.15 |e| in the high spin cases. 11) In spite of the high concentration of the spin density on the atoms around the defect, the residual spin density is far from short ranged. In general (see figure 11) the µ value on the n-th neighbor has opposite sign with respect to the one on the (n+1)-th neighbor of the vacancy, according to a sort of spin density wave. As a consequence also atoms relatively far from the defect region can contribute to the EPR signal.

12) The equilibrium geometry, apart from the very local situation, is not much sensitive to the kind of defect, or to its charge state. The C-C distance in pristine diamond is 1.56 Å. It decreases to 1.50-1.52 Å (see table IX) for the C-C and N-C bonds involving one atom on the border of the vacancy. The next C-C bonds, however, are in between 1.55 and 1.56 Å, confirming that the geometrical perturbation is very localized. 13) The same holds for the b bond populations between the defect atoms and their neighbors, shown in the table: C-C and N-C are large and positive, indicating strong covalent bonds; for the next bonds, the 0.34 |e| value of pristine diamond is recovered. In summary, the simple Mulliken population analysis is able to provide a rich, although not complete, description of the electronic structure around the defect. A similar description is much more cumbersome if a PW basis set is used. Further information on the nature of the unpaired electron, his spatial extent, and the kind of atomic orbitals involved is provided by the spin density maps, like the one of figure 12 which shows the accumulation of spin density on the carbon atoms around the vacancy, and how it propagates to the next neighbors. The figure might be slightly misleading, as the high spin density peaks are truncated at ±0.1 |e| /(a 3 0 ) (hence the large white area on the two sides of these atoms). The shape of the isodensity lines is the one of a p orbital. The two lobes are asymmetric, with the largest one pointing towards the empty space of the vacancy, where however there is competition (or Pauli repulsion) with the other C V atom. Short range repul-sion and spin polarization involve obviously also the N atom, on the other side of the vacancy, and the carbon atoms linked to C V . The comparison of the two figures helps to understand why the singlet state is more stable than the triplet. The presence of an additional nodal plane between the two C V atoms increases the kinetic energy, and destabilizes the singlet state. The different spin on the two C V atoms (α in one case, β in the other), strongly reduces, however, the Pauli repulsion, which is dominant at these short distances. The spin pressure between the two C V atoms, that forces to move part (a very small fraction) of spin density in the direction opposite to the vacancy, and then generates a spin polarization on the atoms linked to C v , is smaller in the singlet case, as the number of isodensity lines (6 for s and more than 10 for t, remember that the density is truncated at the tenth isoline) documents. Note also that the polarization on N essentially disappears (see also the µ values in table IX).

The EPR parameters

Two important parameters can be evaluated with CRYSTAL, that are related to the spin density. The Fermi contact term for an unpaired electron

A n iso = 8π 3 g e µ e g n µ N |ψ(R n )| 2 (12)
arises from the direct interaction of nuclear and electron spins, and is only non-zero for states with finite electron spin density, |ψ(R n )| 2 , at the nuclear site R n (see figure 11). In Eq.12 the parameters g e , µ e , g n and µ N are the free-electron g-factor, Bohr magneton, gyromagnetic ratio of I n and nuclear magneton respectively.

The second is the traceless tensor B n at nucleus R n , defined as:

B n i j = g e µ e g n µ N dr n |ψ(r n )| 2 3r ni r n j -|r n | 2 δ i j |r n | 5 (13) where r n is the electron distance to R n . It samples the longrange properties of the wavefunction as a result of the r -3 scaling. It is essential, therefore, for basis sets to be sufficiently flexible so that these differing dependencies can be adequately described; a quite rich basis set (17 AO/atom), tailored in particular for EPR calculations, has been used for generating the data discussed in this section. 177 The spin density maps, figure 12, are a 2D representation of the 3D spin density function that, when integrated according to equation 13, provides the B n i j values. When the square of the wavefunction magnitude is partitioned according to the Mulliken scheme, the µ values of table IX are obtained. Table X compares the calculated [START_REF] Ferrari | Substitutional nitrogen in diamond: A quantum mechanical investigation of the electronic and spectroscopic properties[END_REF]180 and experimental A iso and B tensor values for two defects, namely the C (substitutional nitrogen) and the VN - 2 center. In the first case there is un unpaired electron on one of the first neighbors of N, hence a doublet state. In the second case there are 3 unpaired electrons, and in principle both a quadruplet and doublet state are Atom A iso B 1 B 2 B 3 reference 82.0 23. 2 -11.6 -11.6 Calc. [START_REF] Ferrari | Substitutional nitrogen in diamond: A quantum mechanical investigation of the electronic and spectroscopic properties[END_REF] 13 C 1 208. 208.9 137.9 -68.9 -68.9 Calc. [START_REF] Ferrari | Substitutional nitrogen in diamond: A quantum mechanical investigation of the electronic and spectroscopic properties[END_REF] 13 C 2 35.5 7.0 -3.5 -3.5 Exp. 188 34. 9 6.4 -3.4 -3.0 Calc. [START_REF] Ferrari | Substitutional nitrogen in diamond: A quantum mechanical investigation of the electronic and spectroscopic properties[END_REF] 15 N +4.0 -0.6 +0.5 +0.1 Exp. 181 +4.1 -0.6 +0.5 +0.1 Calc. possible. The latter is however by far the most stable, and to it refer the experimental data by Green et al.. 181 For the C defect other experimental determinations exist, 182-184 that differ by less than 1% from the ones reported in the table. The agreement between simulation and experiment is excellent. A similar good agreement is observed for other cases, for example the N + 2 center (see references 185 and 186), the positively charged version of the A center (two vicinal substitutional N atoms). The defect is ionized in order to have unpaired electrons, and then to generate an EPR signal. This confirms that, at least for light atoms, the present approach produces excellent results.

Other related properties can be obtained from the simulation (and sometimes from the experiment).

-When the local symmetry is low enough, the symmetric tensor B is not diagonal, and then usually six parameters are produced, the three eigenvalues, and the three eigenvectors, or the three angles defining the orientation of the local frame with respect to the absolute frame. These three angles are obviously a by-product of the calculation.

-The Electric Field Gradient tensor at the nuclear position, that couples with the quadrupolar nuclear momentum when the nuclear spin I is equal or larger than 1 (in atomic units). In many cases it has been computed and compared with experiment since long (1998). For example in reference 189 where the so called trapped holes in MgO, CaO and SrO due to Li and Na substitutions for oxygen are investigated at the Hartree-Fock level, that performs very well in localizing the hole, at variance with other functionals containing a smaller fraction of exact exchange, or no exchange at all. Two very recent applications are in Ref. 186 and 180. The same ingredients can be used for estimating the two important parameters of the Mössbauer effect, namely the chemical shift and the quadrupolar splitting) (see the CRYSTAL manual and Ref. 190).

-At the moment CRYSTAL is on the contrary unable to com-pute the g tensor. A (not completely satisfactory) solution to this limitation consists in obtaining it from a cluster treated with a molecular code, once verified that the periodic and cluster values of A iso , B and P coincide, or are very similar. The amount of data can be much larger if more EPR active atoms around the defect center are considered. From the simulation side this is obviously not a problem (see the A n iso and B values reported in Ref. [START_REF] Ferrari | Substitutional nitrogen in diamond: A quantum mechanical investigation of the electronic and spectroscopic properties[END_REF] for many C atoms at increasing distances from the N defect). The accuracy and resolution of the experimental spectrum can on the contrary matter. In the case of the neutral Li vacancy in LiF (color center), very accurate data produced more that 50 years ago 191 provide partial EPR information up to the seventh neighbors of the vacancy. In a simulation by Mallia et al. 192 , a quite reasonable agreement is observed with the experiment, that might improve if the experimental labelling of some of the farther neighbors is changed.

IX. CONCLUSIONS

To the authors' knowledge, CRYSTAL was the first periodic ab initio code to be distributed publicly (1988). Despite the rapidly evolving computational technology of quantum molecular codes, and even more so of periodic codes, all the major ingredients of the current program were already in place at that time. However, there have been many additions, generalizations and technical improvements. Among the special features that may not be available elsewhere we cite: 1) Extensive use of symmetry to dramatically reduce CPU requirements.

2) Global hybrids, such as B3LYP (available from the nineties) and PBE0, can be used at a cost that is only 3 to 5 times larger than PBE and GGA. This also applies to Hartree Fock. A large set of other functionals spanning Jacob's ladder is also available.

3) The ability to treat large systems. For example, the wavenumbers and IR intensities of systems with a unit cell containing much more than 1000 atoms can readily be obtained. In table VII such calculations for a MgO supercell containing 10648 atoms are documented. The general treatment for systems of that size (total energy, geometry optimization, vibrational frequencies, and the full set of properties) is not a far-off target. On the other hand, ab initio calculations for vibrational properties of systems with, say, more than 200-300 atoms using other codes are extremely rare. 4) Scaling with the number of cores is excellent -up to thousands of cores for sufficiently large systems. 5) Careful attention has been paid to achieving high accuracy. For any new implementation the level of accuracy has always been discussed and documented. Referring again to vibrational frequencies, we mention in this regard the accurate reproduction of experimental spectra in the terahertz region 193 . 6) A very broad array of tensor property calculations is available (elastic, piezolectric, photoelastic, dielectric, static and dynamic first and second-order hyperpolarizability, hyperfine coupling, electric field gradient, IR and Raman intensities, etc.). 7) There are also many tools for characterization of the wavefunction/density in direct or reciprocal space (Mulliken, Bader and Hirshfeld analysis, Wannier localization, x-ray structure factors, electron momentum distributions and Compton profiles). In addition, graphical tools have been developed during the years to help visualization and plotting of computed properties [194][195][196][197][198] . Even though the capabilities of CRYSTAL have been greatly enhanced since its first release in 1988 there are still a number of limitations, some of which we hope to remove in the near future. Amongst the latter are the availability of magnetic and mixed electric-magnetic properties, as well as spinorbit coupling, and a general treatment of anharmonicity. The efficiency of the code must also be improved, for example by parallelizing those parts of the CPHF routines that deal with Raman intensities. At the present stage CRYSTAL represents a reliable, general, accurate computational tool. Its stability over more than 40 years is a strong warranty for any user.

X. ACKNOWLEDGMENTS

A number of people, who have contributed to the formulation of the algorithms and their implementation in the CRYSTAL code over more than 40 years (1976-2020), are no longer with us. Cesare Pisani (he died in 2011) was the first to propose this project. He, together with Carla Roetti (she died in 2010), were two of the three original co-authors (the other being Roberto Dovesi). Long-term collaborators include Victor Ronald Saunders (from 1980 to 2005, when he suffered a transient ischemic attack that stopped his, still very brilliant, scientific activity; he passed away in 2018), Roberto Orlando (1986-2016), Claudio Marcelo Zicovich-Wilson (1994-2016). Finally, we mention Alessio Meyer, who was a second-year PhD student when he died suddenly in 2009. CRYSTAL would simply not have been born, nor would it have achieved its accuracy, generality, and efficiency, without their contributions.

RD would like to thank William Mackrodt and Andreas Savin for accurate reading of some parts of this manuscript, correction of errors, and suggestions for improving it. The data that supports the findings of this study are available within the article and references therein.

FIG. 3 .

 3 FIG. 3. A single-walled chrysotile nanotube (top); two single layers of lizardite with different orientations (center); and a layer of antigorite (bottom). The orientations of the lizardite layers are chosen to match the orientation assumed in the nanotube and in antigorite. Si is coloured in light brown; Mg in green; H in white; O atoms that are bound to Si and coordinated to Mg are coloured in blue; the remaining O atoms are coloured in red.

FIG. 4 .

 4 FIG.4. Variation of the two Raman active modes (A 1 symmetry) as a function of n for a set of (n, -n) chrysotile nanotubes with n ranging from 9 to 59. The continuous curve has been obtained by best fit with the function ν = a+b/n+c/n 2 +d/n 3 . The a coefficient represents then the limiting value for n → ∞. Dotted lines indicate the corresponding lizardite values (3800.9 and 702.1 cm -1 ).

  ) In this case z = 6 (each TM atom is surrounded by 6 TM atoms with opposite spin); and S 2 = 1 (two electrons with spin 1/2) for Ni, and S 2 = 25/4 for Mn (five electrons, S=5/2). The J values obtained with equation 3 are shown in the penultimate column of TableIV B. It turns out that: -The HF values are about 1/3 of the experimental values (see the caption to the table).

FIG. 5 .

 5 FIG.5. The FM and AFM structure of KMnF 3 . In the left panel, all the Mn atoms (at the vertices of the cube) have the same spin (open circles). In the central panel the most stable AFM configuration is shown, with alternating (111) sheets of α and β spin Mn atoms; each Mn atom is surrounded by 6 Mn atoms with opposite spin. In the AFM' structure (to the right), only the four Mn neighbors in the basal plane have opposite spin. The comparison of the energies of AFM and AFM' permits to check the additivity of the superexchange interaction.

  [START_REF] Pisani | Hartree-Fock Ab Initio Treatment of Crystalline Systems[END_REF] expectation values for the KMnF 3 and KNiF 3 . J and Jy are the magnetic coupling constants obtained with equation 3 and 4, respectively. Experimental values span from 3.65 to 3.70 K 114 for KMnF 3 and from 44 to 51 K 115 for KNiF 3 .

  Figure 1 of Ref. 151 the d xyz component of crystalline urea is shown as a function of the wavelength for several Hamiltonians. The d tensor is obtained directly from the hyperpolarizability β via the relation d = π

FIG. 7 .

 7 FIG. 7. Energy distribution of the 1033 SIC (Symmetry Independent Configurations, bottom) and of the 178 SSIC, characterized by a multiplicity lower than 144 (top). The multiplicity is indicated by the colors in the bottom panel, to the right. A configuration with multiplicity, say, 4, is invariant to the action of 144/4 = 36 of the 144 symmetry operators. The zero energy is defined by the most stable configuration.

FIG. 9 .

 9 FIG.9. Scaling with the number of cores as a function of the system size on the SuperMUC system. X16 and X24 are obtained by multiplying the MCM-41 cell by 16 and 24, respectively, along c.

  FIG. 10. The VN 0 2 defect. In panel a) the pristine diamond conventional cell is shown, whereas in panel b) a local cluster of C and N atoms adjacent to the vacancy and including 12 next-nearest neighbors is represented. Mulliken net charges (in |e| units; boldface) are shown near the atoms, for each symmetry irreducible atom. Bond distances (in Å, italics) and Mulliken bond populations (in |e|) are also shown in between the two involved atoms. Data refer to the singlet ground state.

FIG. 11 .

 11 FIG.11. Spin density profiles of the VN 0 2 defect in its singlet (top) and triplet (bottom) spin states along a path connecting the atoms with relatively high spin density. Some of the atomic labels appear in panel b) of figure10, others refer to atoms connected to atoms at the border of the reported cluster. All calculations were done with the S 216 supercell, B3LYP functional and the 6-31G-J* basis set.177

  FIG.11. Spin density profiles of the VN 0 2 defect in its singlet (top) and triplet (bottom) spin states along a path connecting the atoms with relatively high spin density. Some of the atomic labels appear in panel b) of figure10, others refer to atoms connected to atoms at the border of the reported cluster. All calculations were done with the S 216 supercell, B3LYP functional and the 6-31G-J* basis set.177

2 .

 2 The charge and spin density, and the wavefunction analysisOnce the geometry optimized, an overview of what happened upon the insertion of the defect is provided by the local geometry and the charge distribution.

FIG. 12 .

 12 FIG. 12. Spin density maps of the VN 0 2 defect in its singlet (top) and triplet (bottom) spin states. Two carbon (in grey) and one of the nitrogen atoms (in blue) surrounding the vacancy define the plane of each figure. Isodensity lines differ by 0.01 |e| /(a 3 0 ). Spin density is truncated at ±0.1 |e| /(a 3 0 ). Continuous, dot-dashed and dashed lines represent positive, null and negative spin densities, respectively.

  Access to the HPC resources of CINES/IDRIS/TGCC obtained thanks to the grants 2018-[A0040810471] (Fabien Pascale) and 2018-[A0050810537] (Philippe D'Arco) made by GENCI are warmly acknowledged. High Performance Computing resources were partially provided by the EXPLOR centre hosted by the Université de Lorraine.

  Energy E and atomic position vector r a , of atom a;Energy E, electric field and atomic position vector r a ; A

	Property	Tensor rank	Formula	Definitions	Method
	atomic gradient	1	g a t = ∂ E ∂ r a,t		A
	cell gradient polarizability Born charge (IR intensity)	2 2 2	G tu = ∂ E ∂ A tu α tu = -∂ 2 E ∂ ε t ∂ ε u Z a * tu = -∂ 2 E ∂ ε t ∂ r a,u	Energy E and lattice metric matrix A = (a 1 , a 2 , a 3 ); Energy E and electric field ;	A A
	electric field gradient Hessian	2 2	q tu = -∂ ε t ∂ r u H ab tu = ∂ 2 E ∂ r a,t ∂ r b,u	Electric field and position vector r; Energy E and atomic position vector r a ;	A A/N
	direct piezoelectricity	3	e tuv = ∂ P t ∂ η uv	Polarization P and rank-2 strain tensor η;	A/N
	converse piezoelectricity	3	d tuv = ∂ η uv ∂ ε t	Electric field and rank-2 strain tensor η;	A/N
	first hyperpolarizability	3

TABLE II .

 II Properties of perfect bulk MgO evaluated with supercells of different shape, size and symmetry. T summarizes the value of the five thresholds controlling the truncation of the Coulomb and exchange series, where T= T 1 =T 2 =T 3 =T 4 =1/2T 5 . N is the number of atoms per cell.

				.0566 368.8 111.2 189.0 1.71320 6.7442 436.9 4254.9 24 189
	7 40 C(1,1,5) 8 16 B3LYP -0.0087 368.8 111.2 189.0 1.71321 6.7443 437.4 4247.7 24 40
	7 80 C(1,1,10) 8 16 B3LYP -0.0371 368.8 111.2 189.0 1.71321 6.7438 437.2 4247.0 25 75
	7 120 C(1,1,15) 8 16 B3LYP -0.0207 368.8 111.2 189.0 1.71321 6.7411 437.3 4247.7 26 110
	7 2	P	12 48 HF -0.2358 389.5 103.1 200.4 1.51445 5.9081 465.1 4157.7 18	1
	7 8 C(1,1,1) 8 48 HF	0.0000 389.5 103.1 200.4 1.51445 5.9155 464.7 4157.7 20	2
	7 64 C(2,2,2) 8 48 HF -0.0004 389.5 103.2 200.4 1.51445 5.9136 464.8 4157.5 22	2

TABLE III .

 III Wavenumbers (first column, in cm -1 ) and intensities (km mol. -1 ) of the 17 IR active modes of pyrope. The latter were computed with three different algorithms, namely the Berry Phase (BP), Wannier Function (WF) and Coupled Perturbed Hartree-Fock/Kohn Sham (CP) schemes. The absolute and percentage differences (δ and δ %) among the three cases are shown in the last three columns. For example, BP-CP stands for the absolute value of 100*(BP-CP)/BP. I is the mean intensity, |∆| and |∆ max | are the mean and max difference for the full set of vibrations. The shrinking factor for reciprocal space sampling is IS=3; the threshold for truncation of the bielectronic integrals is T = 7; the functional is B3LYP. For more details see reference 53.

		Modes	Intensities			δ (δ %)
		BP	WF	CP	(BP-CP) (CP-WF) (WF-BP)
	1	113.66 3101 3098 3100 1 (0.10)	2 (0.06)	3 (0.03)
	2	137.04 32	32	31	1 (0.00)	1 (3.13)	0 (3.23)
	3	184.89 3303 3299 3298 5 (0.12)	1 (0.03)	4 (0.15)
	4	212.97 0	0	0	0 (0.00)	0 (0.00)	0 (0.00)
	5	258.64 646 646 649	3 (0.00)	3 (0.46)	0 (0.46)
	6	332.64 6384 6372 6368 16 (0.19) 4 (0.06) 12 (0.25)
	7	345.82 45	45	43	2 (0.00)	2 (4.44)	0 (4.65)
	8	379.57 3457 3449 3458 1 (0.23)	9 (0.26)	8 (0.03)
	9	420.81 1367 1359 1370 3 (0.59) 11 (0.81) 8 (0.22)
	10	456.55 13689 13681 13664 25 (0.06) 17 (0.12) 8 (0.18)
	11	480.33 852 847 842 10 (0.59) 5 (0.59)	5 (1.19)
	12	530.44 794 795 798	4 (0.13)	3 (0.38)	1 (0.50)
	13	580.18 1412 1407 1402 10 (0.36) 5 (0.36)	5 (0.71)
	14	672.73 3	3	3	0 (0.00)	0 (0.00)	0 (0.00)
	15	864.83 14100 14056 14043 57 (0.31) 13 (0.09) 44 (0.41)
	16	895.95 5668 5670 5673 5 (0.04)	3 (0.05)	2 (0.09)
	17	970.25 5802 5788 5780 22 (0.24) 8 (0.14) 14 (0.38)
	I	3568 3562 3560 9.71 (0.22) 5.12 (0.05) 6.71 (0.17)
	|∆|%	-			0.17	0.65	0.73
	|∆ max | %	-			0.59	4.44	4.65

  29406 30.0039 4.9996 1.24 3.71 PBE0 1.34328 30.0042 4.9742 5.66 16.95 B3LYP 1.87056 30.0049 4.8799 7.88 23.51 PBE 3.34401 30.0039 4.9639 14.08 42.17 KNiF 3 HF 0.56853 6.0057 2.0030 14.96 44.85 PBE0 2.61530 6.0029 1.9622 68.82 204.38 B3LYP 3.20714 6.0030 1.9453 84.39 249.58 PBE 7.90360 6.0019 1.6996 207.98 580.10

TABLE IV .

 IV

HF, PBE0, B3LYP and PBE energy difference per magnetic center (in mE h ), high and low spin S

  N at N AO N sym N max F IRR F

				ALL Ra t
	Slab	18 236	6	10736 50036 5 1.6
	(9,-9)	324 4248 36 268 29272 987606 34 4.5
	(19,-19) 684 8968 76 268 27519 1954872 71 5.0
	(29,-29) 1044 13688 116 268 27250 2952722 108 5.8
	(39,-39) 1404 18408 156 268 27080 3958734 146 7.2
	(49,-49) 1764 23128 196 268 27042 4947726 183 9.6
	(59,-59) 2124 27848 236 268 27029 5956876 220 11.4
	TABLE V. Computational data for lizardite (slab) and chrysotile nan-
	otubes (n,-n) with n=9,

TABLE VI .

 VI FIG. 8. Distribution of the 1033 SIC (red histogram) and 178 SSIC (blue histogram) per energy interval. Note that each SIC and SSIC is multiplied by its multiplicity, to give the total number of 134596 configurations. The continuous curves represent the integrated quantities, normalized to one (see the scale to the right). System N core N sym t cycle t biel t mono t XC t P t diag Cost in CPU seconds of the various steps of the SCF calculation of three (n,-n) chrysotile nanotubes, when SACOs are used (with P CRYSTAL, P rows) or not (with MPP CRYSTAL, MPP rows). For obtaining the set of frequencies, 61 SCF+G calculations are performed, 60 of which correspond to configurations without symmetry (one atom is displaced). N core is the number of utilized cores. t cycle is the times required by an SCF cycle. t biel , t mono , t XC , t P , t diag are

	Slab 1T P	64	6	1.4 0.4 0.6 0.1 0.0 0.2
	P	64	1	4.2 1.6 1.6 0.5 0.0 0.2
	(19,-19) P	64	76 5.0 2.7 0.6 0.5 0.4 0.6
	(39,-39) P	64 156 7.2 2.8 0.9 0.8 1.6 1.0
	(59,-59) P	64 236 11.4 3.0 1.3 0.9 3.6 2.4
	(19,-19) MPP 64	76 43.1 2.7 0.6 0.6 0.5 38.8
	MPP 64	1 280.7 190.5 24.7 25.6 0.7 38.7
	(39,-39) MPP 64 156 301.9 2.5 1.2 1.0 2.0 295.2
	MPP 128 156 160.4 1.2 1.3 0.9 1.5 153.0
	MPP 64	1 827.8 404.7 61.5 59.5 3.6 296.7
	MPP 128	1 445.9 201.1 30.5 30.4 2.5 150.8
	(59,-59) MPP 128 236 502.1 1.9 0.9 2.2 4.4 490.1
	MPP 256 236 291.4 1.2 0.5 2.1 3.0 284.3
	MPP 128	1 957.5 323.6 54.4 48.9 4.4 488.6
	MPP 256	1 510.5 159.8 33.5 25.3 3.0 286.7
	MPP 512	1 355.9 87.6 15.8 15.3 3.9 229.8

  N at N core t cycle t biel t mono t XC t P t diag

			tt central	tt displ
	1000 32	43.5 5.1 6.6 4.9 1.9 24.5 1609(37)	760(8)
	1000 64	29.6 3.3 3.8 3.1 0.9 18.1	947(32)	552(10)
	1000 128 16.9 2.0 2.6 2.0 0.5 9.0	574(34)	292(10)
	1000 256 11.5 1.5 1.6 2.4 0.3 5.5	391(34)	171(10)
	1000 512	9.2 1.2 1.0 2.4 0.2 4.2	313(34)	113(

TABLE VII .

 VII Relevant CPU times for a frequency calculation of MgO with large supercells. All calculations performed with the MPP (see text) version of the code. N core is the number of cores utilized. t cycle is the time required by an SCF cycle at the equilibrium geometry, characterized by 48 symmetry operators. t biel , t mono , t XC , t P and t diag as in Table VI. t cycle = t biel +t mono +t XC +t P +t diag . tt central and tt displ are the total time for the full SCF of the central and the first displaced point

	30894(20)
	10648 512 6238.0 7.6 43.1 11.7 21.9 6151.4 212092(34) 130580(20)

TABLE IX .

 IX Atomic Mulliken net charges Q, spin densities µ and bond populations b (in |e|)

TABLE X .

 X Experimental and calculated hyperfine coupling constants (MHz) for the P1 (first three rows) and VN - 2 (last two rows) defective center. B tensor components are sorted so that |B 1 | > |B 2 | > |B 3 |.

R. A. Evarestov, Quantum Chemistry of Solids. LCAO Treatment of Crys-