This part of the guide describes the calculation of harmonic frequencies with CRYSTAL, using Cu2O as an example.
During the geometry optimization, we used the first derivatives of the energy with respect to atomic displacements to locate a stationary point on the potential energy surface (derivatives are zero). In a harmonic frequency calculation, we calculate the second derivatives of the enegy with respect to atomic displacements. This yields us the vibrational frequencies of the system and tells us the nature of the stationary point (minimum, maximum, saddle point). In CRYSTAL, the second derivatives are obtained with numerical differentiation of the analytical first derivatives.
Prerequisites
You need the file Cu2O_Pn-3m_PBE0_TZVP_opt_done.d12
created in the geometry optimization of Cu2O (section Create a CIF of the optimized structure). This input file contains the optimized geometry of Cu2O.
1) Create a frequency calculation input from the input with the optimized geometry:
cp Cu2O_Pn-3m_PBE0_TZVP_opt_done.d12 Cu2O_Pn-3m_PBE0_TZVP_fq.d12
2) Replace the FINDSYM and TESTGEOM keywords in the input with:
FREQCALC
NUMDERIV
2
ENDFREQ
NUMDERIV 2
means central differentiation, which is numerically more reliable than simple one-sided differentiation.
3) The default SCF convergence criterion in FREQCALC is 10–10 a.u. which is very tight. For a Γ-point frequency calculation the criterion can typically be relaxed a little bit. Add after TOLINTEG:
TOLDEE
9
The complete input is:
Cu2O (Pn-3m) frequency calculation (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 FREQCALC NUMDERIV 2 ENDFREQ 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 TOLDEE 9 MAXCYCLE 100 FMIXING 80 EXCHSIZE 30000000 BIPOSIZE 30000000 END
Running the frequency calculation
Submit the job normally using jsub:
jsub -np 8 crystal Cu2O_Pn-3m_PBE0_TZVP_fq.d12
Analyzing the results
An example output file is available here: Cu2O_Pn-3m_PBE0_TZVP_fq_only.o.
The vibrational modes can be visualized for example with Jmol or CRYSPLOT (see Aalto Solid State Chemistry Wiki).
The output file lists the vibrational frequencies as follows: :
MODES EIGV FREQUENCIES IRREP IR INTENS RAMAN
(HARTREE**2) (CM**-1) (THZ) (KM/MOL)
1- 3 -0.7906E-21 0.0000 0.0000 (F1u) A ( 0.00) I
4- 6 0.1617E-06 88.2545 2.6458 (F2u) I ( 0.00) I
7- 8 0.2594E-06 111.7855 3.3512 (Eu ) I ( 0.00) I
9- 11 0.4435E-06 146.1578 4.3817 (F1u) A ( 0.00) I
12- 12 0.2447E-05 343.3387 10.2930 (Au ) I ( 0.00) I
13- 15 0.5124E-05 496.7951 14.8935 (F2g) I ( 0.00) A
16- 18 0.7994E-05 620.5167 18.6026 (F1u) A ( 0.00) I
Modes 1–3 with zero frequency at the Γ-point are translational (acoustic) phonon modes. Modes 4–18 with non-zero frequencies are optical modes. For 3D solids, we always have 3 acoustic modes and 3N−3 optical phonon modes. There are many degenerate vibrational modes with the same frequency due to the high symmetry of the Cu2O crystal.
Optical modes 9–11 and 16–18 are IR-active. Modes 13–15 are Raman active. We will calculate the intensities for these modes later.
A thermodynamic analysis is printed under THERMODYNAMIC FUNCTIONS WITH VIBRATIONAL CONTRIBUTIONS
, but normally Γ-point frequencies are not enough for accurate thermodynamics (full phonon dispersion relations are needed). Total Gibbs free energy is the term EL+E0+ET+PV-TS
.