TURBOMOLE manual
v8.0

21.4 Periodic Electrostatic Embedded Cluster Method

21.4.1 General Information

The Periodic Electrostatic Embedded Cluster Method (PEECM) functionality [377] provides electronic embedding of a finite, quantum mechanical cluster in a periodic, infinite array of point charges. It is implemented within HF and DFT energy and gradient TURBOMOLE modules: dscf, grad, ridft, rdgrad, and escf. Unlike embedding within a finite set of point charges the PEEC method always yields the correct electrostatic (Madelung) potential independent of the electrostatic moments of the point charges field. It is also significantly faster than the traditional finite point charges embedding.

21.4.2 Theoretical Background

Generally, the PEEC method divides the entire, periodic and infinite system into two parts, the inner (I) part, or so called cluster, and the outer (O) part which describes its environment. Thus, unlike "true" periodic quantum mechanical methods, PEECM primarily aims at calculations of structure and properties of localized defects in dominantly ionic crystals. The innermost part of the cluster is treated quantum mechanically (QM), whereas in the remaining cluster part cations are replaced by effective core potentials (ECPs) and anions by ECPs or by simply point charges. Such an "isolating" outer ECP shell surrounding the actual QM part is necessary in order to prevent artificial polarization of the electron density by cations which would otherwise be in a direct contact with the QM boundary. The outer part or the environment of the cluster is described by a periodic array of point charges, representing cationic and anionic sites of a perfect ionic crystal.

The electronic Coulomb energy term arising from the periodic field of point charges surrounding the cluster has the following form \[\begin{equation*} J = \sum_{\mu \nu} \sum_k^{N \in \text{UC}} \sum_{\vec L \in \text{O}}^{\infty} \\ D_{\mu \nu} q_k \int \frac{\mu(\vec r)\nu(\vec r)} {|\vec r - \vec R_k - \vec L|} d \vec r \text{,} \end{equation*}\] where UC denotes the unit cell of point charges, \(D_{\mu \nu}\) are elements of the density matrix, \(\mu\), \(\nu\) are basis functions, \(q_k\), \(\vec R_k\) denote charges and positions of point charges, and \(\vec L\) denote direct lattice vectors of the outer part O. It is evaluated using the periodic fast multipole method (PFMM) [378] which, unlike the Ewald method [379], defines the lattice sums entirely in the direct space. In general, PFMM yields a different electrostatic potential than the Ewald method, but the difference is merely a constant shift which depends on the shape of external infinite surface of the solid (i.e. on the way in which the lattice sum converges toward the infinite limit). However, this constant does not influence relative energies which are the same as obtained using the Ewald method, provided that the total charge of the cluster remains constant. Additionally, since the electrostatic potential within a solid is not a well defined quantity, both the absolute total energies and orbital energies have no meaning (i.e. you cannot compare energies of neutral and charged clusters!).

21.4.3 Calculation Setup

There are three key steps in setting up a PEECM calculation. In the first step the periodic field of point charges has to be defined by specifying the point charges unit cell. Next step is the definition of the part infinite of point charges field that will be replaced by the explicit quantum mechanical cluster. Finally, the quantum mechanical cluster together with surrounding ECPs representing cationic sites as well as point charges representing anions is defined and put in place of the point charges. The input preparation steps can be summarized as follows

  1. Dimensionality of the system is specified by the keyword periodic in the $embed section: periodic 3 means a bulk three-dimensional system, periodic 2 denotes a two-dimensional surface with an aperiodic \(z\) direction.

  2. Definition of the unit cell of periodic point charges field is specified in the subsections cell and content of the $embed section.

  3. Definition of the values of the point charges by specifying a charge value per species, using the subsection charges, or a charge value for each point charge, using the subsection ch_list. Note that only one of the subsections can be defined.

  4. Definition of the part of point charges field that will be replaced by the QM cluster together with the isolating shell (ECPs, explicit point charges) is specified in the subsection cluster of the $embed section.

  5. Definition of the quantum mechanical cluster as well as the surrounding ECPs and anionic point charges is included in the usual $coord section.

The following two examples show the definition of the point charges unit cells.

Example 1. Ca\(_4\)F\(_{19}\) cluster embedded in bulk CaF\(_2\)

In this example a QM cluster with the composition Ca\(_4\)F\(_{19}\), surrounded by 212 ECPs and 370 explicit point charges, representing Ca\(^{2+}\) cations and F\(^-\) anions is embedded in a periodic field of point charges (+2 for Ca and -1 for F) corresponding to the CaF\(_2\) fluorite lattice.

