Effective Electronic Hamiltonian for Quantum Subsystem in Hybrid QM/MM

Methods as Derived from APSLG Description of Electronic Structure of Classical Part of Molecular System

A.M. Tokmachev, A.L. Tchougréeff and I.A. Misurkin Karpov Institute of Physical Chemistry Vorontsovo pole 10, Moscow 103064 RUSSIA

September 23, 1999

Abstract

The general formulae representing separation of electronic variables of quantum (reactive) subsystem from those describing electrons in the classical (chemically inert) part of molecular system are specified for the case when the electronic structure of the latter is described by a semi-empirical method based on the trial wave function having the form of antisymmetrized product of strictly localized geminals (APSLG) which leads to a local description of molecular electronic structure in terms of bond functions and lone pair functions. This allowed us to give explicit form of the effective electronic Hamiltonian for the quantum subsystem and by this also to sequentially derive the explicit form of the QM/MM junction between the quantum and classical subsystems. The latter turned out to be a sum of the contributions from different chemical bonds and lone pairs residing in the classical part of the system. Numerical estimates for the effect of the renormalization of the Coulomb interaction of π-electrons due to presence of σ-bonds are performed according to the derived formulae.

1. Introduction

At present research in chemistry and in related areas of science reached the state that require constructing potential energy surfaces (PES) of large systems. This problem can be encountered in the context of chemical reactions of biomolecules, enzymatic reactions, surface reactions and reactions in condensed media. Applying methods of quantum mechanics (QM) to construct PES’s in each point of the nuclear configuration space faces the problem of O(N^{m})-scaling of the QM methods (N is the number of one-electron states involved in the calculation). Practically the exponent m may reach values of 5÷7 for high quality modern QM methods necessary to describe chemical reactivity (bond cleavage and bond formation), which restricts their applicability to molecular systems of rather small size [1]. However, the detailed QM description is necessary for the electronic structure of the reactive site only. The contribution of the rest of molecular system (i.e. of its chemically inert part) to the PES of the system can be calculated by molecular mechanics (MM), which has to reproduce only general features of this part of the molecular system. Thereby the hybrid quantum-mechanical/molecular mechanical (QM/MM) computational schemes become very popular (Refs. [2–7]). This approach significantly reduces the computational costs of PES construction for large systems because only a small part of the latter is considered on the computationally expensive QM level.

There exist several QM/MM schemes implemented in a number of computation packages (see, for example, Refs. [2–7]). The diversity of such approaches is due to variety of both QM and MM methods combined and of the functional form of the junction between them. For example, in Ref. [2] the intersystem Coulomb and exchange integrals are represented as linear combinations of exponential functions with subsequent parameterization; in Ref. [3] ”junction dummy atoms” are introduced; in Ref. [8] the interactions between subsystems are represented also by a sum of electrostatic and van der Waals interactions with adjusted parameters, etc. The common feature of all the mentioned approaches as well as of many others is that the form of the junction parameterized in each of them results from certain ad hoc postulate rather than from a sequential derivation. On the other hand, in Ref. [9] it was proposed to construct a consistent form of the QM/MM junction with use of an explicit procedure of separation of electronic variables. The latter is performed in Ref. [10] by using the Löwdin partition technique [11] and the group function (GF) formalism [12]. In Refs. [9,10] the expression for the effective Hamiltonian for the quantum (reactive) part of a molecular system has been obtained and the form of the QM/MM junction has been represented as a sum of operator averages over the implicit wave function of the classical (inert) subsystem, which is assumed to describe the electronic ground state of the latter. The present work is devoted to derivation of the explicit form of the effective Hamiltonian for the quantum system and, therefore, of the PES of the combined system and of the specific form of the QM/MM junction for a special local form of the wave function of the inert subsystem which has been constructed in Ref. [14] to ensure a ready transition to a description of the MM type which will be described elsewhere [16].

2. Effective Hamiltonian for quantum system

Now we briefly review the main notations and results of Ref. [10], where the general expression for the effective Hamiltonian has been proposed. We denote two subsystems of the whole molecular system by indices R (reactive) and M (inert), where the R-subsystem is considered as a quantum one, whereas the Msubsystem must be finally treated on the classical (MM) level of approximation. The Hamiltonian for the whole system is a sum of Hamiltonians for the subsystems and of the operators for the interactions between the subsystems:

H = H^{R}(q) + H^{M}(q) + V^{c}(q) + V^{r}(q),  (1)

where for the sake of simplicity only the Coulomb V^{c} and the resonance V^{r} (electron transfer) interactions are considered. Further, the Hamiltonian for the M-subsystem is subdivided into the Hamiltonian for the free (without reagents) M-subsystem H_{0}^{M} and the operator describing attraction of electrons of the Msubsystem to the cores of the R-subsystem V^{R}. Analogous subdividing is performed for the R-subsystem. To justify usage of different levels of approximation to different parts of the whole system (specifically, of the MM-like scheme for the M-subsystem) the wave function for the whole system is represented by the antisymmetrized product of electronic wave function for the R-subsystem and that of the ground state for the free M-subsystem (i.e. of the ground state function of H_{0}^{M}):

Ψ_{k} = Φ_{k}^{R} ∧ Φ_{00}^{M}.  (2)

The exact wave function of any electronic state of the whole system can be recast in the form:

Ψ_{k} = ∑_{n_{M} n_{R}} ∑_{i_{M} i_{R}} C_{i_{M} i_{R}}^{k}(n_{M}n_{R})Φ_{i_{M} 0}^{M}(n_{M}) ∧ Φ_{i_{R}}^{R}(n_{R}), (n_{M} + n_{R} = N_{e}) (3)

The transition from the wave function of the general form Eq.(3) to the necessary form of Eq. (2) is made by performing two sequential Löwdin projection procedures: the first one to the subspace of the states with fixed number of electrons in the subsystems (projection operator P and its complementary projection operator Q = 1 − P) and the second one - to the states of the type Eq. 2, i. e. containing the ground state wave function of the free M-subsystem as the multiplier (projection operator P and its complementary projection operator Q = 1 − P). After the first projection we obtain

\begin{matrix}H_{eff}(q, E) = P H^{R}(q)P + P H^{M}(q)P + P V^{c}(q)P + \\ +P V^{rr}(q, E)P + \frac{e^{2}}{2} ∑_{\begin{matrix}A≠B \\ A∈R,B∈M \end{matrix}} Z_{A}^{R}Z_{B}^{M}R_{AB}^{−1}, \end{matrix}  (4)

where

V^{rr}(q, E) = V^{r}(q)Q(E − QHQ)^{−1}QV^{r}(q) = V^{r}(q)QR(q, E)QV^{r}(q). (5)

The second projection and subsequent averaging over the ground state of the M-subsystem gives the effective Hamiltonian for the R-subsystem:

\begin{matrix}H_{eff}^{R}(q, E, ω) = H_{0}^{R}(q) + δV^{M} + ⟨⟨P V^{rr}P ⟩⟩_{M} + ⟨Φ_{00}^{M} | P V^{R}P | Φ_{00}^{M}⟩+ \\ +⟨Φ_{00}^{M} | P W(q, E)P QR(ω)QP W(q, E)P | Φ_{00}^{M}⟩ + \frac{e^{2}}{2} ∑_{\begin{matrix}A≠B \\ A∈R,B∈M \end{matrix}} Z_{A}^{R}Z_{B}^{M}R_{AB}^{−1}, \end{matrix} (6)

where

δV^{M} = V^{M} + ⟨⟨P V^{c}P ⟩⟩_{M} ≈ 0  (7)

and

P WP = P V^{c}P + P V^{rr}P + P V^{R}P.  (8)

The above general form of the effective Hamiltonian was obtained in Ref. [10]. In the present paper we perform the averaging assuming that the wave function of the M-subsystem has a specific form, which gives the local description of the latter.

The form of the wave function of the M-subsystem to be used to perform the averaging has to allow to represent the renormalization of Hamiltonian for the quantum subsystem as a sum of contributions from one or more chemical bonds or lone pairs. This is done to maintain consistency with the adopted MM type of description for the M-subsystem. For this purpose we must use a quantum-chemical method, the energy of which can be presented in a MM-like form. The MM scheme assumes the transferability of the functions representing geometry dependence of different contributions to the molecular energy. Such a transferability is shown to be achieved (Ref. [14]) for the trial wave function in the form of the antisymmetrized product of strictly localized orbitals (APSLG) [15]. The wave function of this method is constructed from two-electron functions (geminals) assigned to chemical bonds and lone pairs:

|Φ⟩ = ∏_{k} g_{k}^{+} | 0⟩,  (9)

where

\begin{matrix}g_{k}^{+} = u_{k}r_{kα}^{+}r_{kβ}^{+} + v_{k}l_{kα}^{+}l_{kβ}^{+} + w_{k}(r_{kα}^{+}l_{kβ}^{+} + l_{kα}^{+}r_{kβ}^{+}), \\ (u_{k}^{2} + 2w_{k}^{2} + v_{k}^{2} = 1) \end{matrix}  (10)

is the electronic pair creation operator for the k-th geminal. Each of the orbitals r_{k} and l_{k} assigned to the k-th chemical bond is a linear combination of the AO’s centered on one atom only, i.e. a hybrid orbital (HO). The unitary matrices of transition from the AO basis to the HO basis and the geminal expansion coefficients u_{k}, v_{k}, and w_{k} are determined variationally in Ref. [14] for a wide range of organic molecules. For the purposes of the present paper it is important to note that the energy of molecule in the APSLG approximation can be represented by a sum of interbond and intrabond (we use the term bond for usual chemical bonds and lone pairs) contributions and that the parameters of these contributions are well transferable. The derivation of the MM description from the QM APSLG method is rather complex and will be published elsewhere [16].

Let us consider the averages renormalizing the Hamiltonian for the R-subsystem, carrying out the summation over the spin projections. The operators of the Coulomb and the resonance interactions between the subsystems can be written as

\begin{matrix}V^{c}(q) = ∑_{\begin{matrix}pp′∈R \\ mm′∈M \end{matrix}}(pp′ ∥ mm′)p^{+}m^{+}m′p′, \\ V^{r}(q) = ∑_{\begin{matrix}p∈R \\ m∈M \end{matrix}} v_{pm}(q)(p^{+}m + m^{+}p), \end{matrix}  (11)

where

(pp′ ∥ mm′) = (pp′ | mm′) − (pm′ | mp′).  (12)

and the indices pp′ and mm′ refer to the one electron states in the R- and Msubsystems, respectively. In the latter case the one-electron states can be taken as the HO’s r_{k} and l_{k} in the M-subsystem. The averaging of the operator of Coulomb interaction between the subsystems yields:

\begin{matrix}⟨⟨P V^{c}P ⟩⟩_{M} = ⟨Φ_{00}^{M}| P V^{c}P |Φ_{00}^{M}⟩ = \\ = ∑_{pp′∈R} p^{+}p′ [ ∑_{mm′∈M} (pp′| |mm′) ⟨⟨m^{+}m′⟩⟩_{M} ] . \end{matrix}  (13)

Since we assume that the M-subsystem is described by the wave function of the APSLG type, all the orbitals in the M-subsystem are either right (r) or left (l) orbitals of geminals. In the APSLG approximation the averages ⟨⟨m^{+}m′⟩⟩_{M} do not vanish only for the spin-orbitals m and m′ belonging to the same geminal and can be written as

\begin{matrix}⟨⟨r_{kσ}^{+}r_{kσ}⟩⟩_{M} = ⟨0| g_{k}r_{kσ}^{+}r_{kσ}g_{k}^{+} |0⟩ = P_{k}^{rr} = u_{k}^{2} + w_{k}^{2}, \\ ⟨⟨l_{kσ}^{+}l_{kσ}⟩⟩_{M} = ⟨0| g_{k}l_{kσ}^{+}l_{kσ}g_{k}^{+} |0⟩ = P_{k}^{ll} = v_{k}^{2} + w_{k}^{2}, \\ ⟨⟨r_{kσ}^{+}l_{kσ}⟩⟩_{M} = ⟨⟨l_{kσ}^{+}r_{kσ}⟩⟩_{M} = P_{k}^{rl} = (u_{k} + v_{k})w_{k}. \end{matrix}  (14)

