3.48 Fragment molecular orbital DFT calculations

The fragment molecular orbital (FMO or FO) scheme allows the efficient calculation of transfer matrix elements for transport calculations. The two reference states are constructed from isolated calculations of the respective fragments (donor and acceptor). This circumvents the electron delocalistion error, but neglects any interactions between the fragments.

Important: This functionality is not yet available for periodic systems.
Important: The $\delta+$ -FODFT scheme is not optimized for large systems.

FODFT flavours

Depending on the approximations done to implement the FODFT method we distinguish between three different flavours of FODFT. Please see [272] for a detailed description and assessment of each. Within FHIaims it is possible to use all different FODFT-schemes.

$H^{2n}@DA$

This is the classic implementation by [277]. The fragment calculations are always done for the neutral fragments and the Hamiltonian is constructed using a wrong number of electrons. $H_{ab}=H_{ba}$ even for hetero dimers.

To use this flavour, simply calculate neutral fragments and continue with the default FODFT options.

$H^{2n-1}@DA$ / $H^{2n+1}@D^{-}A^{-}$

This implementation was first used within the CPMD code[234]. Here, neutral (double anionic) fragment calculations are combined for hole (electron) transfer with a reset of the occupation number in the highest occupied molecular orbital. This leads to the correct number of electrons for the construction of the Hamiltonian.

To use this flavour, calculate neutral (or anionic) fragments and use the fo_flavour keyword with the option reset.

$H^{2n\pm 1}@D^{\pm}A$

This scheme was first implemented within FHIaims and uses charged fragment calculations. This results in slightly different approximations with improved accuracy in the electronic coupling values. Please see [272] for more details.

To use this flavour, calculate appropriately charged fragments and continue with the default FODFT options.

$\delta+$ -FODFT

The implementation of the $\delta+$ -FODFT scheme is not optimized for large systems and considered experimental. In our study[272] we found no significant improvement when using this scheme.

Theoretical Background

The following section intends to give a quick overview over the theory and the implementation within FHIaims. For detailed theory we refer to the before mentioned publications. The following theory is for the charged FODFT flavour ( $H^{2n\pm 1}@D^{\pm}A$ ).

In principle, the transfer matrix element for hole transport between the HOMO and the LUMO of two molecules is given by

\displaystyle H_{AB}=\langle\Psi^{D}_{\text{HOMO}}|\hat{H}_{KS}|\Psi^{A}_{LUMO% }\rangle.

(3.253)

Obtain the fragment wavefunctions

First, two standard-DFT calculations are done to obtain the fragment wavefunctions. Some additional output is needed for the subsequent combination of the fragments to the full system.

Combine the fragment wavefunctions and extract the transfer matrix element

The saved wavefunctions of both fragment calculations are used to obtain the full, non-interacting electron density of the combined system (D+A). Although FHIaims uses atom centered basis functions, it is not generally possible to re-use any restart file. If the ordering of the atoms in the geometry.in files is identical, it is possible to use a restart file for any geometry with translational symmetry. If the geometry is rotated, it is necessary to rotate the wave function.

For this non self-consistent density the Hamiltonian is calculated and used to determine $H_{AB}$ . We therefore use the non-interacting density, but the interacting Hamiltonian with the correct number of electrons. This is in contrast to other implementations, where only the Hamiltonian of the neutral fragment is used.

The first step is a Loewdin-Orthogonalization of the Kohn-Sham-eigenvector.

\displaystyle\mathbf{U^{\top}SU}=\mathbf{s}

(3.254)

\displaystyle\mathbf{X}=\mathbf{S}^{-\frac{1}{2}}=\mathbf{Us}^{-\frac{1}{2}}% \mathbf{U}^{\top}

(3.255)

	$\displaystyle\mathbf{C^{\prime}}$	$\displaystyle=\mathbf{X}^{-1}\mathbf{C}$
	$\displaystyle\mathbf{H^{\prime}}$	$\displaystyle=\mathbf{X^{\top}HX}$		(3.256)

In the next step, the Hamiltonian is transformed,

\displaystyle\mathbf{H^{\prime}}_{\text{st.}}=\mathbf{C^{\prime\top}H^{\prime}% C},

(3.257)

and the overlap of states is calculated:

\displaystyle\mathbf{S}_{\text{st.}}=\mathbf{C^{\top}\cdot C}.

(3.258)

In the final step, the transfer matrix element $H_{AB}$ is calculated with

\displaystyle H_{AB}=\left(1-\mathbf{S}{{}_{\text{st.}}(\text{a,b})}^{2}\right% )^{-1}\left[{\mathbf{H^{\prime}}_{\text{st.}}(\text{a,b})-\mathbf{S}_{\text{st% .}}(\text{a,b})\left(\mathbf{H^{\prime}}_{\text{st.}}(\text{a,a})+\mathbf{H^{% \prime}}_{\text{st.}}(\text{b,b})\right)/2}\right]

