Natural orbitals (#55)

---------

Co-authored-by: Thomas Niehaus <[email protected]>

docs/_archives/recipes/linresp/natural-orbitals/dftb_in.hsd

@@ -0,0 +1,42 @@
+Geometry = GenFormat {
+  <<< "in.gen"
+}
+
+Hamiltonian = DFTB {
+
+  Scc = Yes
+  SccTolerance = 1.0E-10
+  MaxAngularMomentum {
+    N = "p"
+    C = "p"
+    H = "s"
+  }
+  SlaterKosterFiles = Type2FileNames {
+    Prefix = "../../slakos/download/mio/mio-1-1/"
+    Separator = "-"
+    Suffix = ".skf"
+  }
+}
+
+ExcitedState {
+  Casida {
+    NrOfExcitations = 10
+    StateOfInterest = 2	
+    Symmetry = Singlet
+    Diagonaliser = Stratmann {SubSpaceFactor = 30}
+    WriteEigenvectors = Yes		 
+  }
+}
+
+Options {
+  WriteDetailedXml = Yes
+}
+
+Analysis {
+  CalculateForces = Yes    
+  WriteEigenvectors = Yes
+}
+
+ParserOptions {
+  ParserVersion = 10
+}

docs/_archives/recipes/linresp/natural-orbitals/in.gen

@@ -0,0 +1,24 @@
+22 C
+N C H
+     1   1  -2.0382600000E+00   2.2050800000E+00   2.5000000000E-04
+     2   2  -1.4355300000E+00   1.2277400000E+00   1.3000000000E-04
+     3   2  -6.8509500000E-01   0.0000000000E+00   0.0000000000E+00
+     4   2  -1.4354800000E+00  -1.2276900000E+00   8.0000000000E-05
+     5   2   6.8509500000E-01   0.0000000000E+00   0.0000000000E+00
+     6   1  -2.0380300000E+00  -2.2051400000E+00   2.5000000000E-04
+     7   2   1.4355200000E+00  -1.2276900000E+00  -1.1000000000E-04
+     8   2   1.4354800000E+00   1.2277100000E+00  -9.0000000000E-05
+     9   1   2.0382900000E+00  -2.2050000000E+00  -2.3000000000E-04
+    10   1   2.0380500000E+00   2.2051400000E+00  -2.4000000000E-04
+    11   2   0.0000000000E+00   1.3939400000E+00   5.0000000000E+00
+    12   2   1.2071800000E+00   6.9697000000E-01   5.0000000000E+00
+    13   2   1.2071800000E+00  -6.9697000000E-01   5.0000000000E+00
+    14   2   0.0000000000E+00  -1.3939400000E+00   5.0000000000E+00
+    15   2  -1.2071800000E+00  -6.9697000000E-01   5.0000000000E+00
+    16   2  -1.2071800000E+00   6.9697000000E-01   5.0000000000E+00
+    17   3  -2.1463200000E+00   1.2391800000E+00   5.0000000000E+00
+    18   3  -2.1463200000E+00  -1.2391800000E+00   5.0000000000E+00
+    19   3   0.0000000000E+00  -2.4783600000E+00   5.0000000000E+00
+    20   3   2.1463200000E+00  -1.2391800000E+00   5.0000000000E+00
+    21   3   2.1463200000E+00   1.2391800000E+00   5.0000000000E+00
+    22   3   0.0000000000E+00   2.4783600000E+00   5.0000000000E+00

docs/_archives/recipes/linresp/natural-orbitals/waveplot_in.hsd

@@ -0,0 +1,18 @@
+Options {
+  RealComponent = Yes                # Plot real component of the wavefunction
+  PlottedSpins = 1 -1
+  PlottedLevels = 1 -1               # Levels to plot
+  PlottedRegion =  OptimalCuboid {}  # Region to plot
+  NrOfPoints = 50 50 50              # Number of grid points in each direction
+  NrOfCachedGrids = -1               # Nr of cached grids (speeds up things)
+  Verbose = Yes                      # Wanna see a lot of messages?
+}
+
+DetailedXml = "detailed.xml"         # File containing the detailed xml output
+                                     # of DFTB+
+EigenvecBin = "excitedOrbs.bin"      # File cointaining the binary eigenvecs
+
+Basis {
+  Resolution = 0.01
+  <<+ "../../slakos/wfc/wfc.mio-1-1.hsd"
+}

docs/_static/references.bib

@@ -173,3 +173,25 @@ @article{salthammer2022
 year = {2022},
 doi = {10.1021/acs.est.1c06935}
 }
+
+@article{luzanov1976application,
+  title={Application of transition density matrix for analysis of excited states},
+  author={Luzanov, AV and Sukhorukov, AA and Umanskii, VE},
+  journal={Theoretical and Experimental Chemistry},
+  volume={10},
+  pages={354--361},
+  year={1976},
+  publisher={Springer}
+}
+
+@article{martin2003natural,
+  title={Natural transition orbitals},
+  author={Martin, Richard L},
+  journal={The Journal of chemical physics},
+  volume={118},
+  number={11},
+  pages={4775--4777},
+  year={2003},
+  publisher={American Institute of Physics}
+}
+

docs/linresp/diatomic.rst

@@ -24,7 +24,7 @@ excited states with both singlet and triplet multiplicities. The
 solution of the Casida equation is performed after a single-point
 (static) DFTB ground-state calculation for the optimised geometry. All
 this is done in a single DFTB run (one single input file). The
-*dftb_in.hsd* input file should look like this::
+`dftb_in.hsd` input file should look like this::
 
   Geometry = GenFormat {
       <<< "geo_in.gen"
@@ -58,8 +58,8 @@ this is done in a single DFTB run (one single input file). The
   }
 
 The input shows a standard spin-unpolarised static DFTB calculation,
-except for a new block denoted as *ExcitedState* with an embedded
-subblock called *Casida*. Once the SCC loop is converged, this block
+except for a new block denoted as ``ExcitedState`` with an embedded
+subblock called ``Casida``. Once the SCC loop is converged, this block
 instructs DFTB+ to build the response matrix according to Casida,
 using the just-obtained DFTB orbitals and energies, and diagonalise
 it.
@@ -72,9 +72,9 @@ triplet excitations. That is, if we would be only interested in
 singlet transitions, we could just omit the spin constants in the
 input file.
 
-Let us now take a closer look at the *Casida* block. The keyword
-*NrOfExcitations* specifies how many transitions per spin symmetry, or
-multiplicity, we want to compute. In *Symmetry* we specify the
+Let us now take a closer look at the ``Casida`` block. The keyword
+``NrOfExcitations`` specifies how many transitions per spin symmetry, or
+multiplicity, we want to compute. In ``Symmetry`` we specify the
 multiplicity of the transition (either singlet, triplet or both). The
 multiplicity of the transition is the difference between the
 multiplicities of the excited state and the ground state. In our
@@ -85,7 +85,7 @@ means the first 10 singlet-to-singlet and 10 singlet-to-triplet
 transitions.
 
 Once the calculation is finished (it takes a fraction of a second),
-the output file *EXC.DAT* contains the excitation energies and
+the output file `EXC.DAT` contains the excitation energies and
 oscillator strengths, as well as other valuable information. It should
 look like this::
 
@@ -117,20 +117,20 @@ look like this::
 
 
 The triplet transitions are listed first, followed by the singlet
-ones. They can be identified by the letter *T* or *S* in the last
+ones. They can be identified by the letter ``T`` or ``S`` in the last
 column.
 
-The first column *w [eV]* is the excited state energy we are looking
-for, the second one *Osc.Str.* lists the corresponding oscillator
-strength. The column *Transition* reports the indices of the dominant
+The first column ``w [eV]`` is the excited state energy we are looking
+for, the second one ``Osc.Str.`` lists the corresponding oscillator
+strength. The column ``Transition`` reports the indices of the dominant
 molecular orbitals involved in the electronic transition. In our
 example, the singlet state at 12.75 eV features a transition from the
 occupied Kohn-Sham orbital 3 (HOMO-2) to the virtual orbital 6 (the
-LUMO). The next column *Weight* indicates the weight of the
+LUMO). The next column ``Weight`` indicates the weight of the
 corresponding singly excited determinant in the CIS expansion of the
 excited state. Values close to one indicate that the excited state is
 well described by a single electronic excitation, while small values
