Merge pull request #170 from sokolov-group/master

Updates to the ADC module user manual

Author: ssh2

source/user/adc.rst

@@ -1,152 +1,237 @@
 .. _user_adc:
 
 *****************************************
-Algebraic diagrammatic construction (ADC) 
+Algebraic diagrammatic construction (ADC)
 *****************************************
 
 *Modules*: :py:mod:`pyscf.adc`
 
 Introduction
 =============================
-The algebraic diagrammatic construction theory (ADC(n)) is a post-Hartree-Fock
-method used for computing correlated excited states of molecules.
-:cite:`Schirmer1983,Schirmer1998`
-The ADC methods involve a perturbative expansion of a propagator followed
-by truncation at a certain order :math:`n` that defines the ADC(n) approximation.
-Depending upon the property being investigated, propagators can be of different
+Algebraic diagrammatic construction theory (ADC) is a post-Hartree-Fock
+approach for computing excited electronic states.
+:cite:`Schirmer1983,Schirmer1998,Dreuw:2014p82,Banerjee:2023p3037`
+This is done by approximating a linear-response function
+(so-called propagator) to the :math:`n`-th order in perturbation theory,
+defining the hierarchy of ADC(n) methods.
+Depending on the property being investigated, propagators can be of different
 types. Some common examples include the polarization propagator for neutral
 electronic excitations, the one-particle Green's function for charged excitations,
 and the two-electron propagator for Auger electron spectroscopy.
-The different propagators lead to different variants of the ADC method.
-
-At present, the `adc` module in PySCF can be used to calculate the one-particle
-Green's function, that provides access to charged excitations,
-such as ionization potentials (IP-ADC) and
-electron affinities (EA-ADC) with different ADC(n) approximations
-(ADC(2), ADC(2)-X and ADC(3)). :cite:`Trofimov2005,Banerjee2019`
-The ADC(n) methods provide access to the
-excitation energies and their corresponding transition intensities via a
-'one-shot' calculation.
-
-A simple ADC(n) computation involves the calculation of the ground-state energy
-and wavefunction that correspond to those of the :math:`n`-th order
-Møller--Plesset perturbation theory (MPn), followed by the evaluation of the
-ADC(n) excitation energies (e.g., IP or EA) and the corresponding transition
-probabilities.
-
-An example of IP and EA computation using the default ADC method (ADC(2))
-is shown below::
+The different propagators lead to different variants of the ADC methods.
 
-        from pyscf import gto, scf, adc
-        mol = gto.M(atom='H 0 0 0; F 0 0 1', basis='ccpvdz')
-        mf = scf.RHF(mol).run()
-        myadc = adc.ADC(mf)
-        myadc.kernel_gs()
-        e_ip,v_ip,p_ip,x_ip,adc_es = myadc.ip_adc()
-        e_ea,v_ea,p_ea,x_ea,adc_ea = myadc.ea_adc()
+At present, the `adc` module in PySCF can be used to calculate neutral electronic
+excitations from the polarization propagator and charged excitations
+(ionization potentials, electron affinities) from the one-particle
+Green's function.
+For each property (EE, IP, and EA) three ADC approximations are available:
+strict second-order (ADC(2)), extended second-order (ADC(2)-X), and
+strict third-order (ADC(3)) methods. :cite:`Dreuw:2014p82,Banerjee:2023p3037`
+Each method provides access to the excitation energies and transition intensities
+via a 'one-shot' calculation.
 
-In the example shown above, the ground state calculation is performed by
-running the function ``kernel_gs()``, whereas the excited state calculations
-are carried out using the ``ip_adc()`` and ``ea_adc()`` functions.
+A conventional (single-reference) ADC(n) calculation involves computing the
+ground-state energy and wavefunction that correspond to those of the
+:math:`n`-th order Møller--Plesset perturbation theory (MPn), followed by evaluating
+the ADC(n) excitation energies and the corresponding transition probabilities.
 
-Alternatively, both the ground and excited state calculations can be run
-together using the ``kernel()`` function::
+Some examples of ADC calculations for the HF molecule are shown below::
 
         from pyscf import gto, scf, adc
         mol = gto.M(atom='H 0 0 0; F 0 0 1', basis='ccpvdz')
         mf = scf.RHF(mol).run()
-        myadc = adc.ADC(mf)
-        myadc.kernel()
 
-By default, the ``kernel()`` function performs a IP-ADC(2) calculation. One can specify the type of charged
-excitation and order of the desired ADC computation::
+        myadc = adc.ADC(mf)
+        myadc.verbose = 4
 
+        #IP-ADC(3) for 1 root
         myadc.method = "adc(3)"
