In this example, we will optimize the geometry of Cu_{2}O using PBE0 hybrid density functional method and TZVP basis set. In other words, the level of theory is DFT-PBE0/TZVP.

**Table of contents**

# Geometry block

Here is the geometry that was converted into CRYSTAL input format from Cu2O_Pn-3m.cif:

**Cu2O_Pn-3m_PBE0_TZVP_opt.d12**

Cu2O (Pn-3m) geometry optimization (DFT-PBE0/TZVP) CRYSTAL 0 0 0 224 4.267 2 29 0.00 0.00 0.00 8 0.25 0.25 0.25 ENDGEOM

Note that the `TESTGEOM`

and `FINDSYM`

keywords must be removed from the input when you want to run a real calculation.

We can tell CRYSTAL that we want to optimize the geometry by including `OPTGEOM `

block in the geometry block:

**Cu2O_Pn-3m_PBE0_TZVP_opt.d12**

Cu2O (Pn-3m) geometry optimization (DFT-PBE0/TZVP) CRYSTAL 0 0 0 224 4.267 2 29 0.00 0.00 0.00 8 0.25 0.25 0.25 OPTGEOM ENDOPT ENDGEOM

The default settings for the `OPTGEOM`

block are good for most purposes. See Chapter 6 of CRYSTAL manual for full description of the keywords that can be used in the OPTGEOM block.

The OPTGEOM block must be closed with keyword `END`

. It makes sense to use `ENDOPT`

to differentiate from the `ENDGEOM`

keyword that will close the whole geometry block.

# Basis set block

The choice of the quantum chemical method and the Gaussian-type basis set determine the quality and cost of your calculation. The choice of the basis set is a vast topic on its own and here we will simply use basis sets taken from the literature. Hopefully in the near future this guide will also provide a library of basis sets across the periodic table.

In the present case, we will use a TZVP basis sets taken from J. Linnera, A. J. Karttunen, *Phys. Rev. B* **2017**, *96*, 014304 (DOI).

We need to give a basis set for both O and Cu. The basis sets are listed below in the CRYSTAL input format. The O basis set starts from the line 1, where `8 7`

means that the basis set is for *Z* = 8 and it has 7 "shells" (with type *s*, *sp*, *p*, *d*, or *f*). Cu basis set starts from the line 25, where `29 12`

means *Z* = 29 and 12 shells.

**Cu2O_Pn-3m_PBE0_TZVP_opt.d12**

8 7 0 0 6 2.0 1.0 27032.382631 0.21726302465E-03 4052.3871392 0.16838662199E-02 922.32722710 0.87395616265E-02 261.24070989 0.35239968808E-01 85.354641351 0.11153519115 31.035035245 0.25588953961 0 0 2 2.0 1.0 12.271113873 0.39768730901 4.9159842006 0.24627849430 0 0 1 0.0 1.0 0.90086482370 1.0000000000 0 1 1 0.0 1.0 0.25000000000 1.0 1.0 0 2 4 4.0 1.0 75.300554155 0.60685103418E-02 17.743733858 0.41912575824E-01 5.5355828651 0.16153841088 2.0685535103 0.35706951311 0 2 1 0.0 1.0 0.78238772422 1.0000000000 0 3 1 0.0 1.0 1.2000000000 1.0000000000 29 12 0 0 8 2.0 1.0 377518.79923 0.22811766128E-03 56589.984311 0.17688035931E-02 12878.711706 0.91993460227E-02 3645.3752143 0.37411016434E-01 1187.0072945 0.12189873737 426.46421902 0.28983900714 165.70660164 0.41531872174 65.598942707 0.21905799287 0 0 4 2.0 1.0 414.41265811 -0.24682525053E-01 128.32056039 -0.11716827406 20.622089750 0.55301315941 8.7821226045 0.52242718609 0 0 2 2.0 1.0 13.741372006 -0.22736061821 2.2431246833 0.71761210873 0 0 1 1.0 1.0 0.89370549079 1.0000000000 0 0 1 0.0 1.0 0.35 1.0000000000 0 1 1 0.0 1.0 0.14 1.0 1.0 0 2 6 6.0 1.0 2034.7596692 0.23524822298E-02 481.90468106 0.19134070751E-01 154.67482963 0.90171105278E-01 57.740576969 0.26063284735 23.099052811 0.42093485770 9.3882478591 0.24344615121 0 2 3 6.0 1.0 37.596171210 -0.28991094530E-01 5.1240690810 0.54919083831 2.0119996085 0.93793330488 0 2 1 0.0 1.0 0.73860686002 1.0000000000 0 3 4 10.0 1.0 74.129460637 0.14363216676E-01 21.359842587 0.86628177096E-01 7.4995240537 0.25631430541 2.7601394169 0.40374062368 0 3 1 0.0 1.0 0.95362061236 0.39427042447 0 3 1 0.0 1.0 0.28420862520 0.23091146816 99 0 ENDBAS