It is also convenient to introduce the reduced Coulomb integrals

Y_{pp′mm′} = 2 (pp′ | mm′) − (pm′ | mp′) .  (15)

Therefore, we can write the average ⟨⟨P V^{c}P ⟩⟩_{M} as

\begin{matrix}⟨⟨P V^{c}P ⟩⟩_{M} = ∑_{pp′∈R} ∑_{σ} p_{σ}^{+}p′_{σ}× \\ × ∑_{k∈M} [Y_{pp′r_{k}r_{k}}P_{k}^{rr} + Y_{pp′l_{k}l_{k}}P_{k}^{ll} + (Y_{pp′l_{k}r_{k}} + Y_{pp′r_{k}l_{k}})P_{k}^{rl}] . \end{matrix}  (16)

The ZDO approximation assures that p = p′ and m = m′. We denote p ∈ A and m ∈ B. In the case of A ≠ B we obtain the contributions of the type

∑_{σ} p_{σ}^{+}p_{σ}γ_{AB} [ ∑_{\begin{matrix}m_{k}∈B \\ τ \end{matrix}} ⟨⟨m_{kτ}^{+}m_{kτ} ⟩⟩_{M} ] = 2 ∑_{σ} p_{σ}^{+}p_{σ}γ_{AB} ∑_{m_{k}∈B} P_{k}^{mm}.  (17)

The next contribution to δV^{M} is

\begin{matrix}V^{M} = −e^{2} ∑_{B′} \frac{Z_{B′}^{M}}{|r−R_{B′}|} = \\ = − ∑_{pp′∈R} ∑_{B∈M} V_{B}^{pp′}Z_{B}^{M} ∑_{σ} p_{σ}^{+}p′_{σ}. \\ V_{B}^{pp′} = −e^{2} ∫ d^{3}r\frac{ψ_{p}^{∗}(r)ψ_{p′}(r)}{|r−R_{B′}|} \end{matrix}  (18)

Taking the sum of Eq. (17) and (18) and using the ZDO approximation we obtain that the contribution to δV^{M} from the interactions of electrons belonging to different atoms equals

∑_{A∈R} ∑_{p∈A∩R} ∑_{σ} p_{σ}^{+}p_{σ} ∑_{\begin{matrix}B∈M \\ B≠A \end{matrix}} γ_{AB}Q_{B}^{M},  (19)

where Q_{B}^{M} = 2 ∑_{m_{k}∈B} P_{k}^{mm} − Z_{B}^{M} is the effective charge of the atom B. The con-

tribution from the interactions of electrons located on the orbitals belonging to different systems but centered on the same frontier atoms equals to

∑_{pp′∈A} ∑_{σ} p_{σ}^{+}p′_{σ} ∑_{\begin{matrix}m_{k}∈A \\ k∈M \end{matrix}} P_{k}^{mm}Y_{pp′mm}.  (20)

The next contribution to the effective Hamiltonian for the R-subsystem is due to the intersubsystem electron transfers

\begin{matrix}P V^{rr}(q, E)P = ∑_{\begin{matrix}pp′∈R \\ mm′∈M \end{matrix}} v_{pm}(q)v_{p′m′}(q)× \\ × [(p^{+}mR(E)m^{′+}p′) + (m^{+}pR(E)r^{′+}p′)] . \end{matrix}  (21)

The resolvent can be presented as

R(E) = ∑_{i∈ImQ} \frac{|i⟩ ⟨i|}{E − E_{i}} ,  (22)

where |i⟩ are the states with one electron transferred from the M-subsystem to the R-subsystem and vice versa. We assume that every state |i⟩ is an antisymmetrized product of ionized states |μ⟩ and |ρ⟩ of the M- and R-subsystems respectively. Moreover, we construct the states |μ⟩ by removing or adding an electron from or to the Dirac orbitals of the M-subsystem. Therefore, we can assume the energy differences in this equation to be expressed through the ionization potentials (IP) and electron affinities (EA):