First, the program has to know that this is a three-dimensional periodic system. This is specified by the keyword periodic 3, meaning periodicity in three dimensions. The dimensions of the unit cell for bulk CaF\(_2\) are given in the subsection cell of the $embed keyword. By default, the unit cell dimensions are specified in atomic units and can be changed to Å using cell ang. The positions of the point charges in the unit cell are specified in the subsection content. In this example positions are given in fractional crystal coordinates (content frac). You can change this by specifying content for Cartesian coordinates in atomic units or content ang for Cartesian coordinates in Å. The values of point charges for Ca and F are given in the subsection charges.

$embed
periodic  3
cell
  10.47977 10.47977 10.47977 90.0 90.0 90.0
content frac
  F    0.00   0.00   0.00
  Ca  -0.25  -0.75  -0.75
  F    0.50  -0.50   0.00
  F    0.50   0.00  -0.50
  F    0.00  -0.50  -0.50
  F    0.50  -0.50  -0.50
  F    0.00   0.00  -0.50
  F    0.50   0.00   0.00
  F    0.00  -0.50   0.00
  Ca  -0.25  -0.25  -0.25
  Ca   0.25  -0.75  -0.25
  Ca   0.25  -0.25  -0.75
end
...
charges
  F   -1.0
  Ca   2.0
end

The above input defines a periodic, perfect, and infinite three-dimensional lattice of point charges corresponding to the bulk CaF\(_2\) structure. In order to use this lattice for PEECM calculation we have to make “space” for our QM cluster and the isolating shell. This is done by specifying the part of the lattice that is virtually removed from the perfect periodic array of point charges to make space for the cluster. The positions of the removed point charges are specified in the subsection cluster of the $embed keyword. Note, that the position of the QM cluster and the isolating shell must exactly correspond to the removed part of the crystal, otherwise positions of the cluster atoms would overlap with positions of point charges in the periodic lattice, resulting in a “nuclear fusion”.

cluster
  F    0.00000000000000   0.00000000000000   0.00000000000000
  Ca  -2.61994465796043  -2.61994465796043  -2.61994465796043
  Ca   2.61994465796043  -2.61994465796043   2.61994465796043
  Ca   2.61994465796043   2.61994465796043  -2.61994465796043
  Ca  -2.61994465796043   2.61994465796043   2.61994465796043
  F   -5.23988931592086   0.00000000000000   0.00000000000000
  F    0.00000000000000   0.00000000000000  -5.23988931592086
  F    5.23988931592086   0.00000000000000   0.00000000000000
  F    0.00000000000000  -5.23988931592086   0.00000000000000
  F    0.00000000000000   0.00000000000000   5.23988931592086
  F    0.00000000000000   5.23988931592086   0.00000000000000
  F   -5.23988931592086  -5.23988931592086   0.00000000000000
  F   -5.23988931592086   0.00000000000000  -5.23988931592086
  F   -5.23988931592086   0.00000000000000   5.23988931592086
  F   -5.23988931592086   5.23988931592086   0.00000000000000
  F    5.23988931592086  -5.23988931592086   0.00000000000000
  ...

repeated for Ca\(_{216}\)F\(_{389}\)

end

By default, the positions of point charges are specified in atomic units as Cartesian coordinates. You can change this by specifying cluster frac for fractional crystal coordinates or cluster ang for Cartesian coordinates in Å.

Finally, you have to specify the coordinates of the QM cluster along with the surrounding ECPs representing cationic sites and explicit point charges representing anions. This is done in the usual way using the $coord keyword.

$coord
    0.00000000000000      0.00000000000000      0.00000000000000      f
   -2.86167504097169     -2.86167504097169     -2.86167504097169      ca
    2.86167504097169      2.86167504097169     -2.86167504097169      ca
   -2.86167504097169      2.86167504097169      2.86167504097169      ca
    2.86167504097169     -2.86167504097169      2.86167504097169      ca
    0.00000000000000     -5.24009410923923      0.00000000000000      f 
   -5.24009410923923      0.00000000000000      0.00000000000000      f 
    0.00000000000000      5.24009410923923      0.00000000000000      f 
    0.00000000000000      0.00000000000000     -5.24009410923923      f 
    5.24009410923923      0.00000000000000      0.00000000000000      f 
    0.00000000000000      0.00000000000000      5.24009410923923      f 
    0.00000000000000     -5.24009410923923     -5.24009410923923      f 
   -5.24009410923923     -5.24009410923923      0.00000000000000      f 
    5.24009410923923     -5.24009410923923      0.00000000000000      f 
    0.00000000000000     -5.24009410923923      5.24009410923923      f 
    0.00000000000000      5.24009410923923     -5.24009410923923      f 
 ...

