CASEA的奇思妙想

Predicting EPR

发表于 更新于 分类于
本文字数： 10k 阅读时长 ≈ 9 分钟

预测EPR特性和光谱

Predicting EPR

使用ORCA预测EPR特性和光谱，其中甚至包括许多复杂的问题（自旋轨道耦合和其他相对论效应）

The CH3 radical

甲基自由基
甲基自由基是一种平面分子，其平面性部分已被EPR光谱证实。
It has a g of 2.0025 and hyperfine couplings (HFC) Axx=−31.4,Ayy=−93.7, and Azz=−60.30，let’s try to predict that from first principles!

优化结构

1
2
3
4
5
6
7
!B3LYP DEF2-TZVP OPT
* xyz 0 2
C   3.20276693914010      2.91777006829308     -1.46614904250634
H   4.21682401439677      3.22254500904168     -1.25560882817865
H   2.69109288861722      3.31313693290047     -2.33037576806431
H   2.69057615784592      2.24676798976475     -0.79431636125068
*

你将得到优化后的三角平面坐标

下一步是计算EPR性质，通过设置%EPRNMP计算。
重点：这些是ORCA输入的唯一选项，应该几何优化之后
An input to compute the g tensor and HFCs for the hydrogen atoms, using the B3LYP functional could be:

1
2
3
4
5
6
!B3LYP EPR-II AUTOAUX
* XYZFile 0 2 CH3_opt.xyz
%EPRNMR
        GTENSOR   TRUE
        NUCLEI    = ALL H {AISO, ADIP, AORB}
END
  • The GTENSOR flag must be set to TRUE in order to compute it, the default is FALSE.
  • The NUCLEI defines the atoms for the HFC calculation. ALL H calculates the HFC on all hydrogens, or use ALL N, ALL C and so on for different atoms. You can also use NUCLEI = 1,5,8, to give one list per atom type with the atom numbering staring from 1.
  • The AISO,ADIPand AORB flags request the isotropic part of the HFC, its dipolar part and 2nd order contribution from SOC respectively.
  • The AORB term is costly, but contributes very little to organic molecules without heavy atoms and could be neglected in those cases.
  • A large basis is recommended for EPR properties. Here we used the special basis EPR-II, for more info one can check the full documentation.
  • The spin-orbit coupling integrals are calculated using the SOMF(1X)method by default, but more efficient options are also available.

If you are using a list of atoms, like NUCLEI = 1,5,8, there should be one list for each atom type!

The output file

在单点能计算成功后，第一步是计算用于处理SOC的积分

1
2
3
4
5
6
7
8
9
10
11
12
13
------------------------------------------------------------------------------
                      ORCA SPIN-ORBIT COUPLING CALCULATION
------------------------------------------------------------------------------

GBW file                                    ... epr.gbw
Input density                               ... epr.scfp
Output integrals                            ... epr
Operator type                               ... Mean-field/Effective potential
   One-Electron Terms                       ... 1
   Coulomb Contribution                     ... 2
   Exchange Contribution                    ... 3
   Correlation Contribution                 ... 0
   Maximum number of centers                ... 4
1
2
3
4
5
6
7
8
9
10
------------------------------------------------------------------------------
                              ORCA EPR/NMR CALCULATION
------------------------------------------------------------------------------

GBWName                      ... epr.gbw
Electron density file        ... epr.scfp
Spin density file            ... epr.scfr
Spin-orbit integrals         ... epr
Origin for angular momentum  ... Center of electronic charge
Coordinates of the origin    ...    6.05064281    5.51881078   -2.76743570 (bohrs)

Please notice that the “Origin for angular momentum” here is “Center of electronic charge”, and that means these results are not gauge-invariant and can change depending on your coordinates.

The calculation of the integrals involving GIAOs is much more costly than the regular ones, so be prepared for longer calculation times.

Usually, with a sufficiently large basis, the results should be almost invariant, but if one wants the exact results, gauge-invariant atomic orbitals (GIAOs) are needed . These can be constructed from all basis sets, and to use them, just add ORI GIAO to the %EPRNMR section:

1
2
3
4
5
%EPRNMR
        GTENSOR   TRUE
        NUCLEI    = ALL H {AISO, ADIP, AORB}
        ORI       GIAO
END

g Tensor