-speak for a collective excitation. Column *KS [eV]* provides the
+speak for a collective excitation. Column ``KS [eV]`` provides the
 Kohn-Sham transition energy difference :math:`\omega_{ia\sigma} =
 \epsilon_{a\sigma} - \epsilon_{i\sigma}` (see above).
 
@@ -141,7 +141,7 @@ Oxygen molecule
 
 For the |O2| molecule, we will consider its triplet ground state. This
 is specified in the input file through the
-*Hamiltonian/SpinPolarisation* block::
+``Hamiltonian/SpinPolarisation`` block::
 
   SpinPolarisation = Colinear {
       UnpairedElectrons = 2
@@ -162,7 +162,7 @@ our eigenvalue problem into two independent singlet and triplet
 equations, so we have to build and diagonalise the entire response
 matrix in this case. But, how do we know the spin multiplicities of
 the computed transitions? We get this information from the last column
-of the *EXC.DAT* file::
+of the `EXC.DAT` file::
 
   w [eV]       Osc.Str.         Transition         Weight      KS [eV]    D<S*S>
 
@@ -227,7 +227,7 @@ In this case, the first 10 excitations are::
   8.636        0.00000000      3   ->     7        0.657       8.636    -0.000
  11.652        0.49971991      3   ->     6        0.600       8.636    -0.597
 
-Let us pay attention to the last column of the *EXC.DAT*
+Let us pay attention to the last column of the `EXC.DAT`
 file. Contrary to the previous case, here we obtain large non-zero
 :math:`\Delta S^2` values. When :math:`\Delta S^2 = 0`, we are in the
 presence of a doublet-to-doublet transition. Likewise, if
@@ -236,7 +236,7 @@ state. Otherwise, we have some extent of spin contamination in our
 obtained transitions. The last column should help us determine which
 excitations are to be trusted. We can set an arbitrary spin
 contamination threshold to establish which transitions we will
-consider leading to a physical excited state.  In our *NO-TiO2*
+consider leading to a physical excited state.  In our ``NO-TiO2``
 recipe, we will compute the absorption spectrum of a system where
 transitions with a spin contamination beyond an imposed threshold are
 excluded.

docs/linresp/index.rst

@@ -43,3 +43,4 @@ of TD-DFTB will be exploited.
    no-titania.rst
    rangeseparated.rst
    relax.rst
+   naturalorbitals.rst

docs/linresp/macromolecule.rst

@@ -25,12 +25,12 @@ computed as usual::
 
 After obtaining the energy and oscillator strength pairs (takes very
 few minutes), we can use a simple script to plot the absorption
-spectrum. The script (*spectrum.py*) is provided in the input folder
+spectrum. The script (`spectrum.py`) is provided in the input folder
 and requires Python 3.  Simply run::
 
   python spectrum.py > abs.dat
 
-and plot the content of *abs.dat* with your graphing tool of
+and plot the content of `abs.dat` with your graphing tool of
 preference.  The script employs a Lorentzian function for the spectral
 broadening, with a full width at half-maximum parameter of
 :math:`\Gamma` = 0.2 eV. Specifically, what is plotted is the molar
@@ -75,8 +75,8 @@ single-particle transitions and roughly equal to 1 Hartree). As for
 the oscillator strength window, we will use a cutoff of
 :math:`10^{-3}`. These values ensure a negligible loss of accuracy in
 the optical spectra, while reducing the computation time by roughly