Discussing the basis set input in detail is out of the scope of this guide. Basically, the basis set input lists the exponents and (contraction) coefficients for Gaussian functions that describe atomic orbitals. Analogously to molecular orbitals in molecular calculations, crystalline orbitals are formed as linear combinations of the atomic orbitals.

The basis set block ends with line **99 0** and keyword `ENDBAS`

.

# Computational parameters block

The third and final block of the CRYSTAL input determines various computational parameters. For our DFT-PBE0 calculation, the block looks as follows:

**Cu2O_Pn-3m_PBE0_TZVP_opt.d12**

DFT PBE0 ENDDFT SHRINK 0 8 8 8 8 TOLINTEG 8 8 8 8 16 MAXCYCLE 100 FMIXING 80 EXCHSIZE 30000000 BIPOSIZE 30000000 END

Line 1 opens the `DFT`

sub-block, where we tell CRYSTAL that we want to use the PBE0 hybrid density functional. The block is closed with `ENDDFT`

.

Line 4, `SHRINK`

, determines the *k*-mesh used for sampling the reciprocal space according to the Monkhorst-Pack scheme. This is a big sub-topic on its own, but for now, the guidelines are as follows:

- The denser the mesh, the more accurate the sampling of the reciprocal space.
- A denser mesh also means computationally more expensive calculation. We need to balance accuracy and performance.
- The example sets an 8 × 8 × 8
*k*-mesh for the reciprocal space (line 6). - A good starting point is to use
`SHRINK`

, where`SHRINK`

* (real space primitive cell edge) is at least 30 Å. For example, here 8 * 4.267 Å = 34.136 Å. This is on the safe side. - One should always benchmark the
*k*-mesh to see that the targeted feature, such as total energy, is converged with respect to the*k*-mesh.

Line 7, `TOLINTEG`

refers to truncation criteria for bielectronic integrals. This is a crucial, but rather complex setting and for the vast majority of systems, it is enough to use `8 8 8 8 16`

. The default is `6 6 6 6 12`

, but this is too loose.

Line 9, `MAXCYCLE`

, simply tells the maximum number of cycles in the Self-Consistent Field procedure (SCF). Typically 100 is enough for almost every system and you rarely need to increase this. If you have SCF convergence problems, it is better to solve it in some other way than increasing the number of cycles.

Line 11, `FMIXING`

, is a percentage that refers to Fock mixing during SCF (Fock matrix from the previous cycle mixed with the new Fock matrix). This improves the convergence behavior of the SCF. 80% is a high value, but by default this is only used for two SCF cycles before DIIS convergence accelerator kicks in. In the case of **SCF convergence problems**, it is typically a good start to rely longer on FMIXING and turn DIIS on only when the energy change has become relatively small. You can achieve this by adding keyword `THREDIIS `

`0.01`

`EXCHSIZE`

and `BIPOSIZE`

are just memory buffer sizes available for CRYSTAL. Default is ten times smaller than the values in the example. Increasing these keywords to the values here reduces disk traffic, which is generally a good idea because disk is slow. The cost is a larger use of memory. If you run out of memory, these can be reduced.

# Complete input file

Here is the final input file:

**Cu2O_Pn-3m_PBE0_TZVP_opt.d12**