+        myadc.method_type = "ip"
+        eip,vip,pip,xip = myadc.kernel()
+
+        #EA-ADC(2)-x for 1 root
+        myadc.method = "adc(2)-x"
         myadc.method_type = "ea"
-        myadc.kernel()
+        eea,vea,pea,xea = myadc.kernel()
 
-The ADC functions return the ``nroots`` lowest-energy eigenvalues. The
-default value of ``nroots`` is set to 1. More roots can be requested using::
+        #EE-ADC(2) for 3 roots
+        myadc.compute_properties = False
+        myadc.method = "adc(2)"
+        myadc.method_type = "ee"
+        eee,vee,pee,xee = myadc.kernel(nroots = 3)
 
-        myadc.kernel(nroots=3)
+Here, ``method`` specifies the level of ADC approximation,
+``method_type`` defines the type of excitation (EE, IP, and EA), and
+``nroots`` defines the number of excited states to be computed (default is 1).
 
 More examples can be found in
-:source:`examples/adc/01-closed_shell.py`,
-:source:`examples/adc/03-closed_shell_different_setup.py`.
+:source:`examples/adc`.
 
 
 Spin-restricted and spin-unrestricted calculations
 ==========================================================================
-The `adc` module can be used to perform calculations of IP's and EA's of closed- and
-open-shell molecules starting with the RHF and UHF reference
-wavefunctions, leading to the RADC(n) and UADC(n) methods, respectively.
-:cite:`Banerjee2021`
+The `adc` module can be used to compute EE's, IP's, and EA's
+for closed- and open-shell molecules starting with the restricted (RHF) or
+unrestricted (UHF) Hartree-Fock reference wavefunctions using the RADC or UADC
+modules.
+Additionally, open-shell ADC calculations can be carried out using the
+restricted open-shell (ROHF) reference wavefunction, which are supported by
+the UADC module.
 See :ref:`user_scf` to learn more about the different reference wavefunctions.
 
 Shown below is an example of the IP- and EA-UADC(2) calculations for the
 open-shell OH radical::
 
         from pyscf import gto, scf, adc
-        mol = gto.M(atom='O 0 0 0; H 0 0 1', basis='aug-cc-pvdz')
+        mol = gto.Mole()
+        mol.atom = [
+        ['O', ( 0., 0.    , 0.   )],
+        ['H', ( 0., 0.    , 0.969)],]
+        mol.basis = 'aug-cc-pvdz'
         mol.spin  = 1
+        mol.build()
+
         mf = scf.UHF(mol).run()
 
-        myadc = adc.ADC(mf)  # This is a UADC calculation
-        myadc.kernel_gs()
-        eip,vip,pip,xip,ip_es = myadc.ip_adc()
-        eea,vea,pea,xea,ea_es = myadc.ea_adc()
+        #IP-UADC calculation for 4 roots using the UHF reference
+        myadc = adc.ADC(mf)
+        myadc.verbose = 4
+
+        myadc.method = "adc(2)"
+        myadc.method_type = "ip"
+        eip,vip,pip,xip = myadc.kernel(nroots = 4)
+
+        mf = scf.ROHF(mol).run()
+
+        #EA-UADC calculation for 4 roots using the ROHF reference
+        myadc = adc.ADC(mf)
+        myadc.verbose = 4
+
+        myadc.method = "adc(2)"
+        myadc.method_type = "ea"
+        eea,vea,pea,xea = myadc.kernel(nroots = 4)
+
 
 More examples can be found in
