CRYSTAL09 is a major release and the most relevant new features are:
Static polarizability and dielectric tensor through a Coupled Perturbed HF/KS scheme
The static polarizability of crystalline systems, the dielectric
tensor, the refractive indices, the birifrangence and optical axes
(all these quantities can be obtained from the first one)
are computed through a new implementation of the CPHF(KS)
scheme. A tutorial
(High frequency dielectric constants by CPHF(KS))
is available on how to use this new feature.
Alternatively, a saw-tooth finite-field scheme can be applied
to the periodic system; A tutorial
(High frequency dielectric constants by FIELD) is also available for this option.
Limits and merits of this approach with respect to
CPHF(KS) are illustrated.
M. Ferrero, M. Rérat, R. Orlando, R. Dovesi
Coupled perturbed Hartree-Fock for periodic systems:
the role of symmetry and related computational aspects
J. Chem. Phys. 128, 014110 (2008)
M. Ferrero, M. Rérat, P. Orlando, R. Dovesi
The calculation of static polarizabilities of periodic compounds.
The implementation in the CRYSTAL code for 1D, 2D and 3D systems
J. Comput. Chem. 29, 1450-1459 (2008)
M. Ferrero, M. Rérat, B. Kirtman, R. Dovesi
Calculation of first and second static hyperpolarizabilities
of one- to three-dimensional periodic compounds. Implementation
in the CRYSTAL code
J. Chem. Phys. 129, 244110 (2008)
Phonon dispersion using a direct approach
and infrared intensities through a Berry phase approach
Phonon dispersion using a direct approach
The calculation of phonon frequencies at k points different form Gamma is now possible.
This enables a more accurate calculation of the thermodynamic properties and
the comparison of calculated frequencies with experimental Inelastic Neutron Scattering (INS) spectra.
The scheme is based on a direct approach which implies the calculation of the
Hessian matrix on supercells. An efficient exploitation of symmetry helps the reduction
of computational cost, by finding the irreducible set of atomic displacements compatible
with the supercell.
A new section has been added to the “Vibrational Frequencies Calculation” tutorial
to explain the new keywords.
Infrared intensities through a Berry phase approach
A new scheme to compute IR intensities has heen implemented
based on a Berry phase approach.
This scheme parallels the localized Wannier-based one
to which it offers a cost effective alternative.
Transition state search
Several tools that allow molecules, polymers, slabs and crystals to be optimized in valence coordinates
as well as a suitable saddle point optimization technique to search for transition state structures
for this kind of systems have been implemented.
The adoption of these localized coordinate systems largely facilitates the study of chemical processes
in periodic systems with atomic connectivity, as occurs in catalytic reactions on zeolites, clathrates
or oxidic surfaces.
The new features have been illustrated to study the proton jump between oxygen atoms
of the Brønsted site in the H-chabazite zeolite.
A tutorial is also available on how to use these new features
A. Rimola, C. M. Zicovich-Wilson, R. Dovesi, P. Ugliengo
The distinguished reaction coordinate method: search and characterization of
transition state structures in crystalline systems using the periodic CRYSTAL program
J. Chem. Theory Comput. 6 (2010) 1341
Constant pressure geometry optimization of cell constants and atomic positions
The calculation of the stress tensor and related properties have been implemented
in the CRYSTAL code. The stress tensor is obtained from analytical gradients with respect
to the cell parameters. The pressure and the enthalpy are then computed.
The possibility of applying external pressure has also been implemented.
The constant-pressure optimization offers an alternative optimization method
in addition to the already implemented optimization at constant volume.
Analytical stress tensor and pressure calculations with the CRYSTAL code
Mol. Phys. (2010) (DOI: 10.1080/00268970903193028)
Automated calculation of the elastic tensor of crystalline systems
An automated procedure for calculating second-order elastic constants for
crystalline systems of any symmetry has been implemented.
Second derivatives with respect to strain are evaluated numerically
from analytical gradients. The internal co-ordinates are re-optimized with
each applied strain. Point group symmetry is exploited to reduce the number
of needed deformations according to Laue’s classes.
W.F. Perger, J. Cryswell, B. Civalleri and R. Dovesi
Ab initio calculation of elastic constants of cristalline systems with the CRYSTAL code
Comput. Phys. Commun. 180, 1753-1759 (2009)
Automated E vs V calculation for equation of state
E vs V curves can now be computed automatically by just specifying
a range of volumes. A constant volume geometry optimization is perfomed
for each point and then results are fitted to most common equation
of states and polynomial.
Work is in progress to extend this scheme to E vs P calculations.
New DFT methods
New GGA functionals for solids
The set of available GGA functionals in CRYSTAL has been extended
to include the newly developed GGA functionals for solids, namely:
the PBEsol one  for exchange and correlation, the exchange Wu-Cohen
functional  and the exchange SOGGA of Zhao-Truhlar . Also,
the Wilson-Levy correlation functional  has been added.
Old and new functionals have recently been tested for different solids [5-7]
 J. P. Perdew, A. Ruzsinsky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou, K.