Cu2O (Pn-3m) geometry optimization (DFT-PBE0/TZVP) CRYSTAL 0 0 0 224 4.267 2 29 0.00 0.00 0.00 8 0.25 0.25 0.25 OPTGEOM ENDOPT ENDGEOM 8 7 0 0 6 2.0 1.0 27032.382631 0.21726302465E-03 4052.3871392 0.16838662199E-02 922.32722710 0.87395616265E-02 261.24070989 0.35239968808E-01 85.354641351 0.11153519115 31.035035245 0.25588953961 0 0 2 2.0 1.0 12.271113873 0.39768730901 4.9159842006 0.24627849430 0 0 1 0.0 1.0 0.90086482370 1.0000000000 0 1 1 0.0 1.0 0.25000000000 1.0 1.0 0 2 4 4.0 1.0 75.300554155 0.60685103418E-02 17.743733858 0.41912575824E-01 5.5355828651 0.16153841088 2.0685535103 0.35706951311 0 2 1 0.0 1.0 0.78238772422 1.0000000000 0 3 1 0.0 1.0 1.2000000000 1.0000000000 29 12 0 0 8 2.0 1.0 377518.79923 0.22811766128E-03 56589.984311 0.17688035931E-02 12878.711706 0.91993460227E-02 3645.3752143 0.37411016434E-01 1187.0072945 0.12189873737 426.46421902 0.28983900714 165.70660164 0.41531872174 65.598942707 0.21905799287 0 0 4 2.0 1.0 414.41265811 -0.24682525053E-01 128.32056039 -0.11716827406 20.622089750 0.55301315941 8.7821226045 0.52242718609 0 0 2 2.0 1.0 13.741372006 -0.22736061821 2.2431246833 0.71761210873 0 0 1 1.0 1.0 0.89370549079 1.0000000000 0 0 1 0.0 1.0 0.35 1.0000000000 0 1 1 0.0 1.0 0.14 1.0 1.0 0 2 6 6.0 1.0 2034.7596692 0.23524822298E-02 481.90468106 0.19134070751E-01 154.67482963 0.90171105278E-01 57.740576969 0.26063284735 23.099052811 0.42093485770 9.3882478591 0.24344615121 0 2 3 6.0 1.0 37.596171210 -0.28991094530E-01 5.1240690810 0.54919083831 2.0119996085 0.93793330488 0 2 1 0.0 1.0 0.73860686002 1.0000000000 0 3 4 10.0 1.0 74.129460637 0.14363216676E-01 21.359842587 0.86628177096E-01 7.4995240537 0.25631430541 2.7601394169 0.40374062368 0 3 1 0.0 1.0 0.95362061236 0.39427042447 0 3 1 0.0 1.0 0.28420862520 0.23091146816 99 0 ENDBAS DFT PBE0 ENDDFT SHRINK 0 8 8 8 8 TOLINTEG 8 8 8 8 16 MAXCYCLE 100 FMIXING 80 EXCHSIZE 30000000 BIPOSIZE 30000000 END

# Submitting the job

Prerequisites for submitting the job to queue:

- You need to have a terminal connection to the server.
- You need to be in the same directory as the input file.

The job can now be submitted using `jsub`

:

jsub -np 8 crystal Cu2O_Pn-3m_PBE0_TZVP_opt.d12

This will run the job using 8 processors. See the general guidelines for submitting CRYSTAL jobs for more information.

# Short review of the output

Here is an example of the output file from the geometry optimization of Cu_{2}O: Cu2O_Pn-3m_PBE0_TZVP_opt.o.

## Geometry and basis set

1) First, CRYSTAL prints out the same information on the geometry as for a TESTGEOM job (together with a summary of the geometry optimization settings under `COMBINED CELL/ATOM OPTIMIZATION CONTROL`

).

2) Next, CRYSTAL prints out the details of the GTO basis set used in the calculation under the banner

*******************************************************************************

LOCAL ATOMIC FUNCTIONS BASIS SET

*******************************************************************************

This information is sometimes needed for example for Density of States analysis.

3) Under the basis set, CRYSTAL also prints `NUMBER OF AO`

(AO = atomic orbitals = number of basis functions, here 168) and N. OF ELECTRONS PER CELL (here, 132). These are very important values if you ever have problems for example in calculating reaction energies. The number of electrons for the reactants and products must be the same.

4) Next, CRYSTAL prints `NUMBER OF K POINTS IN THE IBZ`

(Irreducible Brillouin Zone) and lists `K POINTS COORDINATES`

. Even though we are using an 8 × 8 × 8 *k*-mesh for sampling the reciprocal space, we only have 35 *k*-points in the final mesh. The reduction of k-points is due to the symmetry of Cu_{2}O (48 symmetry operations). The *k*-point coordinates are sometimes needed for example for the analysis of electronic band structures or crystalline orbitals.

5) Finally, CRYSTAL prints a very important table, NEIGHBORS OF THE NON-EQUIVALENT ATOMS. From the table we can extract for example the following information:

N = NUMBER OF NEIGHBORS AT DISTANCE R

ATOM N R/ANG R/AU NEIGHBORS (ATOM LABELS AND CELL INDICES)

1 CU 2 1.8477 3.4916 5 O 0 0 0 6 O 0 0 0

5 O 4 1.8477 3.4916 1 CU 0 0 0 2 CU 0 0 0 3 CU 0 0 0

4 CU 0 0 0

Cu atom 1 has two neighbors at a distance of 1.8477 Å. Those neighbors are oxygen atoms. Oxygen atom 5 has four neighbors at a distance of 1.8477 Å and those neighbors are Cu atoms. These distances are for the initial geometry. A similar table is printed in the end of the file for the final geometry.

If you encounter SCF convergence problems in your calculation, it is always good to check the neighbor analysis first. Are the distances as you would expect or are there some abnormally short distances? Short interactomic distances are one of the most typical reasons for SCF convergence issues.

## SCF convergence

CRYSTAL starts the Self-consistent field (SCF) procedure from atomic wavefunctions. For the first two SCF cycles, it uses simple Fock mixing (FMIXING) to control the SCF procedure (`CYCLE 1 - MIX 80 %`

). After that, DIIS convergence accelerator kicks in (`DIIS ACTIVE - HISTORY: 2 CYCLES`

). SCF converges in 8 cycles and the final printout is:

CYC 8 ETOT(AU) -6.711744025915E+03 DETOT 5.95E-08 tst 1.49E-10 PX 9.65E-06

TTTTTTTTTTTTTTTTTTTTTTTTTTTTTT DIIS TELAPSE 39.18 TCPU 39.05

DIIS TEST: 0.14893E-10 AT SCF CYCLE 8 - DIIS ACTIVE - HISTORY: 8 CYCLES

TTTTTTTTTTTTTTTTTTTTTTTTTTTTTT FDIK TELAPSE 39.21 TCPU 39.07

INSULATING STATE

TOP OF VALENCE BANDS - BAND 66; K 1; EIG -1.8331020E-01 AU

BOTTOM OF VIRTUAL BANDS - BAND 67; K 1; EIG -9.5306379E-02 AU

DIRECT ENERGY BAND GAP: 2.3947 eV

TTTTTTTTTTTTTTTTTTTTTTTTTTTTTT PDIG TELAPSE 39.25 TCPU 39.11

CHARGE NORMALIZATION FACTOR 1.00000000

TOTAL ATOMIC CHARGES:

28.5306620 28.5306620 28.5306620 28.5306620 8.9386759 8.9386759

The first line shows the current total energy (`ETOT`

) and the change of energy (`DETOT`

). The default SCF convergence criterion during geometry optimization is 1E-07.

CRYSTAL also prints out that the system is in `INSULATING STATE`

, meaning that it has a band gap. The (direct) band gap is 2.3947 eV and also the information on the highest-energy valence band and lowest-energy virtual (conduction) band is printed. Note that the band gap printed here has been obtained for the *k*-mesh used in the calculation. An electronic structure calculation with denser *k*-mesh needs to be carried out to get a final value for the band gap.

The data under `TOTAL ATOMIC CHARGES`

should also be checked, to see that the atomic charges are reasonable. These are based on a Mulliken analysis. In general, Mulliken charges behave better in solid-state calculations compared to gas-phase calculations. Here, Cu atoms are formally Cu(I) and their partial charge is 29 − 28.53 = +0.47 e^{−}. For O atoms, the partial charge is −0.94 e^{−}. These are reasonable values.

## Forces

The optimization next proceeds to FORCE CALCULATION and reports the forces acting on each atom. For Cu_{2}O, these are practically zero (1E-14). The atoms are constrained to their positions by space group symmetry (`THERE ARE NO SYMMETRY ALLOWED DIRECTIONS`

) and only the cell parameter *a* can change.

## Optimization and final optimized geometry

The geometry optimization the proceeds and the geometry is printed for every step until the convergence is reached:

MAX GRADIENT 0.000082 THRESHOLD 0.000450 CONVERGED YES

RMS GRADIENT 0.000082 THRESHOLD 0.000300 CONVERGED YES

MAX DISPLAC. 0.000014 THRESHOLD 0.001800 CONVERGED YES

RMS DISPLAC. 0.000014 THRESHOLD 0.001200 CONVERGED YES