repeated for Ca\(_{216}\)F\(_{389}\)

$end

This is the standard TURBOMOLE syntax for atomic coordinates. The actual distinction between QM cluster, ECP shell, and explicit point charges is made in the $atoms section.

$atoms
f  1,6-23 
   basis =f def-TZVP
ca 2-5   
   basis =ca def-TZVP
ca 24-235 
   basis =none 
   ecp   =ca ecp-18 hay & wadt
f  236-605    
   basis =none
   charge= -1.00000000

In the example above the F atoms 1 and 6-23 as well Ca atoms 2-5 are defined as QM atoms with def-TZVP basis sets. The Ca atoms 24-235 are pure ECPs and have no basis functions (basis =none) and F atoms 236-605 are explicit point charges with charge -1, with no basis functions and no ECP.

This step ends the input definition for the PEECM calculation.

Example 2. Al\(_8\)O\(_{12}\) cluster embedded in \(\alpha\)-Al\(_2\)O\(_3\) (0001) surface

In this example a QM cluster with the composition Al\(_8\)O\(_{12}\), surrounded by 9 ECPs representing Al\(^{3+}\) cations is embedded in a two-dimensional periodic field of point charges (+3 for Al and -2 for O) corresponding to the (0001) surface of \(\alpha\)-Al\(_2\)O\(_3\).

As in the first example, the program has to know that this is a two-dimensional periodic system and this is specified by the keyword periodic 2. The dimensions of the unit cell for the (0001) \(\alpha\)-Al\(_2\)O\(_3\) surface are given in the subsection cell of the $embed keyword. The aperiodic direction is always the \(z\) direction, but you have to specify the unit cell as if it was a 3D periodic system. This means that the third dimension of the unit cell must be large enough to enclose the entire surface in this direction. The unit cell dimensions are specified in Å using cell ang. The positions of the point charges in the unit cell are specified as Cartesian coordinates in Å (content ang). The values of point charges for Al and O are given in the subsection charges.

$embed
periodic 2
cell angs
  4.8043    4.8043   24.0000   90.0000   90.0000  120.0000
content ang
  Al    2.402142286    1.386878848    5.918076515
  Al   -0.000013520   -0.000003382    7.611351967
  Al   -0.000008912    2.773757219    8.064809799
  Al    2.402041674    1.386946321    0.061230399
  Al   -0.000005568   -0.000003223   10.247499466
  Al    2.402137518    1.386872172    9.977437973
  Al    0.000000070    2.773757935    5.390023232
  Al    0.000006283   -0.000005607    3.696748018
  Al    2.402151346    1.386879444    3.243290186
  Al    0.000100868    2.773690462   11.246870041
  Al   -0.000001982   -0.000005796    1.060600400
  Al    0.000004853    2.773764610    1.330662251
  O    -0.731205344    1.496630311    6.749288559
  O     0.743527174    1.296469569    8.957922935
  O     1.588027477    0.104536049   11.127140045
  O     1.471626759    2.779079437    6.749288559
  O     3.309734344   -0.004341011    8.957920074
  O     3.919768333    1.323050499   11.127141953
  O    -0.740424335    4.045563698    6.749289513
  O    -1.651123047    2.868478537    8.957910538
  O     1.698525310    2.733071804   11.127161026
  O     3.133347750    2.664006472    4.558811665
  O     1.658615232    2.864167213    2.350177050
  O     0.814115047    4.056100845    0.180959582
  O     0.930515707    1.381557465    4.558811188
  O     1.494558096    0.004332162    2.350180149
  O    -1.517625928    2.837586403    0.180958077
  O     3.142566681    0.115072958    4.558810234
  O    -0.751034439    1.292158127    2.350189686
  O     0.703617156    1.427564979    0.180938885
end
...
charges
  O   -2.0
  Al   3.0
end

The above input defines a periodic, perfect, and infinite two-dimensional lattice of point charges corresponding to the (0001) \(\alpha\)-Al\(_2\)O\(_3\) surface. In order to use the lattice for PEECM calculation we have to make “space” for our QM cluster and the surrounding ECP shell. This is done by specifying the part of the lattice that is virtually removed from the perfect periodic array of point charges to make space for the cluster. The positions of the removed point charges are specified in the subsection cluster of the $embed keyword. Note, that the position of the QM cluster must exactly correspond to the removed part of the crystal, otherwise positions of the cluster atoms would overlap with positions of point charges in the periodic lattice, resulting in a “nuclear fusion”.