Burke, Phys. Rev. Lett. 100 (2008) 136406
 Z. Wu, R. Cohen, Phys. Rev. B 73 (2006) 235116
 Y. Zhao, D. G. Truhlar J. Chem. Phys. 128 (2008) 184109
 L. C. Wilson, M. Levy, Phys. Rev. B 41 (1990) 12930
 R. Demichelis, B. Civalleri, M. Ferrabone, R. Dovesi, Int. J. Quantum Chem. 110 (2010) 406
 D. I. Bilc, R. Orlando, G. M. Rignanese, J. Iniguez, P. Ghosez, Phys. Rev. B 77 (2008) 165107
 B. Civalleri, D. Middlemiss, R. Orlando, C. Wilson, P. Ugliengo, Chem. Phys. Lett. 451 (2008) 287
London-type empirical correction for dispersion interactions
An empirical correction term to include dispersion interactions in DFT methods
as proposed by Stefan Grimme  is available.
It is a damped pairwise London-type term of the form
The correction is applied to the computed ab initio total energy and gradients.
Therefore, total energy calculation, geometry optimization and vibrational frequency
calculation can be carried out by including the empirical dispersion correction.
The current implementation has been mainly tested and used in combination
with the B3LYP method with application to molecular crystals and layered
 S. Grimme, J. Comput. Chem., 2006, 27, 1787
 B. Civalleri, C.M. Zicovich-Wilson, L. Valenzano, P. Ugliengo, CrystEngComm, 2008, 10, 405; 1693(E)
 P. Ugliengo, C.M. Zicovich-Wilson, S. Tosoni, B. Civalleri, J. Mater. Chem., 2009, 19, 2564
New tools to build 1D structures
Three different options are now
available for the creation of 1D structures, namely:
- POLYMER cases whose max order of the rotational axis is 6
- HELIX (all cases with helical symmetry) with order up to 48
- NANOTUBE (nanotubular systems having helical symmetry)
the input MUST start from a 2D layer, or 3D layered system
Symmetry is fully exploited, so that calculations involving
hundreds of atoms using good-quality basis sets and hybrid HF/DFT functionals
are feasible at a relatively low cost.
The present limit for the roto-translation group is now 48
operators; this limit is going to disappear in the next few months.
Automatic generation of nanotubes from single-layer systems
New options were implemented in the CRYSTAL code, which
permit in a relatively easy way to generate 1D nanotubes
starting from 2D structures, or by cutting a slab from a 3D bulk
and then roll it up into a tube.
This strategy allows users to exploit the full roto-translational
symmetry of nanotubes (helical groups) and to drastically reduce
the computational costs .
(--> Nanotube systems)
shows how to prepare a basic input, describes the related options, how to re-use
the optimized geometry for a given (N1,N2) tube as a
starting point for a new (N1',N2') tube and so on.
Relevant references on recent applications (carbon
nanotubes, inorganic tubes such as imogolite and chrysotile [1,2]),
examples, graphical demonstrations and exercices are also provided.
 Y. Noel, Ph. D'Arco, R. Demichelis, C.M. Zicovich-Wilson and R. Dovesi