(3.259)

Important technical hints

When creating the input geometries for the calculations, the geometry.in must mirror the actual fragmentation. This means, the order of atoms has to be consistent for the dimer and the fragments.

Example:

Calculating the $H_{AB}$ for two methane molecules (CH ${}_{4}^{+}$ – CH₄)

Complete System, CH ${}_{4}^{+}$ –CH₄

#fragment 1
atom    12.7490 4.3034  0.0000  C
atom    13.8413 4.3037  0.0000  H
atom    12.3850 4.7804  -0.9128 H
atom    12.3850 4.8558  0.8693  H
atom    12.3851 3.2750  0.0435  H
#fragment 2
atom    -2.6946 4.4740  0.0002  C
atom    -1.6023 4.4740  0.0000  H
atom    -3.0586 3.6144  -0.5671 H
atom    -3.0586 5.3949  -0.4609 H
atom    -3.0585 4.4127  1.0278  H

Fragment 1, CH ${}_{4}^{+}$

atom    12.7490 4.3034  0.0000  C
atom    13.8413 4.3037  0.0000  H
atom    12.3850 4.7804  -0.9128 H
atom    12.3850 4.8558  0.8693  H
atom    12.3851 3.2750  0.0435  H

Fragment 2, CH₄

atom    -2.6946 4.4740  0.0002  C
atom    -1.6023 4.4740  0.0000  H
atom    -3.0586 3.6144  -0.5671 H
atom    -3.0586 5.3949  -0.4609 H
atom    -3.0585 4.4127  1.0278  H

Folder structure for FODFT-calculation

At least three separate calculations (fragment1, fragment2 and combination step) are needed for a complete FODFT run. To maintain a clean structure and avoid copying of files, a special scheme is enforced. The folders can have arbitrary names, but need to be in the same root directory. This allows easy reuse of fragment calculations for different dimer calculations.

- fodft\_something (main calculation folder)
|
+-- dimer_01
    |
    +-- control.in
    +-- geometry.in
+-- frag1/
    |
    +-- control.in
    +-- geometry.in
+-- frag2/
    |
    +-- control.in
    +-- geometry.in

$\delta^{+}$ FO-DFT - Inclusion of the local potential

In our implementation, an additional option is available to include interactions between the fragments by means of an embedding scheme, using the full local potential ( $V_{loc}$ ) of the embedded fragment. In contrast to the standard implementation, the fragment wavefunction $\Psi_{A}$ is polarized by the potential of the second fragment $\Psi_{B}$ during the SCF cycle.

The general procedure is:

1.
Calculation of fragment wavefunctions
1. (a)
  
  Calculate Fragment A, export potential A
2. (b)
  
  Import potential A, calculate Fragment B, export new potential B
3. (c)
  
  Import potential B, calculate Fragment A^∗∗
2.

Combination of fragment wavefunctions
3.

Determination of $H_{AB}$

**: Please note that this step might not be necessary if fragment A was choosen wisely, depending on the polarizability of the fragments.

The embedding can be used to iteratively converge towards a final solution, but as usual the first steps are the most important ones for the polarization.

IMPORTANT: Since this method requires the reassignment of grid points from the fragment calculations for the local potential, the current implementation is not optimized towards large systems.

Tags for general section of control.in:

Tag: fo_dft(control.in)

Usage: fo_dft type
Purpose: This is the central control keyword for fragment orbital DFT. Its main purpose is the control of the desired calculation step (fragment or combined system) and associated choices. Depending on type, further flags or lines are necessary. Those are explained below.
Options: Possible options are fragment and final.

fo_dft sub-tag: fragment(control.in)

Usage: fo_dft fragment
Purpose: This keyword indicates that a fragment for FODFT will be calculated. The restart information (restart.frag) and an additional output file (info.frag) with information for the recombination of the fragments for the final FODFT run will be written.

fo_dft sub-tag: final(control.in)

Usage: fo_dft final
Purpose: This keyword indicates that this is the final FODFT run, where to wavefunctions of both fragments are combined and the transfer matrix element is calculated.

IMPORTANT: This calculation can only be started after both fragments have been calculated!

Tag: fo_orbitals(control.in)

Usage: fo_orbitals state1 state2 range1 range2 type
Required input for the final FODFT step. Determines the actual states of interest for the calculation of the coupling elements.

•

state1, state2: [Integer] – Isolated fragment states for the determination of the respective matrix element.
•

range1, range2: [Integer] – Default 1. Gives a range for the selected states.
( $\text{state}1=5$ and $\text{range}1=2$ will include state 5 and state 6 of fragment 1.)
•