If the following CP-SCF converges, you will get the g tensor matrix:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
-------------------
ELECTRONIC G-MATRIX
-------------------

      The g-matrix:
                   2.0027537    0.0001534    0.0001065
                   0.0001534    2.0024625   -0.0002479
                   0.0001065   -0.0002479    2.0026476

      gel          2.0023193    2.0023193    2.0023193
      gRMC        -0.0001271   -0.0001271   -0.0001271
      gDSO(tot)    0.0000289    0.0000680    0.0000680
      gPSO(tot)    0.0000034    0.0005594    0.0005595
                  ----------   ----------   ----------
      g(tot)       2.0022245    2.0028196    2.0028196 iso=  2.0026213
      Delta-g     -0.0000948    0.0005003    0.0005004 iso=  0.0003020

Here the components of the g tensor are specified and the total components for the x, y and z directions are printed after g(tot), with the isotropic value afterwards. The Delta-g shown is the difference from the free-electron g value.

As you can see, the predicted value is already quite close to the experimental value, considering we completely ignored any matrix effects and its largest deviation from the free-electron comes due to SOC effects, as suggested by Gordy.

Hyperfine couplings

After the g tensor calculation, comes the HFC for each atom:

1
2
3
4
5
6
7
8
9
10
-----------------------------------------
ELECTRIC AND MAGNETIC HYPERFINE STRUCTURE
-----------------------------------------

-----------------------------------------------------------
Nucleus   1H : A:ISTP=    1 I=  0.5 P=533.5514 MHz/au**3
               Q:ISTP=    2 I=  1.0 Q=  0.0029 barn
-----------------------------------------------------------
Solving for nuclear perturbation:
Forming RHS of the CP-SCF equations         ...

Here “I=” indicates the spin of that given atom and “P=” is the proportionality factor μeμNβeβNμeμNβeβN used. After that, another CP-SCF solution is requested and the HFC information is printed:

1
2
3
4
5
6
7
8
9
10
11
12
13
Raw HFC matrix (all values in MHz):
-----------------------------------
              -30.0552              10.8247               7.4998
               10.8247             -77.0049              20.6090
                7.4998              20.6090             -92.4037

A(FC)         -66.5005             -66.5005             -66.5005
A(SD)          41.0835              -0.8524             -40.2311
A(ORB+DIA)      0.0114              -0.0006               0.0271    A(PC) =   0.0126
A(ORB)          0.0065               0.0006               0.0258    A(PC) =   0.0109
A(DIA)          0.0050              -0.0012               0.0013    A(PC) =   0.0017
            ----------           ----------           ----------
A(Tot)        -25.4056             -67.3536            -106.7045    A(iso)=  -66.4879

Again, the individual components are printed and the A(iso) is shown in the end. As you can see, the results are in quite good agreement with the experimental data.

At the very end, the ORCA EULER output is printed, with information about the relative angle between the g-tensor and the HFCs, which have an arbitrary relationship from the calculation:

1
2
3
4
5
6
7
8
9
10
11
12
-------------------------------------------------------------------------
          Euler rotation of hyperfine tensor to g-tensor
-------------------------------------------------------------------------

-------------------------------------------------------------------
 Atom  |   Alpha    Beta    Gamma   |   Ax       Ay       Az
       |          [degrees]         |           [MHz]
-------------------------------------------------------------------
  1H      -90.1     42.7     89.9     -67.35   -25.41  -106.70
  2H       90.0     17.4    -90.0     -67.36   -25.41  -106.71
  3H      -90.1     12.7     90.0     -67.36  -106.72   -25.42
-------------------------------------------------------------------

You can use the orca_euler software later if you want to rotate any of these for any specific application.

Comparison to experiment

To summarize our results, here is a table comparing the calculated and predicted data:
With these, one can use tools such as EasySpin, which already has a specific function for reading orca outputs and simulate the EPR spectrum. The comparison between the theory prediction and an EPR spectra taken at 4K in a CH4 matrix is:

Of course differences here are expected due to impurities and matrices effects and the graphic was displaced to match the gg.

The HFCs are quite sensitive to the geometry of the molecules! Be sure to have a good one before computing them.

Another example, the naphthalene anion radical

In aprotic solvents and in the presence of alkali metals, polyaromatics such as naphthalene can reduce to form radical anions

These will present nice EPR spectra, with two sets of four protons and thus more complex lines. To predict its spectrum, again one can start by optimizing the geometry of the anion:

1
2
!RI-B2PLYP DEF2-TZVP DEF2-TZVP/C TIGHTOPT
* xyzfile -1 2 nap.xyz

We used a double-hybrid functional to add RI-MP2 correlation and get a better density and TIGHTOPT to obtain a good symmetrical geometry.