-:source:`examples/adc/ 02-open_shell.py`,
-:source:`examples/adc/`04-open_shell_different_setup.py`.
+:source:`examples/adc`.
 
 
-Spectroscopic properties
+Excited-state properties
 =========================
-The `adc` module supports calculation of the spectroscopic factors, which provide
-information about probabilities of transitions in the photoelectron spectra. :cite:`Banerjee2021`
-Computation of spectroscopic factors is performed by default and can be switched
+The `adc` module supports computing transition and excited-state properties.
+For EE, oscillator strengths are available (currently, only in the UADC
+implementation). :cite:`Dreuw:2014p82`
+The IP- and EA-ADC methods provide spectroscopic factors, which measure the
+probabilities of charged excitations in photoelectron spectra. :cite:`Banerjee:2023p3037`
+Computation of transition properties is performed by default and can be switched
 off by setting ``compute_properties = False`` ::
 
         myadc.compute_properties = False
         myadc.method = "adc(3)"
         myadc.method_type = "ip"
         myadc.kernel(nroots = 3)
 
-After the ADC calculation is performed, the `adc` module can be used to compute
-the Dyson orbitals :cite:`Oana2007` corresponding to ionized and electron-attached states::
+Open-shell calculations using EE-UADC further support evaluating the spin
+square operator expectation values for the excited states (<S^2>). 
+The <S^2> values are not computed by default, they can be requested using
+the ``compute_spin_square = True`` flag. See a relevant example
+for more details: :source:`examples/adc/08-open_shell_spin_square.py`.
+
+For the IP- and EA-ADC methods, the `adc` module can be used to compute
+the Dyson orbitals :cite:`Oana2007` visualizing the wavefunction of
+ionized hole or attached electron for states with large
+spectroscopic factors (> 0.5). ::
 
         dyson_orb = myadc.compute_dyson_mo()
 
+Additionally, the ADC codes enable computing the one-particle reduced
+density matrices (1-RDMs) for the ground and excited electronic states.
+To obtain the ground-state 1-RDM, run the ``make_ref_rdm1()`` function
+after a successful ADC calculation::
+
+        myadc.kernel(nroots = 3)
+        rdm1_ref = myadc.make_ref_rdm1()
+
+The excited-state 1-RDMs can be calculated as follows::
+
+        rdm1_exc = myadc.make_rdm1()
+
+The 1-RDMs can be used to compute excited-state one-particle properties,
+such as dipole moments, see the following example:
+:source:`examples/adc/07-closed_shell_1RDMS.py`.
+The 1-RDM functionality is currently limited, see below.
 
-Analysis of spectroscopic properties
+Transition 1-RDMs between the ground and excited states are also available
+for all methods, except EE-RADC.
+They are returned as spectroscopic amplitudes at the end of a successful ADC
+calculation with ``compute_properties = True``, for example::
+
+        from pyscf import gto, scf, adc
+        mol = gto.Mole()
+        mol.atom = [
+        ['O', ( 0., 0.    , 0.   )],
+        ['H', ( 0., 0.    , 0.969)],]
+        mol.basis = 'aug-cc-pvdz'
+        mol.spin  = 1
+        mol.build()
+
+        mf = scf.UHF(mol).run()
+
+        # Transition 1-RDMs (alpha and beta spin) computed using EE-UADC(2)
+        myadc = adc.ADC(mf)
+        myadc.verbose = 4
+        myadc.method = "adc(2)"
+        myadc.method_type = "ee"
+        trdm_a, trdm_b = myadc.kernel(nroots = 4)[3]
+
+
+Analysis of excited-state properties
 =====================================
-The `adc` module allows to perform the analysis of the ADC(n) eigenvectors, that
-can be useful for characterizing the nature of electronic transitions. When
-``compute_properties`` is set to True, this analysis will also display the largest
-contributions to the spectroscopic factors. The analysis of the ADC(n) eigenvectors
+The `adc` module allows to analyze the ADC(n) eigenvectors
+for characterizing the nature of electronic transitions. When
+``compute_properties = True`` is set, this analysis will also display the largest
+contributions to the spectroscopic factors (IP, EA) or amplitudes (EE) for singly
+excited states. The analysis of ADC(n) eigenvectors
 and spectroscopic factors can be invoked using the ``analyze()`` function::
 
         myadc.kernel(nroots = 3)
         myadc.analyze()
 
 
+Core excitation energies and spectra
+=====================================
+The IP-ADC code supports calculating core ionization energies and
+X-ray photoelectron spectra (XPS) using the core-valence separation technique
+(CVS). To invoke CVS, specify the ``ncvs`` parameter 
+of IP-ADC object before running the ``kernel()`` function.
+The ``ncvs`` parameter should be set to an integer defining the index of the core
+orbital that is expected to be ionized first in the XPS spectrum.
+For example, probing the 1s orbital of C in CO can be done by setting
+``ncvs = 2`` since the C 1s orbitals are higher in energy than O 1s.
+
+
 Algorithms and job control
 ===========================
 
-The capabilities of the `adc` module at present are summarized in in the
+The current capabilities of the `adc` module are summarized in the
 following table:
 
-========== ========== ==================== ===============================
- Method     Reference  Spin-adaptation        Properties
----------- ---------- -------------------- -------------------------------
- ADC(2)     RHF, UHF    Yes                IP, EA, spectroscopic factors, Dyson orb
- ADC(2)-X   RHF, UHF    Yes                IP, EA, spectroscopic factors, Dyson orb
- ADC(3)     RHF, UHF    Yes                IP, EA, spectroscopic factors, Dyson orb
-========== ========== ==================== ===============================
-
-The ADC(n) calculations can be performed using different algorithms, depending on
-the available memory controlled by the ``max_memory`` keyword:
+============= ================ =====================================================
+ Method        Reference        Available properties
+------------- ---------------- -----------------------------------------------------
+ EE-ADC(2)     RHF, UHF, ROHF   EEs, osc. strengths (UADC), spin square (UADC)
+ EE-ADC(2)-X   RHF, UHF, ROHF   EEs, osc. strengths (UADC), spin square (UADC)
+ EE-ADC(3)     RHF, UHF, ROHF   EEs, osc. strengths (UADC), spin square (UADC)
+ IP-ADC(2)     RHF, UHF, ROHF   IPs, core IPs, spec. factors, Dyson orb
+ IP-ADC(2)-X   RHF, UHF, ROHF   IPs, core IPs, spec. factors, Dyson orb
+ IP-ADC(3)     RHF, UHF, ROHF   IPs, core IPs, spec. factors, Dyson orb
+ EA-ADC(2)     RHF, UHF, ROHF   EAs, spec. factors, Dyson orb
+ EA-ADC(2)-X   RHF, UHF, ROHF   EAs, spec. factors, Dyson orb
+ EA-ADC(3)     RHF, UHF, ROHF   EAs, spec. factors, Dyson orb
+============= ================ =====================================================
+
+The ADC(n) calculations are performed using one of the three algorithms for
+handling the two-electron integrals:
 
 * In-core
 
@@ -157,23 +242,59 @@ the available memory controlled by the ``max_memory`` keyword:
 
 * Out-of-core
 
-  Use of disk to store the expensive tensors.
-  This algorithm is invoked by setting ``max_memory`` to a small value.
+  Uses disk to store expensive tensors.
+  This algorithm is invoked when storing two-electron integrals requires
+  more memory than specified by the ``max_memory`` keyword (in MB).
   See :source:`examples/adc/05-outcore.py`
 
 
-* Density-fitted (DF) algorithm
+* Density fitting (DF)
+
+  The memory and disk usage can be greatly reduced by approximating the
+  two-electron integrals with density fitting. An example of the ADC(2)
+  calculation with density fitting is provided below::
+
+     from pyscf import gto, scf, adc, df
+     mol = gto.M(atom='H 0 0 0; F 0 0 1', basis='ccpvdz')
+
+     mf = scf.RHF(mol).density_fit('ccpvdz-jkfit').run()
+     myadc = adc.ADC(mf).density_fit('ccpvdz-ri')
+     eip,vip,pip,xip = myadc.kernel()
+
+  Density fitting introduces small errors in excitation energies
+  (~$10^{-3}$ eV), provided that the appropriate auxiliary basis
+  set is used. :cite:`Banerjee2021`
+  For calculations with more than 300 orbitals, using density fitting
+  is strongly recommended.
+
+  More examples can be found in: :source:`examples/adc/06-dfadc.py`.
+
+
+Current limitations
+====================
+
+Some limitations of current implementation are listed below:
+
+* The EE-RADC code does not support calculating oscillator strengths.
+  This property can be computed using the EE-UADC code (i.e., by using
+  the UHF reference) for a closed- or open-shell molecule.
+
+* The EE- and IP/EA-RADC codes compute only the states of lowest spin:
+  S = 0 (singlet) and S = 1/2 (doublet), respectively. Using the
+  corresponding UADC code allows to compute excitations with ΔS = 0,
+  ±1 for EE and ΔS = ±1/2, ±3/2 for IP and EA.
+
+* Computing spin square expectation values is currently only available for
+  EE-UADC.
 
- The memory and disk usage can be greatly reduced by approximating the
- two-electron integrals with density-fitting. A simple example of a
- DF-ADC(2) calculation is::
+* The reference and excited-state 1-RDMs are not implemented for EE-RADC.
+  Also, the reference 1-RDMs are not available for any UADC method.
 
-    from pyscf import gto, scf, adc, df
-    mol = gto.M(atom='H 0 0 0; F 0 0 1', basis='ccpvdz')
+* The EE-UADC(3) calculations of excited-state one-particle reduced density
+  matrices include correlation contributions up to EE-UADC(2)-X, i.e. the
+  third-order terms are missing from the singles-singles and singles-doubles
+  coupling sectors.
 
-    mf = scf.RHF(mol).density_fit('ccpvdz-jkfit').run()
-    myadc = adc.ADC(mf).density_fit('ccpvdz-ri')
-    eip,vip,pip,xip = myadc.kernel()
+* The EE-UADC(3) oscillator strengths do not include the contributions from
+  third-order amplitudes.
 
-More examples can be found in:
-:source:`examples/adc/06-dfadc.py`.