CRYSTAL then prints the final optimized energy, band gap, neighbor analysis, and final geometry

******************************************************************

* OPT END - CONVERGED * E(AU): -6.711745667537E+03 POINTS 8 *

******************************************************************

BAND GAP: 2.3105 eV

[...]

N = NUMBER OF NEIGHBORS AT DISTANCE R

ATOM N R/ANG R/AU NEIGHBORS (ATOM LABELS AND CELL INDICES)

1 CU 2 1.8696 3.5330 5 O 0 0 0 6 O 0 0 0

[...]

FINAL OPTIMIZED GEOMETRY - DIMENSIONALITY OF THE SYSTEM 3

(NON PERIODIC DIRECTION: LATTICE PARAMETER FORMALLY SET TO 500)

*******************************************************************************

LATTICE PARAMETERS (ANGSTROMS AND DEGREES) - BOHR = 0.5291772083 ANGSTROM

PRIMITIVE CELL - CENTRING CODE 1/0 VOLUME= 80.489901 - DENSITY 5.853 g/cm^3

A B C ALPHA BETA GAMMA

4.31764698 4.31764698 4.31764698 90.000000 90.000000 90.000000

*******************************************************************************

ATOMS IN THE ASYMMETRIC UNIT 2 - ATOMS IN THE UNIT CELL: 6

ATOM X/A Y/B Z/C

*******************************************************************************

1 T 29 CU 0.000000000000E+00 0.000000000000E+00 0.000000000000E+00

2 F 29 CU 5.000000000000E-01 -5.000000000000E-01 0.000000000000E+00

3 F 29 CU 0.000000000000E+00 -5.000000000000E-01 -5.000000000000E-01

4 F 29 CU 5.000000000000E-01 0.000000000000E+00 -5.000000000000E-01

5 T 8 O 2.500000000000E-01 2.500000000000E-01 2.500000000000E-01

6 F 8 O -2.500000000000E-01 -2.500000000000E-01 -2.500000000000E-01

The time used for the calculation is reported on the line

TTTTTTTTTTTTTTTTTTTTTTTTTTTTTT END TELAPSE 462.47 TCPU 461.94

TELAPSE is wallclock time in seconds. If it differs a lot from TCPU, the code is probably using a lot of time in writing to disk or communicating between the CPUs.

# Create a CIF of the optimized structure

After the geometry optimization, the CRYSTAL output could be directly viewed in Jmol or CRYSPLOT (Make a plot → Geometry Structure). However, often it is useful to convert the optimized crystal structure to CIF format for visualization in other programs such as VESTA.

1) Create a new input based on the geometry optimization input:

cp Cu2O_Pn-3m_PBE0_TZVP_opt.d12 Cu2O_Pn-3m_PBE0_TZVP_opt_done.d12

2) Use *crygeom* script to get optimized geometry:

crygeom Cu2O_Pn-3m_PBE0_TZVP_opt.o

The command will output the following:

4.31764698 4.31764698 4.31764698 90.000000 90.000000 90.000000

2

29 0.000000000000E+00 0.000000000000E+00 0.000000000000E+00

8 2.500000000000E-01 2.500000000000E-01 2.500000000000E-01

3) Replace the geometry section in the new input *Cu2O_Pn-3m_PBE0_TZVP_opt_done.d12* with these four lines from *crygeom*.

4) Replace `OPTGEOM-ENDGEOM`

block with keywords `TESTGEOM `

and `FINDSYM`

. The geometry block of the final input looks like this:

**Cu2O_Pn-3m_PBE0_TZVP_opt_done.d12**

Cu2O (Pn-3m) optimized geometry (DFT-PBE0/TZVP) CRYSTAL 0 0 0 224 4.31764698 2 29 0.000000000000E+00 0.000000000000E+00 0.000000000000E+00 8 2.500000000000E-01 2.500000000000E-01 2.500000000000E-01 TESTGEOM FINDSYM ENDGEOM

The rest of the input file is not changed at all.

5) You can now run this job with

runcrys Cu2O_Pn-3m_PBE0_TZVP_opt_done.d12

This will produce file `Cu2O_Pn-3m_PBE0_TZVP_opt_done.findsym.cif`

that can be viewed in any CIF viewer. Furthermore, the file `Cu2O_Pn-3m_PBE0_TZVP_opt_done.d12`

can be used for example for a harmonic frequency calculation.