Now we run an EPR calculation, with the specialized EPR-III basis and RI-SOMF(1X) to speed-up the calculation of the integrals for the SOC part:

1
2
3
4
5
6
!RI-B2PLYP EPR-III AUTOAUX RI-SOMF(1X)
* xyzfile -1 2 nap_optimized.xyz
%eprnmr
        gtensor true
        nuclei = all H {aiso, adip, aorb}
end

One can average the two sets of similar HαHα and HβHβ couplings (calculated from the MP2 relaxed density ) and use again some software like EasySpin to compute the predicted spectrum and compare it to the experiment. This will look like:

Visualizing the g tensor orientation

You can use Avogadro to visualize the orientation of the g tensor. Here is how it works: ORCA prints a basename.gori.xyz file, where dummy atoms are used to represent these vectors, He for the origin, Ne for the direction of gz, Ar for gy and Kr for gx

  1. Open this file and select the “Align” option in Avogadro, on the right of the “E” used for optimization.
  2. Now select the He-Ne bond and align the x axis and repeat in order for the y and z axis.
  3. Delete the noble gas atoms.
  4. Click on “Axes”, the first option of the “Display Settings”. Here it is!

You can increase the size of the vectors by clicking on the “Axes” options. The order is red for x, green for y and blue for z (remember the RGB sequence).

Zero-field Splitting (ZFS)

If you have systems with S>1/2, you might be interested in ZFS values. The ORCA_EPRNMR module can also be used to predict those. Let’s try the example of cyclopentadienylidene carbene, which has known D value of 0.408cm−1 and a E/Dratio of 0.046 . First, the geometry is optimized using:

1
2
!B3LYP DEF2-TZVP OPT
* XYZFILE 0 3 carb.xyz
../_images/carbene.png

The ZFS can then be computed on simple DFT level using:

1
2
3
4
5
6
!B3LYP EPR-II AUTOAUX UNO
* XYZFile 0 3 carb_optimized.xyz
%EPRNMR DTENSOR SSANDSO
        DSS     UNO
        DSOC    CP
END

The DTENSOR flag asks to include the spin-spin and the spin-orbit components of the D (SSANDSO), and the DSS and DSOC control the specific algorithms. For more information, consult the ORCA manual.
After the calculation of the SOC-related integrals, one should see the header:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
------------------------------------------------------------------------------
                                ORCA EPR/NMR CALCULATION
------------------------------------------------------------------------------

GBWName                      ... ZFS.gbw
Electron density file        ... ZFS.scfp
Spin density file            ... ZFS.scfr
Spin-orbit integrals         ... ZFS
Multiplicity                 ... 3
g-tensor                     ... 0
NMR chemical shifts          ... 0
D-tensor                     ... spin-spin and spin-orbit
D(SS)-algorithm              ... spin-density from UNOs
D(SOC)-algorithm             ... CP (=Coupled-perturbed)

explaining the details of the calculation. The actual result is below the CP-SCF:

1
2
D   =     0.425943  cm**-1
E/D =     0.015874

which is again close to the experiment, considering that we are using DFT and a single-reference method.

This was intended to be a simple demonstration. In general, much better results for ZFS can be obtained from multi-reference theories such as CASSCF or MRCI. Nonetheless, it is still quite good and is easily applicable to large systems.

Structures

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
# Napthalene
C          2.41963        0.74477        0.00000
C          2.44765       -0.64677       -0.00000
C          1.25295       -1.36872        0.00000
C          0.01423       -0.70692        0.00000
C         -0.01423        0.70692       -0.00000
C          1.19684        1.41804       -0.00000
C         -1.25295        1.36871        0.00000
C         -2.44765        0.64677       -0.00000
C         -2.41963       -0.74477        0.00000
C         -1.19684       -1.41804       -0.00000
H          3.34913        1.30792       -0.00000
H          3.39906       -1.17205        0.00000
H          1.29505       -2.45547       -0.00000
H          1.19516        2.50560        0.00000
H         -1.29505        2.45547       -0.00000
H         -3.39906        1.17205        0.00000
H         -3.34913       -1.30792       -0.00000
H         -1.19516       -2.50561        0.00000

# Cyclopentadienylidene carbene
C         -1.73238        3.09526        0.00000
C         -2.43358        1.94496       -0.00000
C         -0.30953        2.77063        0.00000
C         -0.17373        1.43249       -0.00000
C         -1.48252        0.90001       -0.00000
H         -2.13310        4.09634        0.00000
H         -3.50614        1.82921       -0.00000
H          0.48447        3.50033        0.00000
H          0.74363        0.86480       -0.00000