type: [String] Determines whether the matrix element for the spin-up (up or elec) or spin-down (dn or hole) Hamiltonian is calculated. In most cases, spin down means hole transport (excess hole on fragment 1 to fragment 2) and spin up means electron transport (excess electron on fragment 1 to fragment 2).
Example Zn⁺-Zn: fo_orbitals 15 15 1 1 hole
will output the transfer matrix elements for hole transfer between the LUMO of Zn⁺ (fragment 1) and the HOMO of Zn (fragment 2).

If not deactivated via fo_verbosity 0, the full $H_{ab}$ submatrix will be written to the file full_hab_submatrix and any matrix element can be retrieved from this file after the calculation.

Tag: fo_folders(control.in)

OPTIONAL Usage: fo_folders folder_frag1 folder_frag2
Purpose: This keyword allows the specification of custom names for the folders with the fragment calculation. If not set, standard names (frag1_00 and frag2_00) are used.

Tag: fo_flavour(control.in)

OPTIONAL Usage: fo_flavour flavour
Purpose: This keyword allows the specification of the FO-DFT flavour to be used for the calculation. See FODFT theory section for details. Options:

•

default - no occupations are modified. For $H^{2n}@DA$ and $H^{2n\pm 1}@D^{\pm}A$ flavours. This is the default.
•

reset - selects the $H^{2n-1}@DA/H^{2n+1}@D^{-}A^{-}$ scheme.

Note: There is no option yet to calculate $H_{ab}$ and $H_{ba}$ within one FHIaims run. If the system is heterogeneous, two dimer steps with appropriate fragment ordering and charges are necessary!

Tag: fo_verbosity(control.in)

Usage: fo_verbosity level
Purpose: Control the verbosity of the FO-DFT output. Options:

•

0: Write minimal output to output file. No retrieval of arbitrary couplings without restarting the dimer calculation possible.
•

1: Write calculated couplings to output file and the full $H_{ab}$ matrix to the file full_hab_submatrix. This is the default.
•

2: Write all elements for selected states from eq. 3.259 ( $S_{st},H^{\prime}_{st},\ldots$ ) to aims output.

Tag: fo_deltaplus(control.in)

Usage: fo_deltaplus .true.
Purpose: Enables the polarized $\delta^{+}$ fragment calculations.

With this keyword, the full local potential of the actual fragment (Hartree + nuclei) will be written to files (output_potential_0000*). In addition, if previous potential files are found within the calculation folder, they will be read in and added as an external potential to the current calculation.
Important: In order to use this improved scheme, the (empty) atom positions of the other fragment have to be included in the geometry.in (see 3.48).
Important: The number of cores (determined by mpirun -n or similar) has to be the same for both fragments, or the calculation will fail!
Note: Due to differences in the partitioning of the FHIaims grid between these runs it is necessary to reassign the potential to each grid point. Depending on the system size and the number of MPI processes (communication!), this will take some time.

Additional tags for polarized FODFT in geometry.in and control.in:

The following changes are only necessary for polarized FODFT (fo_deltaplus). In order to allow a polarization of the fragment wavefunction, the atom centered integration grids have to be extended to the region of the embedded fragment.

geometry.in:

The atomic positions of the second, embedded fragment must be included as empty sites with the tag empty (instead of atom). In addition, even the minimal basis functions have to be disabled in the species_defaults for this atom type in the control.in (see next section). To distinguish between real and empty atoms, the species name has to be changed.

The following example for Zn ${}_{2}^{+}$ shows the 3 different geometry.in-files.

Fragment 1, Zn⁺

atom 0.00 0.00 0.00 Zn
initial_charge 1
empty 0.00 0.00 5.00 Zn_empty

Fragment 2, Zn

empty 0.00 0.00 0.00 Zn_empty
#initial_charge 1 #no charge for emtpy atoms
atom 0.00 0.00 5.00 Zn

Complete System, Zn ${}_{2}^{+}$

atom 0.00 0.00 0.00 Zn
initial_charge 1
atom 0.00 0.00 5.00 Zn

control.in:

For every atom type of the embedded system, the species_data part has to be included and the species name must be changed to species_empty. Disable the minimal basis with the keyword species include_min_basis set to .false. and comment all tier basis function lines, as showed in the following example for Zn.


[...]
  species        Zn_empty
#     global species definitions
    nucleus             30
    mass                65.409
#
    l_hartree           4
#
    cut_pot             3.5          1.5  1.0
    basis_dep_cutoff    1e-4
#
    include_min_basis .false.

[...]

#  "First tier" - improvements: -270.82 meV to -12.81 meV
     #hydro 2 p 1.7
     #hydro 3 s 2.9
     #hydro 4 p 5.4
     #hydro 4 f 7.8
     #hydro 3 d 4.5
##  "Second tier" - improvements: -3.35 meV to -0.82 meV
     #hydro 5 g 10.8
     #hydro 2 p 2.4
     #hydro 3 s 6.2
     #hydro 3 d 3