cluster ang
  Al    -0.000012482    5.547518253    9.977437973
  Al     2.402141094    6.934402943    8.064809799
  Al     2.402144432    4.160642624   10.247499466
  Al     4.804287434    5.547518253    9.977437973
  Al     2.402250767    6.934336185   11.246870041
  Al    -0.000005568    8.321288109   10.247499466
  Al     2.402137518    9.708164215    9.977437973
  Al     4.804294586    8.321288109   10.247499466
  O      0.907584429    4.156304836    8.957920074
  O      1.517618299    5.483696461   11.127141953
  O     -0.703624666    6.893717766   11.127161026
  O      3.145677090    5.457115650    8.957922935
  O      3.990177393    4.265182018   11.127140045
  O      0.751026928    7.029124260    8.957910538
  O      4.100675106    6.893717766   11.127161026
  O      0.743527174    9.617761612    8.957922935
  O      1.588027477    8.425827980   11.127140045
  O      3.309734344    8.316950798    8.957920074
  O      3.919768333    9.644342422   11.127141953
  O      5.555326939    7.029124260    8.957910538
  Al     4.804400921   11.094982147   11.246870041
  Al    -0.000008912    2.773757219    8.064809799
  Al    -2.402049065    6.934336185   11.246870041
  Al     4.804400921    2.773690462   11.246870041
  Al     2.402136564    4.160642624    7.611351967
  Al    -0.000013520    8.321288109    7.611351967
  Al    -0.000008912   11.095048904    8.064809799
  Al     7.206440926    6.934402943    8.064809799
  Al     4.804286480    8.321288109    7.611351967
end

The positions of point charges are specified in Å as Cartesian coordinates.

Finally, you have to specify the coordinates of the QM cluster along with the surrounding ECPs. This is done in the usual way using the $coord keyword.

$coord
   -0.00002358760000     10.48329315900000     18.85463057110000      al
    4.53939007480000     13.10412613690000     15.24028611330000      al
    4.53939638280000      7.86247730390000     19.36497297520000      al
    9.07879006320000     10.48329315900000     18.85463057110000      al
    4.53959732680000     13.10399998250000     21.25351019750000      al
   -0.00001052200000     15.72496001430000     19.36497297520000      al
    4.53938331720000     18.34577677080000     18.85463057110000      al
    9.07880357850000     15.72496001430000     19.36497297520000      al
    1.71508649490000      7.85428007030000     16.92802041340000      o
    2.86788376470000     10.36268741690000     21.02725683720000      o
   -1.32965829240000     13.02724227310000     21.02729288000000      o
    5.94446987180000     10.31245694970000     16.92802581990000      o
    7.54034461170000      8.06002818410000     21.02725323160000      o
    1.41923561090000     13.28312353520000     16.92800239300000      o
    7.74915508620000     13.02724227310000     21.02729288000000      o
    1.40506312580000     18.17494056150000     16.92802581990000      o
    3.00093786570000     15.92251179600000     21.02725323160000      o
    6.25449323900000     15.71676368210000     16.92802041340000      o
    7.40729073370000     18.22517102690000     21.02725683720000      o
   10.49804944110000     13.28312353520000     16.92800239300000      o
    9.07900452260000     20.96648359440000     21.25351019750000      al
   -0.00001684120000      5.24164297480000     15.24028611330000      al
   -4.53921616520000     13.10399998250000     21.25351019750000      al
    9.07900452260000      5.24151682240000     21.25351019750000      al
    4.53938151440000      7.86247730390000     14.38337475740000      al
   -0.00002554910000     15.72496001430000     14.38337475740000      al
   -0.00001684120000     20.96660974680000     15.24028611330000      al
   13.61820356690000     13.10412613690000     15.24028611330000      al
    9.07878826040000     15.72496001430000     14.38337475740000      al
$end

This is the standard TURBOMOLE syntax for atomic coordinates. The actual distinction between QM cluster and ECP shell is made in the $atoms section.

$atoms
al 1-8 
   basis =al def-SV(P)
o  9-20
   basis =o def-SV(P)
al 21-29
   basis =none 
   ecp   =al ecp-10 hay & wadt

In the example above the Al atoms 1-8 and O atoms 9-20 are defined as QM atoms with def-SV(P) basis sets. The Al atoms 21-29 are pure ECPs and have no basis functions (basis =none).

This step ends the input definition for the PEECM calculation.