On the use of symmetry in the ab initio quantum mechanical simulation of
nanotubes and related materials
J. Comput. Chem. 31 (2010) 855
 Ph. D'Arco, Y. Noel, R. Demichelis, R. Dovesi
Single-layered chrysotile nanotubes: A quantum mechanical ab initio simulation
J. Chem. Phys. 131 (2009) 204701
Helical symmetry for polymers
Helical polymers can now be modelled by fully exploiting the roto-translational
symmetry of the helix (up to order 48).
This has been recently applied to investigate the conformational behavior
of polyglycine helical infinite chains.
A. M. Ferrari, B. Civalleri, R. Dovesi
Ab initio periodic study of the conformational behavior of glycine helical homopeptides
J. Comput. Chem. 31 (2010) 1777 [DOI:10.1002/jcc.21468]
New tools for initial guess of SCF for d- and f-partly
A new option has been inserted that permits to define the occupation of specific f
in a given shell for open shell systems. Electrons belonging to partially filled shells
can then be assigned to selected AOs.
This allows user to have a better control of the initial guess
and an easier convergence for the SCF process.
New tools treatment of solid solutions
A new option permits to investigate solid solutions
through the following steps:
select an appropriate supercell of the system: the
dimension must be adequate to describe the kind of
disorder you are interested in;
define the number of Y atoms to be substituted: the
symmetry independent configurations
are listed automatically by the code;
- optimize the structure of a few configurations
(quantum mechanical calculations for each configuration);
- define a TWO-BODY model (TBM) involving neighbours up to
a given distance;
- define by best-fit the parameters of the TBM;
- compute with the TBM the energy of all the symmetry
independent configurations, and order them by energy;
- restart from point (3) adding new configurations to the
list of the fully optimized ones;
- when the list in (7) proposes already explored
configurations as the most stable, and
when the parameters of the model are stable with respect to the
number of data in the fitting, perform the statistical
mechanics final calculations.
with more information and examples is available.
A. Meyer, P. D'Arco, R. Orlando, R. Dovesi, J. Phys. Chem. C 113 (2009) 14507
A. Meyer, P. D'Arco, R. Orlando, C. M. Zicovich-Wilson, L. Maschio, R. Dovesi,
J. Chem. Phys. submitted
Revised implementation of Electron Momentum Density analysis and Compton profiles
The CRYSTAL program allows the accurate determination of the one-electron Density Matrix (DM)
of crystalline systems, both at Hartree-Fock level and within the Density Functional Theory;
Since the earlier public versions of the program, many DM-related quantities,
both in direct and momenta space, can be easily evaluated.
The majority of the algorithms which concern these properties has been improved in the CRYSTAL09
version of the program and new algorithms have been implemented; in particular:
- The calculation of the Electron Momentum Density (EMD),
evaluated along a line or within a plane, is generalized to f-type orbitals and to the open-shell case;
- A new algorithm has been implemented for the computation of the EMD along any crystallographic direction
starting from the DM instead of the eigenvectors;
- A new algorithm has been implemented for the computation of the directional Compton Profiles (CP)
as 1D Fourier Transform of the autocorrelation function and not as the integral of the EMD.
With such a scheme CPs can be computed also with f-type orbitals and for open-shell systems.
A tutorial (EMD)
with more information and examples is available.
A. Erba, C. Pisani, S. Casassa, L. Maschio, M. Schutz, D. Usvyat,
Phys. Rev. B 81 (2010) 165108
Enhanced Massive-parallel version (MPPcrystal - distributed memory)
An enhanced Massive parallel version of the code is available that allows
user to reach an improved scaling in parallel execution on super-computing
resources for very large systems.
Here is an example of the scaling of the code with the number of processors
up to 1024 cores.
Note: MPPcrystal does not implement: (i) IR intensities calculation and (ii)
coupled perturbed CPHF/CPKS scheme
P. Ugliengo, M. Sodupe, F. Musso, I.J. Bush, R. Orlando, R. Dovesi
Realistic Models of Hydroxylated Amorphous Silica Surfaces and MCM-41
Mesoporous Material Simulated by Large-scale Periodic B3LYP Calculations
Adv. Mater. 20 (2008) 4579