.. tip:: All input files can be downloaded: :download:`Files `. .. tip:: For more information of this section, please refer to these pages: - :doc:`./tso1` - :doc:`./msdft1` - :doc:`./msdft2` - :doc:`./msdft3` - :doc:`../keywords/nosi` - :doc:`../keywords/msdft` MSDFT (4): Charge Transfer Excitations ================================================================================ .. contents:: :local: This tutorial will lead you step by step to study excited states using Multi-State Density Functional Theory (MSDFT). MSDFT is a powerful method for studying excited states. It optimizes excited states using TSO-DFT method, so it is **free of orbital relaxation problem** like in TDDFT. In this section, you will see that MSDFT can give **much more reasonable** results than TDDFT for excited states. In :doc:`tso1`, we have introduced how to use TSO-DFT to study excited states. In :doc:`msdft1`, we have introduced how to use MSDFT to study excitations generally. In :doc:`msdft2`, we have introduced how to use MSDFT to study core excitations. In :doc:`msdft3`, we have introduced how to use MSDFT to study double excitations. In this tutorial, we will introduce how to use MSDFT to study charge transfer excitations. We strongly recommend you to read these tutorials before reading this tutorial. Actually, MSDFT can be considered as TSO+NOSI. MSDFT is a framework that automaitcally performs TSO-DFT and NOSI for required excited states. Of course, you can also perform TSO-DFT and NOSI separately for special purposes. Example: Charge Transfer Excitations of [TCNE][CH\ :sub:`3`\ OCH=CH\ :sub:`2`] ---------------------------------------------------------------------------------- In this example, we will calculate the charge transfer excitations of the complex [TCNE][CH\ :sub:`3`\ OCH=CH\ :sub:`2`], whiere TCNE is tetracyanoethylene. First, let us calculate its ground state to see its orbitals. .. code-block:: :caption: msdft-4-gs.inp :linenos: basis cc-pVDZ # In product studies, you should use cc-pVTZ at least. end scf charge 0 spin2p1 1 type U end mol C -0.198382 0.676975 0.839378 C -1.315863 0.038931 0.405996 C -2.225365 0.670286 -0.508436 N -2.959552 1.174482 -1.247154 C -1.634267 -1.289698 0.850726 N -1.901673 -2.352749 1.221105 C 0.723624 0.024005 1.725081 N 1.483401 -0.507944 2.417545 C 0.131547 2.005023 0.406927 N 0.406859 3.079229 0.077027 C -0.018608 -1.360064 -2.097094 C 1.086004 -1.277514 -1.352322 H -0.458263 -2.333443 -2.305731 H -0.471114 -0.462748 -2.522620 H 1.590935 -2.157076 -0.934177 O 1.621801 -0.078506 -1.011442 C 2.981195 -0.123945 -0.597217 H 3.629890 -0.409943 -1.438672 H 3.239414 0.887438 -0.263120 H 3.112168 -0.831294 0.237956 end task energy b3lyp end After calculation, we can visualize ``msdft-4-gs.mwfn`` using `Qbics-MolStar `_. First drag ``msdft-4-gs.mwfn`` into explorer, and it will be automatically loaded. Right-click :guilabel:`msdft-4-gs.mwfn` and select :guilabel:`View Molecular Orbitals`, click orbital 48, 49, and 50: .. figure:: figs/a45.png .. figure:: figs/a46.png Obviously, MO48 is the HOMO mainly localized on CH\ :sub:`3`\ OCH=CH\ :sub:`2`, and MO49 and MO50 are the unoccupied π* and σ* orbitals mainly localized on TCNE, so excitation 48(HOMO) → 49 (LUMO) and 48(HOMO) → 50 (LUMO+1) are both close to but not perfect charge transfer. It seems that MO50 is more localized on TCNE but MO49 still has a few components on CH\ :sub:`3`\ OCH=CH\ :sub:`2`. We will use MSDFT to calculate these 2 excitations: .. code-block:: :caption: msdft-4-ct.inp :linenos: basis cc-pVDZ # In product studies, you should use cc-pVTZ at least. end scf charge 0 spin2p1 1 type U end mol C -0.198382 0.676975 0.839378 C -1.315863 0.038931 0.405996 C -2.225365 0.670286 -0.508436 N -2.959552 1.174482 -1.247154 C -1.634267 -1.289698 0.850726 N -1.901673 -2.352749 1.221105 C 0.723624 0.024005 1.725081 N 1.483401 -0.507944 2.417545 C 0.131547 2.005023 0.406927 N 0.406859 3.079229 0.077027 C -0.018608 -1.360064 -2.097094 C 1.086004 -1.277514 -1.352322 H -0.458263 -2.333443 -2.305731 H -0.471114 -0.462748 -2.522620 H 1.590935 -2.157076 -0.934177 O 1.621801 -0.078506 -1.011442 C 2.981195 -0.123945 -0.597217 H 3.629890 -0.409943 -1.438672 H 3.239414 0.887438 -0.263120 H 3.112168 -0.831294 0.237956 end msdft single_ex 48 : 49 50 end task msdft b3lyp end - ``msdft...end`` indicates the MSDFT calculation. Here, ``single_ex 48 : 49 50`` means we want to study the single excitation from orbitals 48 to orbitals 49 and 50. Explicitly, we want to study the following single excitations: - 48 → 49 - 48 → 50 In the output file ``msdft-4-ct.out``, we can find the following information: .. code-block:: :caption: msdft-4-ct.out :linenos: ---- NOSI Results ---- ====================== State NOSI Energies Excited Energy Osc. Str. DX DY DZ (Hartree) (eV) (a.u.) (a.u.) (a.u.) 0 -640.67704670 0.00000000 0.00000000 -1.24736 0.66912 0.44355 1 -640.56771187 2.97500080 0.00000000 -0.00000 -0.00000 0.00000 2 -640.55543283 3.30911339 0.17632767 -1.18048 0.77057 1.54122 3 -640.54664537 3.54822018 0.00000000 0.00000 -0.00000 -0.00000 4 -640.42108402 6.96474459 0.55299075 2.03421 0.55834 1.43209 ---- NOSI State Identification (Coefficients) ---- ================================================== State |0> = -0.999 |msdft-4-ct-gs.mwfn> State |1> = +0.710 |msdft-4-ct-48-to-49-se.mwfn> -0.710 |spin_flip_msdft-4-ct-48-to-49-se.mwfn> State |2> = -0.700 |msdft-4-ct-48-to-49-se.mwfn> -0.700 |spin_flip_msdft-4-ct-48-to-49-se.mwfn> State |3> = -0.708 |msdft-4-ct-48-to-50-se.mwfn> +0.708 |spin_flip_msdft-4-ct-48-to-50-se.mwfn> State |4> = -0.707 |msdft-4-ct-48-to-50-se.mwfn> -0.707 |spin_flip_msdft-4-ct-48-to-50-se.mwfn> State 2 and 4 are singlet charge transfer excitations. Example: Distance-Dependence of Charge Transfer Excitations ---------------------------------------------------------------------------------- Theoretical analysis reveals that the charge transfer excitation energies :math:`\omega` and the distance between the acceptor and donor :math:`d` has such a relationship: .. math:: \omega=A+\frac{B}{d} However, conventional TDDFT **cannot** repoduce this relationship because **TDDFT has systematic error for charge transfer excitations!** Let's see if MSDFT can reproduce this relationship for the above example. Generate Structures ^^^^^^^^^^^^^^^^^^^^^^^^^ To study the distance dependence, we will generate a series of [TCNE][CH\ :sub:`3`\ OCH=CH\ :sub:`2`] structures that the donor and acceptor departs. First, we save the structure in ``msdft-4-gs.inp`` in an XYZ file called ``r0.xyz``: .. code-block:: :caption: r0.xyz :linenos: 20 d = 3.21 Angstrom C -0.198382 0.676975 0.839378 C -1.315863 0.038931 0.405996 C -2.225365 0.670286 -0.508436 N -2.959552 1.174482 -1.247154 C -1.634267 -1.289698 0.850726 N -1.901673 -2.352749 1.221105 C 0.723624 0.024005 1.725081 N 1.483401 -0.507944 2.417545 C 0.131547 2.005023 0.406927 N 0.406859 3.079229 0.077027 C -0.018608 -1.360064 -2.097094 C 1.086004 -1.277514 -1.352322 H -0.458263 -2.333443 -2.305731 H -0.471114 -0.462748 -2.522620 H 1.590935 -2.157076 -0.934177 O 1.621801 -0.078506 -1.011442 C 2.981195 -0.123945 -0.597217 H 3.629890 -0.409943 -1.438672 H 3.239414 0.887438 -0.263120 H 3.112168 -0.831294 0.237956 Now, we write a python script to generate a series of structures: .. code-block:: python :caption: generate_structure.py :linenos: #!/usr/bin/python3 def read_xyz(path): with open(path, 'r') as f: lines = f.readlines() natoms = int(lines[0].strip()) title = lines[1].rstrip("\n") coords = [line.rstrip("\n") for line in lines[2:]] if len(coords) != natoms: raise ValueError("The number of atoms is inconsistent with the first line") return natoms, title, coords def write_xyz(path, natoms, title, coords): with open(path, 'w') as f: f.write(f"{natoms}\n") f.write(f"{title}\n") for line in coords: f.write(f"{line}\n") def translate_atoms(coords, vector, start_idx, end_idx): dx, dy, dz = vector for i in range(start_idx, end_idx + 1): if i >= len(coords): break parts = coords[i].split() if len(parts) < 4: raise ValueError(f"Line {i+3} format error: {coords[i]}") sym, x, y, z = parts[0], float(parts[1]), float(parts[2]), float(parts[3]) x += dx y += dy z += dz coords[i] = f"{sym:>2s} {x:15.8f} {y:15.8f} {z:15.8f}" return coords if __name__ == "__main__": dx = (-0.198382)-1.086004 dy = (0.676975)-(-1.277514) dz = (0.839378)-(-1.352322) r = (dx*dx+dy*dy+dz*dz)**0.5 ex = dx/r ey = dy/r ez = dz/r in_file = "r0.xyz" for i in range(1,9): natoms, title, coords = read_xyz(in_file) dr = 0.5*i d = 3.21+dr vector = (ex*dr, ey*dr, ez*dr) coords = translate_atoms(coords, vector, start_idx=0, end_idx=9) out_file = f"r{i:d}.xyz" title = f"d = {d:.2f} Angstrom" write_xyz(out_file, natoms, title, coords) print(f"Write {out_file}.") Run this python script by ``./generate_structure.py`` and we will generate 14 structures, with TCNE moving along axis atom 1-12 by 0.5 Angstorm every time. See below: .. figure:: figs/a47.jpg Calculate Energies ^^^^^^^^^^^^^^^^^^^^^^^^^ Now, prepare the input file: .. code-block:: :caption: msdft-4-rd.inp :linenos: basis cc-pVDZ # In product studies, you should use cc-pVTZ at least. end scf charge 0 spin2p1 1 type U end mol r0.xyz end msdft single_ex 48 : 49 50 end task msdft b3lyp end Now, we write a shell script to run the 15 calculations by replacing ``r0.xyz`` to ``r1.xyz``, etc.: .. code-block:: bash :caption: run.sh :linenos: inp=msdft-4-rd.inp for i in {0..9}; do sed "s/r0\.xyz/r${i}.xyz/g" "$inp" > "msdft-4-r${i}.inp" # Replace qbics-linux-cpu to the suitable one for your system qbics-linux-cpu "msdft-4-r${i}.inp" -n 24 > "msdft-4-r${i}.out" done When all calculations are completed, we can extract excitation energies to ``energies.txt``: .. code-block:: bash :caption: get_energies.sh :linenos: out="energies.txt" echo "Distance E1 E2" > ${out} for i in {0..9}; do d=$(echo "3.21 + $i * 0.5" | bc -l) read -r e1 e2 < <(awk 'NR==3{printf "%s ",$3} NR==5{printf "%s",$3}' "msdft-4-r${i}-spectrum.txt") echo "$d $e1 $e2" >> ${out} done Fit and Plot Data ^^^^^^^^^^^^^^^^^^^^^^^^^ Now, we write a python script to read the distances and excitation energies, fit the energies with respect to :math:`\omega=A+\frac{B}{d}`, and plot them: .. code-block:: python :caption: plot_fit.py :linenos: #!/usr/bin/python3 import matplotlib.pyplot as plt import numpy as np import matplotlib.gridspec as gridspec from scipy.stats import linregress # Fit data. def fit(a, b): res = linregress(a, b) return a*res.slope+res.intercept, res # Get data. fn = "energies.txt" data = np.loadtxt(fn, skiprows = 1).transpose() d = data[0] e1 = data[1] e2 = data[2] idx1 = 4 idx2 = 2 d1_fit = d[idx1:] d2_fit = d[idx2:] e1_fit, res1 = fit(1/d1_fit, e1[idx1:]) e2_fit, res2 = fit(1/d2_fit, e2[idx2:]) # Draw figures. plt.subplot(1, 2, 1) plt.plot(d, e1, "ro", label="Excitation: 48->49") plt.plot(d1_fit, e1_fit, "r-", label="Fit") plt.text(4, 3.35, f"$\omega = {res1.intercept:.3f}{res1.slope:.3f}/d$", fontsize=12, color='red') plt.text(4, 3.33, f"$R^2$ = {res1.rvalue**2:.3f}", fontsize=12, color='red') plt.xlabel("$d$ (Angstrom)", fontsize=12) plt.ylabel("$\omega$ (eV)", fontsize=12) plt.legend() plt.subplot(1, 2, 2) plt.plot(d, e2, "bo", label="Excitation: 48->50") plt.plot(d2_fit, e2_fit, "b-", label="Fit") plt.text(6, 7.025, f"$\omega = {res2.intercept:.3f}{res2.slope:.3f}/d$", fontsize=12, color='blue') plt.text(6, 7.015, f"$R^2$ = {res2.rvalue**2:.3f}", fontsize=12, color='blue') plt.xlabel("$d$ (Angstrom)", fontsize=12) plt.ylabel("$\omega$ (eV)", fontsize=12) plt.legend() # Finally, plot. plt.show() The plot is shown below: .. figure:: figs/a48.jpg Interstingly, at the begining, the excitation energies strongly fluctuate, but after 4 or 5 Angstrom, the excitation energy increases following the :math:`\omega=A+\frac{B}{d}` rule with very good R-value (0.969 and 0.994, respectively). It seems that 48 (HOMO) → 50 (LUMO+1) singlet excitation is more charge transfer-like!