In this example, we will optimize the geometry of Cu2O 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.

Key points concerning geometry optimizations

  • By default, CRYSTAL will optimize both the lattice vectors and the atomic coordinates. See section 6.3.1 in the manual for other options such as ATOMONLY, CELLONLY, or CVOLOPT.
  • CRYSTAL will always honor the space group symmetry during the optimization, unless the user deliberately breaks the symmetry before the optimization (for example with SYMMREMO keyword).
  • The optimization is a local optimization within the constraints of the space group symmetry. If you are interested in crystal structure prediction and finding global minima, have a look at the USPEX interface for CRYSTAL.

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. 

SHRINK is always given for the primitive cell

Note that SHRINK is given for the primitive cell. This is yet another reason, why you should first do a TESTGEOM run to determine the lattice vectors of the primitive cell. Only then you can determine a reasonable value for the SHRINK. It is possible that the cell edge lengths of a centered crystallographic are anisotropic, but the cell edge lengths of the primitive cell are isotropic!

Give the SHRINK keyword as follows:

SHRINK
0 SHRINK1
SHRINK1 SHRINK2 SHRINK3

The second line only needs to change if we are dealing with metals with partially occupied bands and need more accurate k-mesh for the determination of the Fermi energy (see section 2.3 in CRYSTAL17 manual).

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:

  1. You need to have a terminal connection to the server.
  2. 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 Cu2O: 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 Cu2O (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 Cu2O, 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.

You MUST manually remove unnecessary lattice parameters from the new input! crygeom prints out all lattice parameters, but for a cubic system, we remove all lattice parameters except a. For a hexagonal system, we would keep a and c, etc. 

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.


  • No labels