E_{i} − E = \{ \begin{matrix}I_{μ} − A_{ρ} − g_{μρ} \\ I_{ρ} − A_{μ} − g_{ρμ} \end{matrix} .  (23)

Now we specify the approximate form of the states |i⟩ in the resolvent. First, we notice that the correlations and bonding can be accounted to the same extent as they are in the geminal Eq. (10) if one employs bonding (b) and antibonding (a) bond orbitals (BO) for the k-th geminal (which are also the Dirac orbitals for this geminal):

\{ \begin{matrix}b_{kσ} = x_{k}l_{kσ} + y_{k}r_{kσ} \\ a_{kσ} = −y_{k}l_{kσ} + x_{k}r_{kσ} \end{matrix} ; (x_{k}^{2} + y_{k}^{2} = 1).  (24)

for constructing geminals:

g_{k}^{+} = U_{k}b_{kα}^{+}b_{kβ}^{+} + V_{k}a_{kα}^{+}a_{kβ}^{+}; (U_{k}^{2} + V_{k}^{2} = 1).  (25)

The APSLG wave function remains unchanged since the coefficient sets (U_{k}, V_{k}; x_{k}, y_{k}) and (u_{k}, w_{k}, v_{k}) are uniquely related:

\{ \begin{matrix}u_{m} = U_{m}y_{m}^{2} + V_{m}x_{m}^{2} \\ v_{m} = U_{m}x_{m}^{2} + V_{m}y_{m}^{2} \\ w_{m} = (U_{m} − V_{m}) x_{m}y_{m} \end{matrix} .  (26)

Using these bond orbitals we construct the M-multiplier of the charge transfer states |i⟩ in the form:

\begin{matrix}b_{kσ}^{+} ∏_{l≠k} g_{l}^{+} | 0⟩, a_{kσ}^{+} ∏_{l≠k} g_{l}^{+} | 0⟩, \\ b_{kσ}^{+}r_{k−σ}^{+}l_{k−σ}^{+} ∏_{l≠k} g_{l}^{+} | 0⟩, a_{kσ}^{+}r_{k−σ}^{+}l_{k−σ}^{+} ∏_{l≠k} g_{l}^{+} | 0⟩. \end{matrix}  (27)

The IP’s and the EA’s for the bond (lone pair) states within the APSLG-MINDO/3 approximation have the form:

\begin{matrix}I_{kb} = W_{1k}^{r} (y_{k}^{2} − 2P_{k}^{rr}) + W_{1k}^{l} (x_{k}^{2} − P_{k}^{ll}) + 2W_{1k}^{rl} (x_{k}y_{k} − 2P_{k}^{rl}) − \\ −W_{2k}^{r}u_{k}^{2} − W_{2k}^{l}v_{k}^{2} − 2W_{2k}^{rl}w_{k}^{2}, \\ I_{ka} = W_{1k}^{r} (x_{k}^{2} − 2P_{k}^{rr}) + W_{1k}^{l} (y_{k}^{2} − P_{k}^{ll}) − 2W_{1k}^{rl} (x_{k}y_{k} + 2P_{k}^{rl}) − \\ −W_{2k}^{r}u_{k}^{2} − W_{2k}^{l}v_{k}^{2} − 2W_{2k}^{rl}w_{k}^{2}, \\ −A_{kb} = W_{1k}^{r} (1 + y_{k}^{2} − 2P_{k}^{rr}) + W_{1k}^{l} (1 + x_{k}^{2} − 2P_{k}^{ll}) + \\ +2W_{1k}^{rl} (x_{k}y_{k} − 2P_{k}^{rl}) + W_{2k}^{r} (y_{k}^{2} − u_{k}^{2}) + \\ +W_{2k}^{l} (x_{k}^{2} − v_{k}^{2}) + 2W_{2k}^{rl} (1 − w_{k}^{2}) , \\ −A_{ka} = W_{1k}^{r} (1 + x_{k}^{2} − 2P_{k}^{rr}) + W_{1k}^{l} (1 + y_{k}^{2} − 2P_{k}^{ll}) − \\ −2W_{1k}^{rl} (x_{k}y_{k} + 2P_{k}^{rl}) + W_{2k}^{r} (x_{k}^{2} − u_{k}^{2}) + \\ +W_{2k}^{l} (y_{k}^{2} − v_{k}^{2}) + 2W_{2k}^{rl} (1 − w_{k}^{2}) . \end{matrix}  (28)

where these quantities are expressed through the parameters of the effective Hamiltonian for the k-th bond in the APSLG-MINDO/3 approximation, which can be written as

\begin{matrix}H_{k}^{eff} = W_{1k}^{r} ∑_{σ} r_{kσ}^{+}r_{kσ} + W_{1k}^{l} ∑_{σ} l_{kσ}^{+}l_{kσ} + W_{1k}^{rl} ∑_{σ} (r_{kσ}^{+}l_{kσ} + l_{kσ}^{+}r_{kσ}) + \\ +W_{2k}^{r}r_{kα}^{+}r_{kβ}^{+}r_{kβ}r_{kα} + W_{2k}^{l}l_{kα}^{+}l_{kβ}^{+}l_{kβ}l_{kα} + W_{2k}^{rl} ∑_{σ} r_{kσ}^{+}l_{k−σ}^{+}l_{k−σ}r_{kσ}, \end{matrix}  (29)

with

\begin{matrix}W_{1k}^{r} = (U_{k}^{A_{k}} − ∑_{B′≠A_{k}} γ_{A_{k}B′}Z_{B′}) + \\ + ∑_{t_{n}∈A_{k},n≠k} Y_{r_{k}r_{k}t_{n}t_{n}}^{A_{k}}P_{n}^{tt} + 2 ∑_{B′≠A_{k}} γ_{A_{k}B′} ∑_{t_{n}∈B′,n≠k} P_{n}^{tt}, \\ W_{1k}^{l} = (U_{k}^{B_{k}} − ∑_{B′≠B_{k}} γ_{B_{k}B′}Z_{B′}) + \\ + ∑_{t_{n}∈A_{k},n≠k} Y_{l_{k}l_{k}t_{n}t_{n}}^{B_{k}}P_{n}^{tt} + 2 ∑_{B′≠B_{k}} γ_{B_{k}B′} ∑_{t_{n}∈B′,n≠k} P_{n}^{tt}, \\ W_{1k}^{rl} = −β_{k}, \\ W_{2k}^{r} = (r_{k}r_{k} | r_{k}r_{k})^{A_{k}} , \\ W_{2k}^{l} = (l_{k}l_{k} | l_{k}l_{k})^{B_{k}} , \\ W_{2k}^{rl} = γ_{A_{k}B_{k}}. \end{matrix}  (30)

The resolvent contribution from the first projection to the effective Hamiltonian of the R-subsystem can be written in the form including the ionized states of the R-subsystem:

\begin{matrix}⟨⟨P V^{rr}(q, E)P ⟩⟩_{M} = − ∑_{pp′∈R} ∑_{σ} ∑_{k∈M} ∑_{ij∈\{r,l\}} ∑_{f∈\{a,b\}} v_{pi_{k}}v_{p′j_{k}}× \\ (  ∑_{ρ∈ImO_{R}(N_{R}−1)} p_{σ}^{+} |ρ⟩ ⟨ρ| p′_{σ} \frac{θ_{if}^{k} θ_{jf}^{k}}{I_{ρ}−A_{kf} −g_{ρf_{k}}} + \\ + ∑_{ρ∈ImO_{R}(N_{R}+1)} p_{σ} |ρ⟩ ⟨ρ| p_{σ}^{′+} \frac{η_{if}^{k} η_{jf}^{k}}{I_{kf} −A_{ρ}−g_{f_{k}ρ}}), \end{matrix}  (31)

where we use the vacuum averages

\begin{matrix}θ_{if}^{k} = ⟨0| g_{k}i_{kα}f_{kβ}^{+}r_{kα}^{+}l_{kα}^{+} |0⟩ , \\ η_{if}^{k} = ⟨0| g_{k}i_{kα}^{+}f_{kβ}^{+} |0⟩ . \end{matrix}  (32)

which can be easily expressed in terms of the parameters of the APSLG wave function in the representations Eq. (10) and Eqs. (25), (26):

\begin{matrix}θ_{rb}^{k} = y_{k}w_{k} + x_{k}v_{k}, θ_{ra}^{k} = x_{k}w_{k} − y_{k}v_{k}, \\ θ_{lb}^{k} = −x_{k}w_{k} − y_{k}u_{k}, θ_{la}^{k} = y_{k}w_{k} − x_{k}u_{k}; \\ η_{rb}^{k} = y_{k}u_{k} + x_{k}w_{k}, η_{ra}^{k} = x_{k}u_{k} − y_{k}w_{k}, \\ η_{lb}^{k} = y_{k}w_{k} + x_{k}v_{k}, η_{la}^{k} = x_{k}w_{k} − y_{k}v_{k}. \end{matrix}  (33)

Now we consider rather cumbersome contribution which arise during the second projection:

⟨Φ_{00}^{M} | P W(q, E)P QR(ω)QP W(q, E)P | Φ_{00}^{M}⟩.  (34)

In order to do that we reconsider the substance of the notion of the quantum character of the R-subsystem and of the classical one of the M-subsystem. As it is mentioned in Ref. [13] the quantum character of a part of a molecular system manifests itself in its spectrum which possesses excited states in a narrow energy range close to its ground state. It makes it possible to observe several quantum states in experiment at least in principle. By contrast the characteristic of a classical part of molecular system is that its properties are determined by its ground state only so that the energies of its excited states are very high as compared to the energy range probed experimentally. In the present derivation we are interested in obtaining the effective Hamiltonian for the states of the R- (quantum) subsystem close to its ground state. We assume that the dependence of the resolvent on ω is weak and the ω values of interest are much smaller than the resolvent poles which are all higher that the first excitation energy in the M- (classical) subsystem which in its turn is higher than the excitation energies of interest in the R-subsystem. Therefore we can turn to the limit ω → 0 in the expression (34). The resolvent can be represented as

R(0) = lim_{ω→0} R(ω) = lim_{ω→0} ∑_{\begin{matrix}ρ \\ μ≠0 \end{matrix}} \frac{|ρ, μ⟩ ⟨ρ, μ|}{ω − ω_{ρμ}} = − ∑_{\begin{matrix}ρ \\ μ≠0 \end{matrix}} \frac{|ρ, μ⟩ ⟨ρ, μ|}{ϵ_{ρ} + ϵ_{μ}},  (35)

where the sum excludes the states having the ground state of the M-subsystem as the multiplier and ϵ_{ρ} and ϵ_{μ} are the energies of the excitations in the R- and M-subsystems, respectively. The contribution Eq.(34) can be expressed as:

⟨Φ_{00}^{M} | P W(q, E)P QR(0)QP W(q, E)P | Φ_{00}^{M}⟩ =
= ⟨Φ_{00}^{M} | P V^{c}P QR(0)QP V^{c}P | Φ_{00}^{M}⟩+
+⟨Φ_{00}^{M} | P V^{rr}P QR(0)QP V^{rr}P | Φ_{00}^{M}⟩+
+⟨Φ_{00}^{M} | P V^{R}P QR(0)QP V^{R}P | Φ_{00}^{M}⟩+
+2⟨Φ_{00}^{M} | P V^{c}P QR(0)QP V^{rr}P | Φ_{00}^{M}⟩+
+2⟨Φ_{00}^{M} | P V^{c}P QR(0)QP V^{R}P | Φ_{00}^{M}⟩+
+2⟨Φ_{00}^{M} | P V^{rr}P QR(0)QP V^{R}P | Φ_{00}^{M}⟩.

(36)

We will write these averages explicitly. Since the R-subsystem is quantum and thus has low-lying excited states but the M-subsystem is classical and its excitation energies are high and thus ϵ_{ρ} << ϵ_{μ} and ϵ_{ρ} can be neglected in Eq. (35) as compared to ϵ_{μ}. Using the symmetry with respect to the spin indices we obtain:

\begin{matrix}⟨Φ_{00}^{M} | P V^{c}P QR(0)QP V^{c}P | Φ_{00}^{M}⟩ = \\ = ∑_{pp′qq′∈R} ∑_{kn∈M} ∑_{ii′jj′∈\{r,l\}}[ ∑_{\begin{matrix}στ \\ σ′τ′ \end{matrix}} p_{σ}^{+}p′_{σ} (1 − |0_{R}⟩ ⟨0_{R}|) q_{τ}^{+}q′_{τ}Π_{i_{kσ′}i′_{nσ′}j_{nτ′}j′_{kτ′}}^{M}× \\ ×((pp′ | i_{k}i′_{n}) − δ_{σσ′}(pi′_{n} | i_{k}p′))((qq′ | j_{n}j′_{k}) − δ_{ττ′}(qj′_{k} | j_{n}q′))+ \\ + ∑_{σ} p_{σ}^{+}p′_{−σ} (1 − |0_{R}⟩ ⟨0_{R}|) q_{−σ}^{+}q′_{σ}(pi′_{n} | i_{k}p′)(qj′_{k} | j_{n}q′)Π_{i_{k−σ}i′_{nσ}j_{nσ}j′_{k−σ}}^{M}], \end{matrix}  (37)

where the state |0_{R}⟩ is the ground state of the free R-subsystem. In this expression we use the zero frequency polarization propagators of the M-subsystem [12], which are defined by the expression

Π_{mm′nn′}^{M} = − ∑_{μ≠0} \frac{⟨Φ_{00}^{M} |m^{+}m′| μ⟩ ⟨μ |n^{+}n′| Φ_{00}^{M}⟩}{ϵ_{μ}}.  (38)

to which the excitations in the M-subsystem contribute. The APSLG form of the ground state wave function implies specific classification for these excitations: they are either intrabond singlet - singlet or singlet - triplet excitations or the interbond one- and two-electron transfers. We present the explicit expressions for the contributions to the polarization propagator from the excitations of different types.

1. Π_{i_{kσ}i′_{kσ}j_{kτ} j′_{kτ}}^{M}. The excited state |μ⟩ is one of the states |μ_{1,2}⟩ = g_{k}^{(1,2)+} ∏_{k′≠k} g_{k′}^{+} |0⟩

(these are excited singlet configurations of the k-th geminal with others unchanged) and |μ_{3}⟩ = \frac{1}{\sqrt{2}}(r_{kα}^{+}l_{kβ}^{+} − l_{kα}^{+}r_{kβ}^{+}) ∏_{k′≠k} g_{k′}^{+} |0⟩ (the triplet configuration with

s_{z} = 0 of the k-th geminal with others unchanged). Therefore,

Π_{i_{kσ}i′_{kσ}j_{kτ} j′_{kτ}}^{M} = − ∑_{s∈\{1÷3\}} \frac{⟨Φ_{00}^{M} |i_{kσ}^{+}i′_{kσ}| μ_{s}⟩ ⟨μ_{s} |j_{kτ}^{+}j′_{kτ} | Φ_{00}^{M}⟩}{ϵ_{μ_{s}}}.  (39)

The energies of excitations ϵ_{μ_{s}} in the APSLG-MINDO/3 approximation Ref. [14] are:

\begin{matrix}ϵ_{μ_{s=1,2}} = 2W_{1k}^{r} (u_{k}^{(s)2} + w_{k}^{(s)2} − P_{k}^{rr}) + \\ +2W_{1k}^{l} (v_{k}^{(s)2} + w_{k}^{(s)2} − P_{k}^{ll}) + \\ +W_{2k}^{r} (u_{k}^{(s)2} − u_{k}^{2}) + W_{2k}^{l} (v_{k}^{(s)2} − v_{k}^{2}) + \\ +4W_{1k}^{rl} [(u_{k}^{(s)} + v_{k}^{(s)}) w_{k}^{(s)} − P_{k}^{rl}] + 2W_{2k}^{rl} (w_{k}^{(s)2} − w_{k}^{2}) , \end{matrix}  (40)

\begin{matrix}ϵ_{μ_{3}} = (W_{1k}^{r} − W_{1k}^{l}) (v_{k}^{2} − u_{k}^{2}) − \\ −W_{2k}^{r}u_{k}^{2} − W_{2k}^{l}v_{k}^{2} − 4W_{1k}^{rl}P_{k}^{rl} + W_{2k}^{rl} (u_{k}^{2} + v_{k}^{2}) . \end{matrix}  (41)

The matrix elements entering the expression for the polarization propagator are

\begin{matrix}⟨Φ_{00}^{M} |r_{kσ}^{+}r_{kσ}| μ_{s=1,2}⟩ = u_{k}u_{k}^{(s)} + w_{k}w_{k}^{(s)}, \\ ⟨Φ_{00}^{M} |l_{kσ}^{+}l_{kσ}| μ_{s=1,2}⟩ = v_{k}v_{k}^{(s)} + w_{k}w_{k}^{(s)}, \\ ⟨Φ_{00}^{M} |r_{kσ}^{+}l_{kσ}| μ_{s=1,2}⟩ = u_{k}w_{k}^{(s)} + w_{k}v_{k}^{(s)}, \\ ⟨Φ_{00}^{M} |l_{kσ}^{+}r_{kσ}| μ_{s=1,2}⟩ = w_{k}u_{k}^{(s)} + v_{k}w_{k}^{(s)}; \end{matrix}  (42)

\begin{matrix}⟨Φ_{00}^{M} |r_{kα}^{+}r_{kα}| μ_{3}⟩ = \frac{w_{k}}{\sqrt{2}}, ⟨Φ_{00}^{M} |l_{kα}^{+}l_{kα}| μ_{3}⟩ = −\frac{w_{k}}{\sqrt{2}}, \\ ⟨Φ_{00}^{M} |r_{kα}^{+}l_{kα}| μ_{3}⟩ = \frac{−u_{k}}{\sqrt{2}} , ⟨Φ_{00}^{M} |l_{kα}^{+}r_{kα}| μ_{3}⟩ = \frac{v_{k}}{\sqrt{2}} \\ ⟨Φ_{00}^{M} |i_{kβ}^{+}j_{kβ}| μ_{3}⟩ = − ⟨Φ_{00}^{M} |i_{kα}^{+}j_{kα}| μ_{3}⟩ . \end{matrix}  (43)

For real orbitals the following ⟨μ_{s} |i^{+}j| Φ_{00}^{M}⟩ = ⟨Φ_{00}^{M} |j^{+}i| μ_{s}⟩ holds.

2. Π_{i_{k−σ}i′_{kσ}j_{kσ}j′_{k−σ}}^{M}. The excited state |μ_{4}⟩ = r_{kσ}^{+}l_{kσ}^{+} ∏_{k′≠k} g_{k′}^{+} |0⟩ (the triplet con-

figuration with s_{z} = 1 for σ = α and s_{z} = −1 for σ = β in the k-th geminal with others unchanged).

⟨Φ_{00}^{M} |i_{k−σ}^{+}i′_{kσ}| μ_{4}⟩ = ⟨0 |g_{k}i_{k−σ}^{+}i′_{kσ}r_{kσ}^{+}l_{kσ}^{+}| 0⟩ .  (44)

The energies of such excitations obviously equal to respective energies of excitations to the triplet states with s_{z} = 0, i.e. ϵ_{μ_{4}} = ϵ_{μ_{3}}. The matrix elements equal to

\begin{matrix}⟨Φ_{00}^{M} |r_{kα}^{+}r_{kβ}| μ_{4}⟩ = w_{k}, ⟨Φ_{00}^{M} |l_{kα}^{+}l_{kβ}| μ_{4}⟩ = −w_{k}, \\ ⟨Φ_{00}^{M} |l_{kα}^{+}r_{kβ}| μ_{4}⟩ = v_{k}, ⟨Φ_{00}^{M} |r_{kα}^{+}l_{kβ}| μ_{4}⟩ = −u_{k} \\ ⟨Φ_{00}^{M} |i_{kβ}^{+}j_{kα}| μ_{4}⟩ = − ⟨Φ_{00}^{M} |i_{kα}^{+}j_{kβ}| μ_{4}⟩ . \end{matrix}  (45)

3. Π_{i_{kσ}i′_{nσ}j_{nτ} j′_{kτ}}^{M} (n ≠ k). This polarization propagator differs from zero only for

σ = τ. The excited state involved is |μ_{5}⟩ = h_{n−σ}^{+}r_{nσ}^{+}l_{nσ}^{+}f_{k−σ}^{+} ∏_{k′≠k,n} g_{k′}^{+} |0⟩, where f

and h are either b or a BO’s. The energy of the excitation is estimated as

ϵ_{5} = I_{kf} − A_{nh} − g_{f_{k}h_{n}}.  (46)

The IP’s and the EA’s (I_{kf}, A_{nh}) are given above. Therefore, we specify here only g_{f_{k}h_{l}}. It is convenient to introduce new quantities c_{fi}, which are amplitudes of the i-th HO (r or l) in the BO f (a or b) given by Eq. (24). Then

\begin{matrix}g_{f_{k}h_{n}} = ∑_{ij∈\{r,l\}}[Y_{i_{k}i_{k}j_{n}j_{n}} (2P_{k}^{ii} (1 + c_{f_{k}i_{k}}^{2}) + 2P_{n}^{jj}c_{h_{n}j_{n}}^{2}) − \\ − (i_{k}i_{k} | j_{n}j_{n}) (1 + c_{f_{k}i_{k}}^{2}) c_{h_{n}j_{n}}^{2} + (i_{k}j_{n} | i_{k}j_{n}) c_{f_{k}i_{k}}^{2}c_{h_{n}j_{n}}^{2}]. \end{matrix}  (47)

The required matrix elements are

⟨Φ_{00}^{M} |i_{kσ}^{+}i′_{nσ}| μ_{5}⟩ = η_{if}^{k}θ_{i′h}^{n}.  (48)

4. Π_{i_{k−σ}i′_{nσ}j_{nσ}j′_{k−σ}}^{M}(n ≠ k). The excited state is |μ_{6}⟩ = h_{n−σ}^{+}r_{nσ}^{+}l_{nσ}^{+}f_{kσ}^{+} ∏_{k′≠k,n} g_{k′}^{+} |0⟩,

where f and h are either b or a. The energy of the excitation is estimated as

ϵ_{6} = I_{kf} − A_{nh} − g′_{f_{k}h_{n}},  (49)

where

\begin{matrix}g′_{f_{k}h_{n}} = ∑_{ij∈\{r,l\}}[Y_{i_{k}i_{k}j_{n}j_{n}} (2P_{k}^{ii} (1 + c_{f_{k}i_{k}}^{2}) + 2P_{n}^{jj}c_{h_{n}j_{n}}^{2}) − \\ − (i_{k}i_{k} | j_{n}j_{n}) (1 + c_{f_{k}i_{k}}^{2}) c_{h_{n}j_{n}}^{2} + (i_{k}j_{n} | i_{k}j_{n}) c_{h_{n}j_{n}}^{2}]. \end{matrix}  (50)

The required matrix elements are

⟨Φ_{00}^{M} |i_{k−σ}^{+}i′_{nσ}| μ_{6}⟩ = −η_{if}^{k}θ_{i′h}^{n}.  (51)

If the ZDO approximation is used for the two-center Coulomb interaction parameters the contribution Eq. (34) into the effective Hamiltonian for the R-subsystem simplifies significantly:

\begin{matrix}⟨Φ_{00}^{M} | P V^{c}P QR(0)QP V^{c}P | Φ_{00}^{M}⟩ = \\ = ∑_{pq∈R} ∑_{k∈M} ∑_{ij∈\{r,l\}}[ ∑_{\begin{matrix}στ \\ σ′τ′ \end{matrix}} p_{σ}^{+}p_{σ} (1 − |0_{R}⟩ ⟨0_{R}|) q_{τ}^{+}q_{τ}Π_{i_{kσ′}i_{kσ′}j_{kτ′}j_{kτ′}}^{M}× \\ ×((pp | i_{k}i_{k}) − δ_{σσ′}(pi_{k} | i_{k}p))((qq | j_{k}j_{k}) − δ_{ττ′}(qj_{k} | j_{k}q))+ \\ + ∑_{σ} p_{σ}^{+}p_{−σ} (1 − |0_{R}⟩ ⟨0_{R}|) q_{−σ}^{+}q_{σ}(pi_{k} | i_{k}p)(qj_{k} | j_{k}q)Π_{i_{k−σ}i_{kσ}j_{kσ}j_{k−σ}}^{M}]. \end{matrix}  (52)

The next contribution to the effective Hamiltonian can be presented as a sum of one-geminal and two-geminal contributions:

\begin{matrix}⟨Φ_{00}^{M} | P V^{rr}P QR(0)QP V^{rr}P | Φ_{00}^{M}⟩ = \\ ∑_{k}⟨Φ_{00}^{M} | P V^{rr}P QR(0)QP V^{rr}P | Φ_{00}^{M}⟩_{k}+ \\ + ∑_{k≠n}⟨Φ_{00}^{M} | P V^{rr}P QR(0)QP V^{rr}P | Φ_{00}^{M}⟩_{kn}. \end{matrix}  (53)

The one-geminal contribution has the form

\begin{matrix}⟨Φ_{00}^{M} | P V^{rr}P QR(0)QP V^{rr}P | Φ_{00}^{M}⟩_{k} = \\ = ∑_{\begin{matrix}pp′∈R \\ qq′∈R \end{matrix}} ∑_{\begin{matrix}ii′jj′∈\{r,l\} \\ fh∈\{a,b\} \end{matrix}} ∑_{\begin{matrix}στ \\ ςς′ \end{matrix}} v_{pi_{k}}v_{p′i′_{k}}v_{qj_{k}}v_{q′j′_{k}}× \\ ×[  ∑_{ρ_{1},ρ_{2}∈ImO_{R}(N_{R}+1)} ^{\frac{A_{pp′qq′}^{στςς′}(ρ_{1}ρ_{2})Ξ_{i_{kσ}i{′ \atop kς} j_{kς′} j{′ \atop kτ}}^{f_{k−σ}h_{k−τ}}}{(I_{kf} −A_{ρ_{1}}−g_{f_{k}ρ_{1}})(I_{kh}−A_{ρ_{2}}−g_{h_{k}ρ_{2}})}}− \\ − ∑_{\begin{matrix}ρ_{1}∈ImO_{R}(N_{R}+1) \\ ρ_{2}∈ImO_{R}(N_{R}−1) \end{matrix}} (−1)^{δ_{jj′}} ^{\frac{B_{pp′qq′}^{στςς′}(ρ_{1}ρ_{2})Ξ_{i_{kσ}i{′ \atop kς}_{j_{kς′}} j{′ \atop kτ}}^{f_{k−σ}h_{k−τ}}}{(I_{kf} −A_{ρ_{1}}−g_{f_{k}ρ_{1}})(I_{ρ_{2}}−A_{kh}−g_{ρ_{2}h_{k}})}}− \\ − ∑_{\begin{matrix}ρ_{1}∈ImO_{R}(N_{R}−1) \\ ρ_{2}∈ImO_{R}(N_{R}+1) \end{matrix}} (−1)^{δ_{ii′}} ^{\frac{C_{pp′qq′}^{στςς′}(ρ_{1}ρ_{2})Ξ_{_{i_{kσ}}i{′ \atop kς} j_{kς′} j{′ \atop kτ}}^{f_{k−σ}h_{k−τ}}}{(I_{ρ_{1}}−A_{kf} −g_{ρ_{1}f_{k}})(I_{kh}−A_{ρ_{2}}−g_{h_{k}ρ_{2}})}}+ \\ + ∑_{ρ_{1},ρ_{2}∈ImO_{R}(N_{R}−1)}(−1)^{δ_{ii′}+δ_{jj′}} ^{\frac{D_{pp′qq′}^{στςς′}(ρ_{1}ρ_{2})Ξ_{_{i_{kσ}}i{′ \atop kς}_{j_{kς′}} j{′ \atop kτ}}^{f_{k−σ}h_{k−τ}}}{(I_{ρ_{1}}−A_{kf} −g_{ρ_{1}f_{k}})(I_{ρ_{2}}−A_{kh}−g_{ρ_{2}h_{k}})}}], \end{matrix}  (54)

while the two-geminal contributions have the form

\begin{matrix}⟨Φ_{00}^{M} | P V^{rr}P QR(0)QP V^{rr}P | Φ_{00}^{M}⟩_{kn} = \\ = ∑_{\begin{matrix}pp′∈R \\ qq′∈R \end{matrix}} ∑_{\begin{matrix}ii′jj′∈\{r,l\} \\ f∈\{a,b\} \end{matrix}} ∑_{σς} v_{pi_{k}}v_{p′i′_{n}}v_{qj_{n}}v_{q′j′_{k}}× \\ ×[  ∑_{ρ_{1},ρ_{2}∈ImO_{R}(N_{R}+1)} ^{\frac{A_{pp′qq′}^{σσςς}(ρ_{1}ρ_{2})Ξ_{i_{kσ}i{′ \atop nς} jnς j{′ \atop kσ}}^{f_{k−σ}f_{k−σ}}}{(I_{kf} −A_{ρ_{1}}−g_{f_{k}ρ_{1}})(I_{kf} −A_{ρ_{2}}−g_{f_{k}ρ_{2}})}}+ \\ + ∑_{ρ_{1},ρ_{2}∈ImO_{R}(N_{R}−1)} ^{\frac{D_{pp′qq′}^{σσςς}(ρ_{1}ρ_{2})Ξ_{i{′ \atop nς} i_{kσ}j{′ \atop kσ}jnς}^{f_{k−σ}f_{k−σ}}}{(I_{ρ_{1}}−A_{kf} −g_{ρ_{1}f_{k}})(I_{ρ_{2}}−A_{kf} −g_{ρ_{2}f_{k}})}}]− \\ + ∑_{\begin{matrix}pp′∈R \\ qq′∈R \end{matrix}} ∑_{\begin{matrix}ii′jj′∈\{r,l\} \\ fh∈\{a,b\} \end{matrix}} ∑_{σς} v_{pi_{k}}v_{p′i′_{n}}v_{qj_{k}}v_{q′j′_{n}}(−1)^{δ_{σς}}× \\ ×[  ∑_{\begin{matrix}ρ_{1}∈ImO_{R}(N_{R}+1) \\ ρ_{2}∈ImO_{R}(N_{R}−1) \end{matrix}} ^{\frac{B_{pp′qq′}^{σσςς}(ρ_{1}ρ_{2})Ξ_{i_{kσ}i{′ \atop nς} j{′ \atop nς} j_{kσ}}^{f_{k−σ}h_{n−ς}}}{(I_{kf} −A_{ρ_{1}}−g_{f_{k}ρ_{1}})(I_{ρ_{2}}−A_{nh}−g_{ρ_{2}hn})}}+ \\ + ∑_{\begin{matrix}ρ_{1}∈ImO_{R}(N_{R}−1) \\ ρ_{2}∈ImO_{R}(N_{R}+1) \end{matrix}} ^{\frac{C_{pp′qq′}^{σσςς}(ρ_{1}ρ_{2})Ξ_{i{′ \atop nς} i_{kσ}j_{kσ}j{′ \atop nς}}^{f_{k−σ}h_{n−ς}}}{(I_{ρ_{1}}−A_{kf} −g_{ρ_{1}f_{k}})(I_{nh}−A_{ρ_{2}}−g_{hnρ_{2}})}}], \end{matrix}  (55)

where ĩ = r for i = l and ĩ = l for i = r. Operators can be represented as

\begin{matrix}A_{pp′qq′}^{στςς′} (ρ_{1}ρ_{2}) = δ_{ςσ}δ_{ς′τ}p_{σ} |ρ_{1}⟩ ⟨ρ_{1}| p_{σ}^{′+} (1 − |0_{R}⟩ ⟨0_{R}|) q_{τ} |ρ_{2}⟩ ⟨ρ_{2}| q_{τ}^{′+}+ \\ +δ_{ςς′}δ_{στ}δ_{ς−σ}p_{σ} |ρ_{1}⟩ ⟨ρ_{1}| p_{−σ}^{′+} (1 − |0_{R}⟩ ⟨0_{R}|) q_{−σ} |ρ_{2}⟩ ⟨ρ_{2}| q_{σ}^{′+}, \\ B_{pp′qq′}^{στςς′} (ρ_{1}ρ_{2}) = δ_{ςσ}δ_{ς′τ}p_{σ} |ρ_{1}⟩ ⟨ρ_{1}| p_{σ}^{′+} (1 − |0_{R}⟩ ⟨0_{R}|) q_{τ}^{+} |ρ_{2}⟩ ⟨ρ_{2}| q′_{τ}− \\ −δ_{ςς′}δ_{στ}δ_{ς−σ}p_{σ} |ρ_{1}⟩ ⟨ρ_{1}| p_{−σ}^{′+} (1 − |0_{R}⟩ ⟨0_{R}|) q_{−σ}^{+} |ρ_{2}⟩ ⟨ρ_{2}| q′_{σ}, \\ C_{pp′qq′}^{στςς′} (ρ_{1}ρ_{2}) = δ_{ςσ}δ_{ς′τ}p_{σ}^{+} |ρ_{1}⟩ ⟨ρ_{1}| p′_{σ} (1 − |0_{R}⟩ ⟨0_{R}|) q_{τ} |ρ_{2}⟩ ⟨ρ_{2}| q_{τ}^{′+}− \\ −δ_{ςς′}δ_{στ}δ_{ς−σ}p_{σ}^{+} |ρ_{1}⟩ ⟨ρ_{1}| p′_{−σ} (1 − |0_{R}⟩ ⟨0_{R}|) q_{−σ} |ρ_{2}⟩ ⟨ρ_{2}| q_{σ}^{′+}, \\ D_{pp′qq′}^{στςς′} (ρ_{1}ρ_{2}) = δ_{ςσ}δ_{ς′τ}p_{σ}^{+} |ρ_{1}⟩ ⟨ρ_{1}| p′_{σ} (1 − |0_{R}⟩ ⟨0_{R}|) q_{τ}^{+} |ρ_{2}⟩ ⟨ρ_{2}| q′_{τ}+ \\ +δ_{ςς′}δ_{στ}δ_{ς−σ}p_{σ}^{+} |ρ_{1}⟩ ⟨ρ_{1}| p′_{−σ} (1 − |0_{R}⟩ ⟨0_{R}|) q_{−σ}^{+} |ρ_{2}⟩ ⟨ρ_{2}| q′_{σ}. \end{matrix}  (56)

Also we introduce the quantities

Ξ_{ii′jj′}^{fh} = − ∑_{μ≠0} \frac{⟨Φ_{00}^{M} |i^{+}f^{+}fi′| μ⟩ ⟨μ |j^{+}h^{+}hj′| Φ_{00}^{M}⟩}{ϵ_{μ}},  (57)

which in general case can not be reduced to the polarization propagators. Now we consider different cases (The energies of excitations were determined above).

1. Ξ_{i_{kσ}i′_{kσ}j_{kσ}j′_{kσ}}^{f_{k−σ}h_{k−σ}}. The excited state |μ⟩ is one of |μ_{1,2,3}⟩, defined above. We have to determine only matrix elements.

⟨Φ_{00}^{M} |i_{kσ}^{+}f_{k−σ}^{+}f_{k−σ}i′_{kσ}| μ_{s=1,2}⟩ = θ_{if}^{k}θ_{i′f}^{k}(s),  (58)

where θ_{if}^{k}(s) coincides with θ_{if}^{k} defined above (see Eq. (32)), but with u_{k}, v_{k}, and w_{k} changed to u_{k}^{(s)}, v_{k}^{(s)}, and w_{k}^{(s)}.

⟨Φ_{00}^{M} |i_{kσ}^{+}f_{k−σ}^{+}f_{k−σ}i′_{kσ}| μ_{3}⟩ = \frac{1}{\sqrt{2}}θ_{if}^{k} (δ_{i′r}c_{f_{k}l_{k}} − δ_{i′l}c_{f_{k}r_{k}}) .  (59)

2. Ξ_{i_{k−σ}i′_{kσ}j_{kσ}j′_{k−σ}}^{f_{kσ}h_{kσ}}. The excited state is |μ_{4}⟩. The matrix elements are

⟨Φ_{00}^{M} |i_{k−σ}^{+}f_{kσ}^{+}f_{kσ}i′_{kσ}| μ_{4}⟩ = (−1)^{δ_{σα}}θ_{if}^{k} (δ_{i′r}c_{f_{k}l_{k}} − δ_{i′l}c_{f_{k}r_{k}}) .  (60)

3. Ξ_{i_{kσ}i′_{nσ}j_{nσ}j′_{kσ}}^{f_{k−σ}f_{k−σ}}. The excited state is of |μ_{5}⟩ type. The required matrix elements are

⟨Φ_{00}^{M} |i_{kσ}^{+}f_{k−σ}^{+}f_{k−σ}i′_{nσ}| μ_{5}⟩ = η_{if}^{k}θ_{i′h}^{n}.  (61)

4. Ξ_{i′_{nσ}i_{kσ}j′_{kσ}j_{nσ}}^{f_{k−σ}f_{k−σ}}. The excited state is of |μ_{5}⟩ type: f_{k−σ}^{+}r_{kσ}^{+}l_{kσ}^{+}h_{n−σ}^{+} ∏_{k′≠k,n} g_{k′}^{+} |0⟩.

The required matrix elements are

⟨Φ_{00}^{M} |i_{nσ}^{′+}f_{k−σ}^{+}f_{k−σ}i_{kσ}| μ_{5}⟩ = θ_{if}^{k}η_{i′f}^{n}.  (62)

5. Ξ_{i_{k−σ}i′_{nσ}j_{nσ}j′_{k−σ}}^{f_{kσ}f_{kσ}}. The excited state is of |μ_{6}⟩ type. The required matrix elements are

⟨Φ_{00}^{M} |i_{k−σ}^{+}f_{kσ}^{+}f_{kσ}i′_{nσ}| μ_{6}⟩ = −η_{if}^{k}θ_{i′h}^{n}.  (63)

6. Ξ_{i′_{nσ}i_{k−σ}j′_{k−σ}j_{nσ}}^{f_{kσ}f_{kσ}}. The excited state is of |μ_{6}⟩ type: f_{kσ}^{+}r_{k−σ}^{+}l_{k−σ}^{+}h_{n−σ}^{+} ∏_{k′≠k,n} g_{k′}^{+} |0⟩.

The required matrix elements are

⟨Φ_{00}^{M} |i_{nσ}^{′+}f_{kσ}^{+}f_{kσ}i_{k−σ}| μ_{6}⟩ = −θ_{if}^{k}η_{i′h}^{n}.  (64)

Other quantities Ξ do not contain new averages and can be obtained using expressions written above.

Also we can note that

Π_{i_{kσ}i′_{kσ}j_{kσ}j′_{kσ}}^{M} = ∑_{fh∈{r,l}} Ξ_{i_{kσ}i′_{kσ}j_{kσ}j′_{kσ}}^{f_{k−σ}h_{k−σ}}  (65)

and analogous expressions can be obtained for other types of Ξ, i.e., the summation over upper indices of Ξ’s gives the polarization propagator with respective lower indices. Therefore, if we use the approximation that the energies of electron transfer from (or to) the k-th bond of the M-subsystem to (or from) R-subsystem do not depend on the type of the resultant state of the bond (or equivalently consider only one excited state of each type for ) we obtain this contribution to the effective Hamiltonian for the R-subsystem to be expressed in terms of the polarization propagators:

\begin{matrix}⟨Φ_{00}^{M} | P V^{rr}P QR(0)QP V^{rr}P | Φ_{00}^{M}⟩_{k} = \\ = ∑_{\begin{matrix}pp′∈R \\ qq′∈R \end{matrix}} ∑_{ii′jj′∈\{r,l\}} ∑_{\begin{matrix}στ \\ ςς′ \end{matrix}} v_{pi_{k}}v_{p′i′_{k}}v_{qj_{k}}v_{q′j′_{k}}× \\ ×[  ∑_{ρ_{1},ρ_{2}∈ImO_{R}(N_{R}+1)} ^{\frac{A_{pp′qq′}^{στςς′}(ρ_{1}ρ_{2})Π_{i_{kσ}i{′ \atop kς} j_{kς′} j{′ \atop kτ}}^{M}}{(I_{k}−A_{ρ_{1}}−g_{kρ_{1}})(I_{k}−A_{ρ_{2}}−g_{kρ_{2}})}}− \\ − ∑_{\begin{matrix}ρ_{1}∈ImO_{R}(N_{R}+1) \\ ρ_{2}∈ImO_{R}(N_{R}−1) \end{matrix}} (−1)^{δ_{jj′}} ^{\frac{B_{pp′qq′}^{στςς′}(ρ_{1}ρ_{2})Π_{i_{kσ}i{′ \atop kς}_{j_{kς′}} j{′ \atop kτ}}^{M}}{(I_{k}−A_{ρ_{1}}−g_{kρ_{1}})(I_{ρ_{2}}−A_{k}−g_{ρ_{2}k})}}− \\ − ∑_{\begin{matrix}ρ_{1}∈ImO_{R}(N_{R}−1) \\ ρ_{2}∈ImO_{R}(N_{R}+1) \end{matrix}} (−1)^{δ_{ii′}} ^{\frac{C_{pp′qq′}^{στςς′}(ρ_{1}ρ_{2})Π_{_{i_{kσ}}i{′ \atop kς} j_{kς′} j{′ \atop kτ}}^{M}}{(I_{ρ_{1}}−A_{k}−g_{ρ_{1}k})(I_{k}−A_{ρ_{2}}−g_{kρ_{2}})}}+ \\ + ∑_{ρ_{1},ρ_{2}∈ImO_{R}(N_{R}−1)}(−1)^{δ_{ii′}+δ_{jj′}} ^{\frac{D_{pp′qq′}^{στςς′}(ρ_{1}ρ_{2})Π_{_{i_{kσ}}i{′ \atop kς}_{j_{kς′}} j{′ \atop kτ}}^{M}}{(I_{ρ_{1}}−A_{k}−g_{ρ_{1}k})(I_{ρ_{2}}−A_{k}−g_{ρ_{2}k})}}], \end{matrix}  (66)

and

\begin{matrix}⟨Φ_{00}^{M} | P V^{rr}P QR(0)QP V^{rr}P | Φ_{00}^{M}⟩_{kn} = \\ = ∑_{\begin{matrix}pp′∈R \\ qq′∈R \end{matrix}} ∑_{ii′jj′∈\{r,l\}} ∑_{σς} v_{pi_{k}}v_{p′i′_{n}}v_{qj_{n}}v_{q′j′_{k}}× \\ ×[  ∑_{ρ_{1},ρ_{2}∈ImO_{R}(N_{R}+1)} \frac{A_{pp′qq′}^{σσςς}(ρ_{1}ρ_{2})Π_{i_{kσ}i{′ \atop nς} jnς j{′ \atop kσ}}^{M}}{(I_{k}−A_{ρ_{1}}−g_{kρ_{1}})(I_{k}−A_{ρ_{2}}−g_{kρ_{2}})}+ \\ + ∑_{ρ_{1},ρ_{2}∈ImO_{R}(N_{R}−1)} \frac{D_{pp′qq′}^{σσςς}(ρ_{1}ρ_{2})Π_{i{′ \atop nς} i_{kσ}j{′ \atop kσ}jnς}^{M}}{(I_{ρ_{1}}−A_{k}−g_{ρ_{1}k})(I_{ρ_{2}}−A_{k}−g_{ρ_{2}k})}]− \\ + ∑_{\begin{matrix}pp′∈R \\ qq′∈R \end{matrix}} ∑_{\begin{matrix}ii′jj′∈\{r,l\} \\ fh∈\{a,b\} \end{matrix}} ∑_{σς} v_{pi_{k}}v_{p′i′_{n}}v_{qj_{k}}v_{q′j′_{n}}(−1)^{δ_{σς}}× \\ ×[  ∑_{\begin{matrix}ρ_{1}∈ImO_{R}(N_{R}+1) \\ ρ_{2}∈ImO_{R}(N_{R}−1) \end{matrix}} \frac{B_{pp′qq′}^{σσςς}(ρ_{1}ρ_{2})Π_{i_{kσ}i{′ \atop nς} j{′ \atop nς} j_{kσ}}^{M}}{(I_{k}−A_{ρ_{1}}−g_{kρ_{1}})(I_{ρ_{2}}−A_{n}−g_{ρ_{2}n})}+ \\ + ∑_{\begin{matrix}ρ_{1}∈ImO_{R}(N_{R}−1) \\ ρ_{2}∈ImO_{R}(N_{R}+1) \end{matrix}} \frac{C_{pp′qq′}^{σσςς}(ρ_{1}ρ_{2})Π_{i{′ \atop nς} i_{kσ}j_{kσ}j{′ \atop nς}}^{M}}{(I_{ρ_{1}}−A_{k}−g_{ρ_{1}k})(I_{n}−A_{ρ_{2}}−g_{nρ_{2}})}], \end{matrix}  (67)

The next contribution to the effective Hamiltonian for the R-subsystem is a cnumber. The operator V^{R} can be written as

\begin{matrix}V^{R} = − ∑_{mm′∈M} ∑_{σ} m_{σ}^{+}m′_{σ} ∑_{B∈R} V_{B}^{mm′}Z_{B}^{R} \\ ≈ − ∑_{A∈M} ∑_{m∈A∩M} ∑_{σ} m_{σ}^{+}m_{σ} ∑_{\begin{matrix}B∈R \\ B≠A \end{matrix}} γ_{AB}Z_{B}^{R} \end{matrix}  (68)

Respective contribution is

\begin{matrix}⟨Φ_{00}^{M} | P V^{R}P QR(0)QP V^{R}P | Φ_{00}^{M}⟩ = \\ = ∑_{στ} ∑_{kn∈M} ∑_{ii′jj′∈\{r,l\}} ∑_{BB′∈R} V_{B}^{i_{k}i′_{n}}V_{B′}^{j_{n}j′_{k}}Z_{B}^{R}Z_{B′}^{R}Π_{i_{kσ}i′_{nσ}j_{nτ} j′_{kτ}}^{M} ≈ \\ ≈ ∑_{στ} ∑_{k∈M} Π_{r_{kσ}r_{kσ}r_{kτ} r_{kτ}}^{M} ∑_{BB′∈R} Z_{B}^{R}Z_{B′}^{R}× \\ × ((1 − δ_{BA_{k}})γ_{BA_{k}} − (1 − δ_{BB_{k}})γ_{BB_{k}}) × \\ × ((1 − δ_{B′A_{k}})γ_{B′A_{k}} − (1 − δ_{B′B_{k}})γ_{B′B_{k}}) . \end{matrix}  (69)

From this approximate expression we can portion out the contribution which includes BB′ ∉ {A_{k}, B_{k}}:

∑_{στ} ∑_{k∈M} Π_{r_{kσ}r_{kσ}r_{kτ} r_{kτ}}^{M} ∑_{\begin{matrix}BB′∈R \\ BB′∉\{A_{k},B_{k}\} \end{matrix}} Z_{B}^{R}Z_{B′}^{R} (γ_{BA_{k}} − γ_{BB_{k}}) (γ_{B′A_{k}} − γ_{B′B_{k}}) . (70)

This contribution can be expressed through the bond polarizabilities. The two center integrals γ_{AB} can be approximated by their values calculated at the center of the k-th bond with the corrections linear in interatomic vectors:

γ_{BA_{k}} = e(V (R_{B}^{(k)}) − \frac{1}{2}(∇V (R_{B}^{(k)})R_{A_{k}B_{k}})),  (71)

where V (R_{B}^{(k)}) is the potential induced by a unit charge placed on the atom B at the center of the k-th bond; the correction contains the gradient of this potential. Substituting these expressions into Eq. (70) we get:

\begin{matrix}∑_{στ} ∑_{\begin{matrix}BB′∈R \\ BB′∉\{A_{k},B_{k}\} \end{matrix}} Z_{B}^{R}Z_{B′}^{R} ∑_{k∈M}(∇V (R_{B}^{(k)})μ_{A_{k}B_{k}}^{(k)})× \\ ×(∇V (R_{B′}^{(k)})μ_{A_{k}B_{k}}^{(k)})Π_{r_{kσ}r_{kσ}r_{kτ} r_{kτ}}^{M} = \\ = − ∑_{\begin{matrix}BB′∈R \\ BB′∉\{A_{k},B_{k}\} \end{matrix}} Z_{B}^{R}Z_{B′}^{R} ∑_{k∈M}(∇V (R_{B′}^{(k)}) |α̂^{(k)}(0)| ∇V (R_{B}^{(k)})), \end{matrix}  (72)

where the standard expressions reviewed in Ref. [17] for the bond polarizability tensors α̂^{(k)}(ω) for the k-th bond through the polarization propagator for the corresponding geminal and the bond dipole vector μ_{A_{k}B_{k}}^{(k)} of the k-th bond. The bond polarizabilities α̂^{(k)}(0) are tabulated, for example, in Ref. [18].

The cross term between the P V^{c}P and P V^{rr}P is a sum of one-geminal and two-geminal contributions:

\begin{matrix}⟨Φ_{00}^{M} | P V^{c}P QR(0)QP V^{rr}P | Φ_{00}^{M}⟩ = \\ ∑_{k}⟨Φ_{00}^{M} | P V^{c}P QR(0)QP V^{rr}P | Φ_{00}^{M}⟩_{k}+ \\ + ∑_{k≠n}⟨Φ_{00}^{M} | P V^{c}P QR(0)QP V^{rr}P | Φ_{00}^{M}⟩_{kn}. \end{matrix}  (73)

The one-geminal contribution equals to

\begin{matrix}⟨Φ_{00}^{M} | P V^{c}P QR(0)QP V^{rr}P | Φ_{00}^{M}⟩_{k} = ∑_{\begin{matrix}pp′∈R \\ qq′∈R \end{matrix}} ∑_{\begin{matrix}ii′jj′∈\{r,l\} \\ f∈\{a,b\} \end{matrix}} ∑_{σσ′τ} v_{qj_{k}}v_{q′j′_{k}}× \\ ×[  ∑_{ρ_{1}∈ImO_{R}(N_{R}+1)} ^{\frac{F_{pp′qq′}^{σσττ}(ρ_{1})((pp′|i_{k}i′_{k})−δ_{σσ′}(pi′_{k}|i_{k}p′))Ω_{i_{kσ′} i{′ \atop kσ′} j_{kτ} j{′ \atop kτ}}^{f_{k−τ}}}{(I_{kf} −A_{ρ_{1}}−g_{f_{k}ρ_{1}})}} − \\ − ∑_{ρ_{1}∈ImO_{R}(N_{R}+1)} ^{\frac{δ_{σσ′}δ_{στ} F_{pp′qq′}^{σ−σ−σσ}(ρ_{1})(pi′_{k}|i_{k}p′)Ω_{i_{k−σ}i{′ \atop kσ}j_{kσ}j{′ \atop k−σ}}^{f_{kσ}}}{(I_{kf} −A_{ρ_{1}}−g_{f_{k}ρ_{1}})}} − \\ − ∑_{ρ_{1}∈ImO_{R}(N_{R}−1)}(−1)^{δ_{jj′}} ^{\frac{G_{pp′qq′}^{σσττ}(ρ_{1})((pp′|i_{k}i′_{k})−δ_{σσ′}(pi′_{k}|i_{k}p′))Ω_{i_{kσ′} i{′ \atop kσ′}_{j_{kτ}} j{′ \atop kτ}}^{f_{k−τ}}}{(I_{ρ_{1}}−A_{kf} −g_{ρ_{1}f_{k}})}} − \\ − ∑_{ρ_{1}∈ImO_{R}(N_{R}−1)}(−1)^{δ_{jj′}} ^{\frac{δ_{σσ′}δ_{στ} G_{pp′qq′}^{σ−σ−σσ}(ρ_{1})(pi′_{k}|i_{k}p′)Ω_{i_{k−σ}i{′ \atop kσ}_{j_{kσ}}j{′ \atop k−σ}}^{f_{kσ}}}{(I_{ρ_{1}}−A_{kf} −g_{ρ_{1}f_{k}})}} ], \end{matrix}  (74)

while the two-geminal contribution is

\begin{matrix}⟨Φ_{00}^{M} | P V^{c}P QR(0)QP V^{rr}P | Φ_{00}^{M}⟩_{kn} = \\ = ∑_{\begin{matrix}pp′∈R \\ qq′∈R \end{matrix}} ∑_{\begin{matrix}ii′jj′∈\{r,l\} \\ f∈\{a,b\} \end{matrix}} ∑_{στ} ∑_{ρ_{1}∈ImO_{R}(N_{R}+1)} v_{qj_{n}}v_{q′j′_{k}}× \\ ×[^{\frac{F_{pp′qq′}^{σσττ}(ρ_{1})((pp′|i_{k}i′_{n})−δ_{στ} (pi′_{n}|i_{k}p′))Ω_{i_{kτ} i{′ \atop nτ} jnτ j{′ \atop kτ}}^{f_{k−τ}}}{(I_{kf} −A_{ρ_{1}}−g_{f_{k}ρ_{1}})}} − \\ −^{\frac{δ_{στ} F_{pp′qq′}^{σ−σ−σσ}(ρ_{1})(pi′_{n}|i_{k}p′)Ω_{i_{k−σ}i{′ \atop nσ}jnσj{′ \atop k−σ}}^{f_{kσ}}}{(I_{kf} −A_{ρ_{1}}−g_{f_{k}ρ_{1}})}} ]− \\ − ∑_{\begin{matrix}pp′∈R \\ qq′∈R \end{matrix}} ∑_{\begin{matrix}ii′jj′∈\{r,l\} \\ f∈\{a,b\} \end{matrix}} ∑_{στ} ∑_{ρ_{1}∈ImO_{R}(N_{R}−1)} v_{qj_{k}}v_{q′j′_{n}}× \\ ×[\frac{G_{pp′qq′}^{σσττ}(ρ_{1})((pp′|i_{k}i′_{n})−δ_{στ} (pi′_{n}|i_{k}p′))Ω_{i_{kτ} i{′ \atop nτ} j{′ \atop nτ} j_{kτ}}^{f_{n−τ}}}{(I_{ρ_{1}}−A_{nf} −g_{ρ_{1}fn})} − \\ −^{\frac{δ_{στ} G_{pp′qq′}^{σ−σ−σσ}(ρ_{1})(pi′_{n}|i_{k}p′)Ω_{i_{k−σ}i{′ \atop nσ}j{′ \atop nσ}j_{k−σ}}^{f_{n−σ}}}{(I_{ρ_{1}}−A_{nf} −g_{ρ_{1}fn})}} ], \end{matrix}  (75)

where the operator multipliers are:

\begin{matrix}F_{pp′qq′}^{σσ′ττ′} (ρ_{1}) = p_{σ}^{+}p′_{σ′} (1 − |0_{R}⟩ ⟨0_{R}|) q_{τ} |ρ_{1}⟩ ⟨ρ_{1}| q_{τ′}^{′+}, \\ G_{pp′qq′}^{σσ′ττ′} (ρ_{1}) = p_{σ}^{+}p′_{σ′} (1 − |0_{R}⟩ ⟨0_{R}|) q_{τ}^{+} |ρ_{1}⟩ ⟨ρ_{1}| q′_{τ′}, \end{matrix}  (76)

and new quantity which is somewhat intermediate between Π and Ξ are introduced:

Ω_{ii′jj′}^{f} = − ∑_{μ≠0} \frac{⟨Φ_{00}^{M} |i^{+}i′| μ⟩ ⟨μ |j^{+}f^{+}fj′| Φ_{00}^{M}⟩}{ϵ_{μ}}.  (77)

All the excited states, their relative energies and matrix averages are determined above. As for the Ξ’s, the summation over upper indices of the Ω’s gives the polarization propagator. Therefore, assuming that ionization energies for the different BO’s of the k-th geminal are equal we obtain

\begin{matrix}⟨Φ_{00}^{M} | P V^{c}P QR(0)QP V^{rr}P | Φ_{00}^{M}⟩_{k} = ∑_{\begin{matrix}pp′∈R \\ qq′∈R \end{matrix}} ∑_{ii′jj′∈\{r,l\}} ∑_{σσ′τ} v_{qj_{k}}v_{q′j′_{k}}× \\ ×[  ∑_{ρ_{1}∈ImO_{R}(N_{R}+1)} \frac{F_{pp′qq′}^{σσττ}(ρ_{1})((pp′|i_{k}i′_{k})−δ_{σσ′}(pi′_{k}|i_{k}p′))Π_{i_{kσ′} i{′ \atop kσ′} j_{kτ} j{′ \atop kτ}}^{M}}{(I_{k}−A_{ρ_{1}}−g_{kρ_{1}})} − \\ − ∑_{ρ_{1}∈ImO_{R}(N_{R}+1)} ^{\frac{δ_{σσ′}δ_{στ} F_{pp′qq′}^{σ−σ−σσ}(ρ_{1})(pi′_{k}|i_{k}p′)Π_{i_{k−σ}i{′ \atop kσ}j_{kσ}j{′ \atop k−σ}}^{M}}{(I_{k}−A_{ρ_{1}}−g_{kρ_{1}})}} − \\ − ∑_{ρ_{1}∈ImO_{R}(N_{R}−1)}(−1)^{δ_{jj′}} ^{\frac{G_{pp′qq′}^{σσττ}(ρ_{1})((pp′|i_{k}i′_{k})−δ_{σσ′}(pi′_{k}|i_{k}p′))Π_{i_{kσ′} i{′ \atop kσ′}_{j_{kτ}} j{′ \atop kτ}}^{M}}{(I_{ρ_{1}}−A_{k}−g_{ρ_{1}k})}} − \\ − ∑_{ρ_{1}∈ImO_{R}(N_{R}−1)}(−1)^{δ_{jj′}} ^{\frac{δ_{σσ′}δ_{στ} G_{pp′qq′}^{σ−σ−σσ}(ρ_{1})(pi′_{k}|i_{k}p′)Π_{i_{k−σ}i{′ \atop kσ}_{j_{kσ}}j{′ \atop k−σ}}^{M}}{(I_{ρ_{1}}−A_{k}−g_{ρ_{1}k})}} ], \end{matrix}  (78)

for the one-geminal contributions and

\begin{matrix}⟨Φ_{00}^{M} | P V^{c}P QR(0)QP V^{rr}P | Φ_{00}^{M}⟩_{kn} = \\ = ∑_{\begin{matrix}pp′∈R \\ qq′∈R \end{matrix}} ∑_{ii′jj′∈\{r,l\}} ∑_{στ} ∑_{ρ_{1}∈ImO_{R}(N_{R}+1)} v_{qj_{n}}v_{q′j′_{k}}× \\ ×[\frac{F_{pp′qq′}^{σσττ}(ρ_{1})((pp′|i_{k}i′_{n})−δ_{στ} (pi′_{n}|i_{k}p′))Π_{i_{kτ} i{′ \atop nτ} jnτ j{′ \atop kτ}}^{M}}{(I_{k}−A_{ρ_{1}}−g_{kρ_{1}})} − \\ −^{\frac{δ_{στ} F_{pp′qq′}^{σ−σ−σσ}(ρ_{1})(pi′_{n}|i_{k}p′)Π_{i_{k−σ}i{′ \atop nσ}jnσj{′ \atop k−σ}}^{M}}{(I_{k}−A_{ρ_{1}}−g_{kρ_{1}})}} ]− \\ − ∑_{\begin{matrix}pp′∈R \\ qq′∈R \end{matrix}} ∑_{ii′jj′∈\{r,l\}} ∑_{στ} ∑_{ρ_{1}∈ImO_{R}(N_{R}−1)} v_{qj_{k}}v_{q′j′_{n}}× \\ ×[\frac{G_{pp′qq′}^{σσττ}(ρ_{1})((pp′|i_{k}i′_{n})−δ_{στ} (pi′_{n}|i_{k}p′))Π_{i_{kτ} i{′ \atop nτ} j{′ \atop nτ} j_{kτ}}^{M}}{(I_{ρ_{1}}−A_{n}−g_{ρ_{1}n})} − \\ −\frac{δ_{στ} G_{pp′qq′}^{σ−σ−σσ}(ρ_{1})(pi′_{n}|i_{k}p′)Π_{i_{k−σ}i{′ \atop nσ}j{′ \atop nσ}j_{k−σ}}^{M}}{(I_{ρ_{1}}−A_{n}−g_{ρ_{1}n})} ], \end{matrix}  (79)

for the two-geminal contributions.

The cross term between the P V^{c}P and P V^{R}P equals

\begin{matrix}⟨Φ_{00}^{M} | P V^{c}P QR(0)QP V^{R}P | Φ_{00}^{M}⟩ = \\ = − ∑_{pp′∈R} ∑_{B∈R} Z_{B}^{R} ∑_{σσ′τ} p_{σ}^{+}p′_{σ} (1 − |0_{R}⟩ ⟨0_{R}|) × \\ × ∑_{kn∈M} ∑_{ii′jj′∈\{r,l\}}((pp′ | i_{k}i′_{n}) − δ_{σσ′}(pi′_{n} | i_{k}p′))V_{B}^{j_{n}j′_{k}}Π_{i_{kσ′}i′_{nσ′}j_{nτ} j′_{kτ}}^{M} . \end{matrix}  (80)

In the case when the ZDO approximation is employed for the two-electron integrals we have:

\begin{matrix}⟨Φ_{00}^{M} | P V^{c}P QR(0)QP V^{R}P | Φ_{00}^{M}⟩ = \\ = − ∑_{p∈R} ∑_{B∈R} Z_{B}^{R} ∑_{k∈M} ∑_{i∈\{r,l\}} ∑_{σττ′} p_{σ}^{+}p_{σ} (1 − |0_{R}⟩ ⟨0_{R}|) × \\ ×((pp | i_{k}i_{k}) − δ_{σσ′}(pi_{k} | i_{k}p)) (γ_{BA_{k}} − γ_{BB_{k}}) Π_{i_{kσ′}i_{kσ′}r_{kτ} r_{kτ}}^{M} . \end{matrix}  (81)

The cross term between P V^{rr}P and P V^{R}P equals

\begin{matrix}⟨Φ_{00}^{M} | P V^{R}P QR(0)QP V^{rr}P | Φ_{00}^{M}⟩ = \\ = ∑_{k}⟨Φ_{00}^{M} | P V^{R}P QR(0)QP V^{rr}P | Φ_{00}^{M}⟩_{k}+ \\ + ∑_{k≠n}⟨Φ_{00}^{M} | P V^{R}P QR(0)QP V^{rr}P | Φ_{00}^{M}⟩_{kn}. \end{matrix}  (82)

The one-geminal contribution equals

\begin{matrix}⟨Φ_{00}^{M} | P V^{R}P QR(0)QP V^{rr}P | Φ_{00}^{M}⟩_{k} = \\ = − ∑_{pp′∈R} ∑_{στ} ∑_{\begin{matrix}ii′jj′∈\{r,l\} \\ f∈\{a,b\} \end{matrix}} ∑_{B∈R} Z_{B}^{R}V_{B}^{i_{k}i′_{k}}v_{pj_{k}}v_{p′j′_{k}} (1 − |0_{R}⟩ ⟨0_{R}|) × \\ ×[  ∑_{ρ_{1}∈ImO_{R}(N_{R}+1)} ^{\frac{p_{τ} |ρ_{1}⟩⟨ρ_{1}|p_{τ}^{′+}Ω_{i_{kσ}i{′ \atop kσ}j_{kτ} j{′ \atop kτ}}^{f_{k−τ}}}{I_{kf} −A_{ρ_{1}}−g_{f_{k}ρ_{1}}}} − \\ − ∑_{ρ_{1}∈ImO_{R}(N_{R}−1)}(−1)^{δ_{jj′}} ^{\frac{p_{τ}^{+}|ρ_{1}⟩⟨ρ_{1}|p′_{τ} Ω_{i_{kσ}i{′ \atop kσ}_{j_{kτ}} j{′ \atop kτ}}^{f_{k−τ}}}{I_{ρ_{1}}−A_{kf} −g_{ρ_{1}f_{k}}}}], \end{matrix}  (83)

while the two-geminal contribution equals

\begin{matrix}⟨Φ_{00}^{M} | P V^{R}P QR(0)QP V^{rr}P | Φ_{00}^{M}⟩_{kn} = \\ = − ∑_{pp′∈R} ∑_{σ} ∑_{\begin{matrix}ii′jj′∈\{r,l\} \\ f∈\{a,b\} \end{matrix}} ∑_{B∈R} Z_{B}^{R}V_{B}^{i_{k}i′_{n}} (1 − |0_{R}⟩ ⟨0_{R}|) × \\ ×[  ∑_{ρ_{1}∈ImO_{R}(N_{R}+1)} ^{\frac{p_{σ}|ρ_{1}⟩⟨ρ_{1}|p_{σ}^{′+}Ω_{i_{kσ}i{′ \atop nσ}jnσj{′ \atop kσ}}^{f_{k−σ}} v_{pjn}v_{p′j′_{k}}}{I_{kf} −A_{ρ_{1}}−g_{f_{k}ρ_{1}}}} − \\ − ∑_{ρ_{1}∈ImO_{R}(N_{R}−1)}(−1)^{δ_{jj′}} \frac{p_{σ}^{+}|ρ_{1}⟩⟨ρ_{1}|p′_{σ}Ω_{i_{kσ}i{′ \atop nσ}j{′ \atop nσ}j_{kσ}}^{f_{n−σ}} v_{pj_{k}} v_{p′j′_{n}}}{I_{ρ_{1}}−A_{nf} −g_{ρ_{1}fn}}]. \end{matrix}  (84)

This term, in the case when the ZDO approximation is used for two-electron integrals, equals:

\begin{matrix}⟨Φ_{00}^{M} | P V^{R}P QR(0)QP V^{rr}P | Φ_{00}^{M}⟩ = \\ = − ∑_{pq∈R} ∑_{B∈R} Z_{B}^{R} (1 − |0_{R}⟩ ⟨0_{R}|) ∑_{k∈M} ∑_{\begin{matrix}jj′∈\{r,l\} \\ f∈\{a,b\} \end{matrix}} ∑_{στ} v_{pj_{k}}v_{qj′_{k}}× \\ ×[  ∑_{ρ_{1}∈ImO_{R}(N_{R}+1)} ^{\frac{∑_{σ} p_{τ} |ρ_{1}⟩⟨ρ_{1}|q_{τ}^{+}(γ_{BA_{k}} Ω_{r_{kσ}r_{kσ}j_{kτ} j{′ \atop kτ}}^{f_{k−τ}} +γ_{BB_{k}} Ω_{l_{kσ}l_{kσ}j_{kτ} j{′ \atop kτ}}^{f_{k−τ}} )}{I_{kf} −A_{ρ_{1}}−g_{f_{k}ρ_{1}}}}− \\ − ∑_{ρ_{1}∈ImO_{R}(N_{R}−1)}(−1)^{δ_{jj′}} ^{\frac{∑_{σ} p_{τ}^{+}|ρ_{1}⟩⟨ρ_{1}|q_{τ}(γ_{BA_{k}} Ω_{r_{kσ}r_{kσ}_{j_{kτ}} j{′ \atop kτ}}^{f_{k−τ}} +γ_{BB_{k}} Ω_{l_{kσ}l_{kσ}_{j_{kτ}} j{′ \atop kτ}}^{f_{k−τ}} )}{I_{ρ_{1}}−A_{kf} −g_{ρ_{1}f_{k}}}}]. \end{matrix} (85)

If the energies of electron transfer from (or to) the k-th geminal does not depend on the type of ionized orbital we obtain that the Eq. (85) transforms to

\begin{matrix}⟨Φ_{00}^{M} | P V^{R}P QR(ω)QP V^{rr}P | Φ_{00}^{M}⟩ = \\ = − ∑_{pq∈R} ∑_{B∈R} Z_{B}^{R} (1 − |0_{R}⟩ ⟨0_{R}|) ∑_{k∈M} ∑_{jj′∈\{r,l\}} ∑_{στ} v_{pj_{k}}v_{qj′_{k}}× \\ ×[  ∑_{ρ_{1}∈ImO_{R}(N_{R}+1)} ^{\frac{∑_{σ} p_{τ} |ρ_{1}⟩⟨ρ_{1}|q_{τ}^{+}(γ_{BA_{k}} Π_{r_{kσ}r_{kσ}j_{kτ} j{′ \atop kτ}}^{M} +γ_{BB_{k}} Π_{l_{kσ}l_{kσ}j_{kτ} j{′ \atop kτ}}^{M} )}{I_{k}−A_{ρ_{1}}−g_{kρ_{1}}}}− \\ − ∑_{ρ_{1}∈ImO_{R}(N_{R}−1)}(−1)^{δ_{jj′}} ^{\frac{∑_{σ} p_{τ}^{+}|ρ_{1}⟩⟨ρ_{1}|q_{τ}(γ_{BA_{k}} Π_{r_{kσ}r_{kσ}_{j_{kτ}} j{′ \atop kτ}}^{M} +γ_{BB_{k}} Π_{l_{kσ}l_{kσ}_{j_{kτ}} j{′ \atop kτ}}^{M} )}{I_{ρ_{1}}−A_{k}−g_{ρ_{1}k}}}]. \end{matrix} (86)

The last contribution to the effective Hamiltonian of the quantum subsystem is a c-number and can be written as

⟨Φ_{00}^{M} | P V^{R}P | Φ_{00}^{M}⟩ = −2 ∑_{A∈M} ∑_{\begin{matrix}B∈R \\ B≠A \end{matrix}} γ_{AB}Z_{B}^{R} ∑_{\begin{matrix}t_{n}∈A \\ n∈M \end{matrix}} P_{n}^{tt}.  (87)

Thereby, we obtained the explicit form of the effective Hamiltonian for the Rsubsystem. It allows not only to determine the wave function of the R-subsystem but also to obtain the electronic energy of the whole molecule, which can be written [10] as

E_{0}(q) = E_{0}^{R}(q) + E_{00}^{M}(q),

where E_{0}^{R}(q) is the lowest eigenvalue of the effective Hamiltonian Eq. (6) of the Rsubsystem and E_{00}^{M}(q) is the energy of the M-subsystem which is parameterized in the MM form (The detailed transition from the local QM description of molecular electronic structure to the MM will be published elsewhere [16]).

To demonstrate the importance of transition from the bare Hamiltonian for the R-subsystem to its effective Hamiltonian we have estimated using the above formulae the renormalization of the two-electron parameters of Coulomb interaction (γ_{11} and γ_{12}) in the PPP Hamiltonian due to interaction of the π-system with the σ-core. In fact the σ-π-separation is one of the first examples of separation of electron variables. Incidentally, in this case the contribution of the intersubsystem resonance vanishes for the symmetry reasons. We want to emphasize that the values of correction are independent on the bare values of parameters. We start from the bare value of ∼7.6 eV for the γ_{12} parameter for the double C-C bond calculated by the formulae accepted in the standard MINDO/3 method and the bare value of ∼11.8 eV for the γ_{11} parameter which is obtained from atomic spectra [19]. We accepted the bond length value 1.339 Åwhich corresponds to that in ethylene. The correction to the γ_{11} and γ_{12} parameters is due to the contribution Eq. (37) to the effective Hamiltonian only. In fact in the ethylene molecule the geminals corresponding to the C-C and C-H bonds contribute to the correction. The numeric estimate of the correction to γ_{11} is negative (the renormalized γ_{11} is smaller than its bare value) and equals to ∼0.45 eV. This value is not very large which justifies the application of the perturbation theory to correct the bare Hamiltonian. At the same time this value suffices to be important for constructing the Hamiltonian for the π-subsystem. Also we obtained that the renormalization of the γ_{12} parameter is very small (less than 0.3%) and thus can be neglected.

3. Discussion.

The QM/MM methods become more and more popular in theoretical modeling of large molecular systems. It is caused by rapidly growing needs of biological and organometallic chemistry. The application of the QM/MM schemes to calculations on molecular electronic structure and chemical reactivity is very advantageous. However, the results of up-to-day QM/MM calculations depend on the form of the junction between QM and MM subsystems and, therefore, are ambiguous. All the junction forms proposed in the literature are empirical. It seems to be important to substantiate the form of the junction on the basis of quantum chemical consideration of interaction between the parts of the molecular system treated on the quantum and classical level. In the present work we performed a sequential derivation of the hybrid QM/MM scheme on the basis of a previously developed local description of molecular electronic structure [14]. This allows to express finally the renormalizations of the bare Hamiltonian in terms of the local quantities characterizing the two-electron bonds in the chemically inert i.e. classical subsystem, which can be either calculated or estimated from experimental data on the bond polarizabilities [18].

In order to construct the effective Hamiltonian for the quantum system or equivalently the QM/MM junction we have taken into account both electron transfers between the subsystems and the excitations in the classical part of the molecular system. We have obtained the explicit expressions for this renormalization, which is the function of well-defined characteristics like the polarization propagators of the M-subsystem. Our consideration also required to introduce new quantities Ξ and Ω which appear due to entanglements of one-electron transfers between the subsystems and excitations leaving the number of electrons in two subsystems unchanged. If the energies of the states with one electron transfers between the two subsystems are replaced with an average value, the quantities Ξ and Ω in the expressions for the effective Hamiltonian for the quantum subsystem are reduced to the polarization propagators. This allows to express the contributions renormalizing the Coulomb interaction of electrons in the quantum part of the system through such observable quantities as the bond polarizabilities are. The same applies to the corrections renormalizing the one-electron terms of the Hamiltonian for the quantum part. The possibility to express the corrections to the Hamiltonian for the R-subsystem through the experimentally observable quantities such as bond polarizabilities and the ionization potentials of the bonds allows to eliminate the calculations of the electronic structure of the classical part of the whole molecular system and to parameterize them. The above reduction of the QM/MM junction to the sum of transferable contributions from the chemical bonds remaining in the classical part was possible due to the special form for the implicit wave function of electrons underlying the description of the classical part of the system with use of the MM type schemes. We use the APSLG form for this wave function. This provides a route to a possible derivation of the classical (MM) description of the M-subsystem.

With use of the formulae derived above for the effective Hamiltonian for the quantum subsystem we addressed the old problem of substantiation of the parameters of the PPP Hamiltonian (see, for example, [20]). The PPP Hamiltonian for the π-electrons is historically one of the earliest examples of successful separation of electron variables with those in the π-subsystem treated on the QM level, while those in the σ-subsystem left to a classical (in fact, MM-like) description. The value of parameter γ_{11} in the PPP Hamiltonian equals approximately to 11.1 eV (the numerical values of parameters vary by about 0.5 eV in different works on the π-Hamiltonian). The value obtained from atomic spectra of carbon [19] is larger by ∼0.7 eV. This value is effective and implicitly includes the multiparticle corrections [19,21] which appear due to reduction of the complete one-electron basis to the valence one. The difference between the PPP and experimentally estimated atomic values can be explained by the influence of the σ-core on the π-subsystem. The correction value (0.45 eV) obtained using the formulae written above and the MINDO/3 approximation for the molecular integrals is smaller than the required value of 0.7 eV by absolute value but the sign (negative) and the order of magnitude for this increment are correct. The bare value of the γ_{12} parameter strongly depends on the type of approximation. The correction to this value is small (see above). In Ref. [21] an attempt to obtain the corrections to π-Hamiltonian of benzene from the σ-core is mentioned. The result was a very small correction to the γ_{12} parameter which is in accordance with our result. However, the correction to the γ_{11} parameter is positive and exceeds 1.2 eV. Thus the value of the γ_{11} parameter obtained in that work is by 10% larger than the PPP one. Moreover, the positive sign of the correction contradicts the physical picture of screening of electronic interactions in the π-subsystem due to interaction with the polarization fluctuations in the σ-subsystem (see, for example, [22]). More advanced attempt to obtain parameters of the π-electron Hamiltonian was made in Ref. [23] by using the ab initio effective valence shell Hamiltonian formalism. Starting from the ab initio Hamiltonian the parameters of the π-Hamiltonian were estimated perturbatively. The value of the parameter γ_{12} was obtained to be ∼7.8 eV that is close to the standard value of the PPP Hamiltonian. However, also in Ref. [23] the value of the parameter γ_{11} was obtained to be ∼12 eV. This is significantly larger than the value of the PPP Hamiltonian (which is sufficiently well established from the experimental data) and is close to the value obtained from analysis of the excitation spectra of free carbon atoms and ions [19]. Therefore, we conclude that the screening of electronic interactions due to σ-core can be estimated with use of the formalism introduced in the present paper. It noticeably reduces the value of the one-center electron-electron repulsion γ_{11} parameter obtained from the atomic spectra, but the correction to the two-center electron-electron repulsion γ_{12} parameter is negligibly small.

References

[1] J.A. Pople, Angew. Chem. Int. Ed. 38 (1999) 1894.

[2] F. Bernardi, M. Olivucci, and M.A. Robb, J. Am. Chem. Soc. 114 (1992) 1606.

[3] K. Morokuma, R.D.J. Froese, S. Humbel, M. Svensson, T. Matsubara, and S. Sieber, J. Comp. Chem. 16 (1995) 1170; K. Morokuma, R.D.J. Froese, S. Humbel, M. Svensson, T. Matsubara, and S. Sieber, J. Phys. Chem. 100 (1996) 2573.

[4] V.V. Vasilyev, A.A. Bliznyuk, and A.A. Voityuk, Int. J. Quant. Chem. 44 (1992) 897.

[5] V. Théry, D. Rinaldi, J.-L. Rivail, B. Maigret, and G. Ferenczy, J. Comput. Chem. 15 (1994) 269.

[6] I.B. Bersuker, M.K. Leong, J.E. Boggs, and R.S. Pearlman, Int. J. Quantum Chem. 63 (1997) 1051.

[7] CECAM-NSF Meeting on QC/MM methods, Int. J. Quant. Chem 60 (1996) No 6, Special Issue.

[8] T.N. Truong, E.V. Stefanovich, Chem. Phys. Lett. 256 (1996) 348.

[9] A.L. Tchougréeff, Khim. Fiz. 16 (1997) No 6, 62 (in Russian); Chem. Phys. Reps. 16 (1997) 1035 (in English).

[10] A.L. Tchougréeff, Phys. Chem. Chem. Phys. 1 (1999) 1051.

[11] P.-O. Löwdin, J. Math. Phys. 3 (1962) 969; P.-O. Löwdin, in Perturbation Theory and its Applications in Quantum Mechanics, ed. C.H. Wilcox, J. Wiley, NY, 1966.

[12] R.McWeeny, Methods of Molecular Quantum Mechanics, 2nd Edition, AP, London, 1992.

[13] M.B. Darkhovskii, A.L.Tchougréeff, Khim. Fiz. 18 (1999) No 1, 73 (in Russian).

[14] A.M. Tokmachev, A.L. Tchougréeff, Zh. Fiz. Khim. 73 (1999) 347 (in Russian); Russ. J. Phys. Chem. 73 (1999) 320 (in English); A.M. Tokmachev, A.L. Tchougréeff, J. Comp. Chem., submitted.

[15] P.R. Surján, The Concept of the Chemical Bond, in Theoretical Models of Chemical Bonding, Part 2. Ed. Z.B. Maksić, Springer, Heidelberg, 1989, p. 205.

[16] A.M. Tokmachev, A.L. Tchougréeff, in preparation.

[17] P. Claverie, in Intermolecular Interactions: From Diatomics to Biopolymers, ed. B. Pullman, Wiley-Interscience, NY, 1978.

[18] R.J.W. Le Fevre, Molecular Refractivity and Polarizability, in Advances in Physical Organic Chemistry, Vol. 3, Ed. V. Gold, Academic Press, London and New York, 1965, p. 282.

[19] I.A. Misurkin, A.A. Ovchinnikov. Opt. and Spectr. 30 (1971) 616 (in Russian).

[20] T.E. Peacock, Electronic Properties of Aromatic and Heterocyclic Molecules, Academic Press, London and New York, 1965.

[21] K.F. Freed, Theoretical Basis for Semiempirical Theories, in Semiempirical Methods of Electronic Structure Calculation, Part A. Techniques Ed. G.A. Segal, Plenum Press, New York and London, 1977.

[22] J. Linderberg, Y. Öhrn, Propagators in Quantum Chemistry, Academic Press, London and New York, 1973.

[23] C.H. Martin, K.F. Freed, J. Chem. Phys. 101 (1994) 4011.