-half. Those are specified in the *Casida* block using the
-*EnergyWindow* and *OscillatorWindow* keywords::
+half. Those are specified in the ``Casida`` block using the
+``EnergyWindow`` and ``OscillatorWindow`` keywords::
 
   ExcitedState {
       Casida {

docs/linresp/naturalorbitals.rst

@@ -0,0 +1,72 @@
+.. highlight:: none
+
+***********************************************
+Visualization of natural orbitals with waveplot
+***********************************************
+
+[Input: `recipes/linresp/natural-orbitals`]
+
+The majority of excited states are not well described by the
+transition of an electron from a single Kohn-Sham orbital to another
+single Kohn-Sham orbital. In other words, the excited state is in
+general a superposition of many determinants. In such a situation, the
+excited state weight in the file EXC.DAT will be smaller than
+one. Natural transition orbitals (NTO) allow for a compact
+representation of the excitation in such a scenario :cite:`luzanov1976application`. They are obtained
+from a singular value decomposition of the transition density matrix
+and deliver `hole` and `electron` orbitals for a given excited
+state. We discuss their visualization in the context of the
+Benzene-Tetracyanoethylene dimer we looked at :ref:`earlier <benz_dimer>`. Strictly
+speaking, natural orbitals are not really required for this system,
+since for example the lowest excited state is well described by a HOMO
+to LUMO transition, without contributions from other single-particle
+transitions. We still discuss the NTO feature for this example, because
+it allows us to compare the results with the earlier calculations.
+
+We first run DFTB+ with the following modified input:
+
+.. literalinclude:: ../_archives/recipes/linresp/natural-orbitals/dftb_in.hsd
+   :lines: 21-38
+
+Here ``CalculateForces`` requests to compute excited forces for state
+``StateOfInterest``. During this calculation also the NTO for that state
+are created and written to a file (always called `excitedOrbs.bin`), if
+``WriteEigenvectors`` in the ``Casida`` block is set to ``Yes``. In the
+``Options`` block, the keyword ``WriteDetailedXml`` is required for the
+plot. In the ``Analysis`` block, ``WriteEigenvectors`` advises the code to
+print out the ground state molecular orbitals also.
+
+Let us look at an excerpt of the generated file `detailed.xml`::
+
+ <excitedoccupations>
+  <spin1>
+   <k1>
+    -1.00000341372109 -5.201183337079720E-002 -5.170515789354190E-002 .....
+    .....................
+    5.170515780471528E-002 5.201183337096685E-002 1.00638172624593
+
+   </k1>
+  </spin1>
+ </excitedoccupations>
+
+We see the state occupations :math:`n_{i}` of the NTO. Negative values
+indicate hole orbitals and positive values electron orbitals. Most
+excited states (unless you have degeneracy) have only one important
+(:math:`n_{i}\approx -1.0`) hole and one important
+(:math:`n_{i}\approx 1.0`) electron NTO, even though the CI expansion
+might include a large number of single-particle transitions.
+
+As seen in the section  :ref:`here <sec-basics-waveplot>`, we can now visualize these orbitals using the waveplot code. The input file `waveplot_in.hsd` has this form:
+
+.. literalinclude:: ../_archives/recipes/linresp/natural-orbitals/waveplot_in.hsd
+   :emphasize-lines: 4,13
+
+Important is here the line ``EigenvecBin`` which reads the created file
+with natural transition orbitals instead of the ground state MO. The
+keyword ``PlottedLevels`` in this example requests the first and last
+orbital to be plotted (note that the curly bracket contains the
+indices of the NO, not the occupation). These are exactly the natural
+transition orbitals with highest occupation. If you do the plot, you
+will realize that the hole orbital is essentially Kohn-Sham state 37
+(the HOMO), while the electron orbital is Kohn-Sham state 38 (the
+LUMO).

docs/linresp/no-titania.rst

@@ -29,7 +29,7 @@ the UV-A region, so that we can observe both near band-gap excitations
 of TiO\ :sub:`2` and ligand-to-substrate charge-transfer (CT)
 transitions in the visible range. The challenge here is in the high
 number of excitations that need to be solved to cover the desired
-energy region. Specifically, we have to set *NrOfExcitations* to 1200
+energy region. Specifically, we have to set ``NrOfExcitations`` to 1200
 to cover a region of up to 4 eV::
 
   ExcitedState {
@@ -40,18 +40,18 @@ to cover a region of up to 4 eV::
       }
   }
 
-we have set *EnergyWindow = 0.01* to save computational time. We can
+we have set ``EnergyWindow = 0.01`` to save computational time. We can
 choose such a small energy cutoff as titania band-to-band transitions
 have a rather low collective character. The visible range of the
 spectrum, which is of most interest, is negligibly affected by this
 truncation, as the CT absorption bands are away from the energy window
 frontier. This constraint speeds up the calculation roughly 7 times!
 
-After obtaining the *EXC.DAT* output file (we provide it in the recipe
+After obtaining the `EXC.DAT` output file (we provide it in the recipe
 folder as the calculation will take some time to complete), we can run
-the *spectrum.py* script to obtain the absorption spectrum. Only those
+the `spectrum.py` script to obtain the absorption spectrum. Only those
 transition with a spin contamination smaller than 0.5 will be
-considered. This value is hard-coded in *spectrum.py*, and the reader
+considered. This value is hard-coded in `spectrum.py`, and the reader
 can modify it at his/her convenience. The absorption spectrum should
 finally look like this: