Potential Energy Surfaces in Hybrid QM/MM Methods

Content-type: text/HTML

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

Abstract

The problem of conjunction between quantum and classical parts in the hybrid QM/MM methods is considered. The form of the junction is deduced with use of the perturbation expansion in the assumption that the wave function underlying molecular mechanics is of the APSLG type. The renormalization of the potential energy surfaces of the combined system due to interaction between quantum and classical subsystems is represented as a sum of contributions with transparent physical sense.

Key words: hybrid QM/MM methods, perturbation expansion, effective Hamiltonian, potential energy surfaces, strictly localized geminals.

Introduction.

The hybrid quantum mechanical/molecular mechanical (QM/MM) methods became very popular in modern research on large molecular systems (particularly, biological ones) []. These approaches are based on the well justified experimental observation: chemical transformations are usually local, i.e. restricted to a small region of the entire molecular system called reaction center (RC). For that reason it seems to be practical to use approximations of different level in order to describe different parts of the system when its potential energy surface (PES) is considered. The smaller part of the molecular system with actually transforming electronic structure (RC) must be described by a thorough (preferrably, highly correlated) quantum chemical approach while the rest of the molecular system can be described either by molecular mechanics or by a faster quantum chemical method capable to reproduce only general features of the PES. The physical reason of that difference is the relative importance of electron correlations for correct description of the processes of the bond formation and bond cleavage - those where the uncorrelated SCF approach gives wrong asymptotic behavior of the electronic wave function of the system.

Two important questions immediately arise for such hybrid schemes. The first one is about the margin between the QM and MM subsystems. The important restriction must be set on the subsystems. Elaborated quantum chemical approaches (as well as the molecular mechanical schemes) work properly only being applied to the objects with an integer number of electrons. Therefore, distributing electrons between the subsystems is restricted by the condition of small fluctuations of numbers of electrons in the subsystems.

The second question arises from the fact that the classical part of the molecule modifies the PES of the reactive (quantum) subsystem. The interaction between subsystems can not be neglected. It leads to various forms of the junction between quantum and classical subsystems accepted in many hybrid schemes implemented by now in the literature Aqvist,Bak,Carmer,Field,Field1,Gao,Gao1,Gao2,Luzhkov,Maseras,Stanton,Thompson,Tunon,Truong,Morokuma,Bernardi. At the same time it is noteworthy that the form of the junction in these approaches is usually taken ad hoc, without proper substantiation. For example, authors of Ref. [] represent the junction between subsystems as a sum of electrostatic and van der Waals interactions with adjusted parameters; in Ref. [] ''junction dummy atoms'' are introduced; in Ref. [] the intersystem Coulomb and exchange integrals are taken as linear combinations of exponential functions with subsequent parameterization.

From out point of view the correct form of the junction between the quantum and classical parts of the molecular system must be constructed on the basis of the quantum chemical description of the whole system and of the sequential separation of electronic variables related to respective subsystems. This approach was originally proposed in Ref. []. It is based on the perturbation expansion and uses effective electron Hamiltonian technique. The margin between quantum and classical subsystems is chosen so that the chemical bonds are not broken and can be assigned only to one subsystem. The construction of hybrid schemes in Ref. [] supposed that the parameters of the MM method must be renormalized due to interaction with quantum subsystem. Another construction of hybrid schemes was proposed in Ref. []: MM method remains constant and interactions modify only Hamiltonian of quantum subsystem. The effective Hamiltonian of quantum subsystem was explicitly constructed in Ref. TheoChem. The main idea of all these works is to single out the subset of electronic variables responsible for a chemical transformation to be described with use of a QM approach and to construct (using the Löwdin partition technique []) the effective Hamiltonian for the reactive subsystem which takes into account the interactions between subsystems perturbatively. To perform this program i.e. to construct the correction to the bare Hamiltonian for the RC it is necessary to assume some form of the electronic wave function of the inert environment underlying its MM description. This function has been chosen to be of the form of antisymmetrized product of strictly localized geminals (APSLG) Surjan,Weeny,ZhFizKhim,JCC. The APSLG based approach operates by local quantities such as chemical bonds and electron lone pairs, describes the ground states of the organic molecules with closed electron shells and high excitation energies, and provides the natural starting point for a transition from the quantum mechanics to the molecular mechanics (to be described elsewhere []). The effective Hamiltonian of the reactive subsystem in the presence of ''classical'' subsystem described by the APSLG method has been explicitly constructed in Ref. []. The present work is aimed to construct the PES of the molecular system described by the hybrid method with the APSLG wave function underlying the MM description of the classical subsystem.

PES of composite quantum-classical system.

Let us consider the application of the effective electron Hamiltonian method to the separation of the electronic variables which is aimed to constructing the PES for molecular systems in the hybrid QM/MM methods. We denote the reactive subsystem (described by quantum mechanics) by subscript R and the inert subsystem (described by MM) by subscript M. The Hamiltonian for the whole molecular system can be divided into parts corresponding to separate subsystems and to the interaction between them:

H = H^R(q)+H^M(q)+V^c(q)+V^r(q)+ e²
2
å
\Sb A ¹ BA Î R B Î M\endSb Z_A^RZ_B^MR_AB^-1.
(0)
The distinction between the R- and M-subsystems is formalized by ascribing each one-electron basis state of the entire system either to the R - or to the M-subsystem. The operators H^R(q) and H^M(q) are then polynomials with respect to the fermion operators creating or annihilating electrons in the one-electron states ascribed to the respective subsystems. On the other hand the operators containing the products of the fermion operators ascribed to the different subsystems are the interaction operators. We restrict ourselves by only the Coulomb V^c and the resonance V^r (electron transfer) interactions. The core charges are also subdivided between two subsystems. The contributions to the Hamiltonian can be further specified. The Hamiltonian for the R-subsystem is a sum of the Hamiltonian for the free R-subsystem (without environment) H₀^R and of the operator V^M of attraction of electrons of the R-subsystem to the cores of the M-subsystem. Analogously, the Hamiltonian for the M-subsystem H^M can be presented as a sum of the Hamiltonian for free M-subsystem H₀^M and of the attraction of electrons of the M-subsystem to the cores of the R-subsystem V^R.

The wave function for the whole system is then approximated by the antisymmetrized product of the wave function of the reactants (which can be taken as a linear combination of configurations) F_k^R (where k is the number of the electronic state of the R-subsystem) and the ground state wave function of the inert environment to be calculated without reactants F₀₀^M:

Y_k^appr = F_k^RÙF₀₀^M.
(0)

The wave function F₀₀^M is the approximation to the lowest eigenstate of the Hamiltonian H₀^M. It is assumed to be obtained in the framework of the APSLG method [,,,]. Therefore, it has the form:

| F₀₀^M ñ =
Õ
m
g_m⁺|0 ñ ,
(0)
where each geminal (two-electron wave function) g_m⁺|0 ñ is assigned either to a chemical bond or to an electron lone pair.

To construct these geminals a unitary transformation must be performed from the initial atomic basis set to the hybrid one for each heavy (non hydrogen) atom:

t_ms =
å
i Î { s,x,y,z}
h_mi^Aa_is,
(0)
where | t_ms ñ is a spin-orbital with the spin projection s, assigned to the m-th bond and located on the right or the left end of the bond (| r_ms ñ or |l_ms ñ ); h^A Î SO(4) is a matrix of transforming the AO basis orbitals s, p_x, p_y, and p_z of the atom A to the hybrid orbitals(HO's), and a_is are operators annihilating an electron on a spin-orbital of the AO basis (the HO basis set is thus orthonormal).

Each geminal is a linear combination with varying amplitudes of three possible singlet states (two ionic and one covalent - Heitler-London type) constructed from the hybrid spin-orbitals | r_ms ñ and | l_ms ñ :

g_m⁺ = u_mr_ma⁺r_mb⁺+v_ml_ma⁺l_mb⁺+w_m(r_ma⁺l_mb⁺+l_ma⁺r_mb⁺)
(0)
with the normalization condition u_m²+v_m²+2w_m² = 1. These functions are introduced in the theory by S. Weinbaum [].

The approximate wave function Eq. (2) can be obtained from the exact wave function of the whole system

Y_k^exact =
å
n_Mn_R

å
i_Mi_R
C_{i_Mi_R}^k(n_Mn_R)F_{i_M0}^M(n_M)ÙF_{i_R}^R(n_R), (n_M+n_R = N_e)
(0)
by two subsequent projections: first one projects the space where the exact wave function Eq. (6) resides onto the space of functions with fixed number of electrons in each of the subsystems (provided that n_M+n_R = N_e) eliminating by this the fluctuations of numbers of electrons in the subsystems. The corresponding projection operator is denoted as P and Q = 1-P is the corresponding complementary projection operator. The second projection takes into account the fact that the chemically inert environment is described by its ground state wave function (it is also true when the wave function underlying the MM is used). Therefore, the second projection is onto the space of states with the ground state of the M-subsystem taken as a multiplier. The corresponding projection operator P is 1_RÄ| F₀₀^M ñ á F₀₀^M| . Analogously we introduce the complementary projection operator Q = 1-P.

The projections imply the transition from the exact Hamiltonian for the entire molecular system to the relevant effective Hamiltonians. The construction of the effective Hamiltonian is performed to fulfill the requirement that the eigenvalues of the effective Hamiltonian acting in the restricted subspace coincide with those of the original Hamiltonian acting in the full space. First projection yields the effective Hamiltonian of the form:

H_eff(q,E) = PH^R(q)P+PH^M(q)P+PV^c(q)P+

+PV^rr(q,E)P+ e²
2
å
\Sb A ¹ B

A Î R,B Î M\endSbZ_A^RZ_B^MR_AB^-1,

(0)
where

V^rr(q,E) = V^r(q)QR(q,E)QV^r(q)
(0)
with the resolvent operator

R(q,E) = (E-QHQ)^-1.
(0)
The second projection gives the operator acting in the space of functions of the type Eq. (2):

H_eff(q,E,w) = e²
2
å
\Sb A ¹ B

A Î R,B Î M\endSb Z_A^RZ_B^MR_AB^-1+ PP[H₀^R(q)+H₀^M(q)+

+W(q,E)+V^M(q)+W(q,E)PQR(q,w)QPW(q,E)]PP,

(0)
where

W(q,E) = V^c(q)+V^rr(q,E)+V^R(q)
(0)
and the resolvent operator

Â(q,w) = (w-QHQ)^-1.
(0)

The energy of the k-th state of the combined system can be obtain by averaging the effective Hamiltonian over the approximate wave function Eq. (2):

E_k = á Y_k^appr| H_eff(q,E,w)| Y_k^appr ñ .
(0)
This energy obviously differs from that obtained by using ab initio Hamiltonian for R-subsystem and unrenormalized MM description for M-subsystem:

E_k = (E_k0^R+E₀^M)+DE_k,
(0)
where the quantities E_k0^R and E₀^M are the averages of the free subsystem Hamiltonians over the functions F_k^R and F₀₀^M, respectively. The correction to the ab initio energy (junction between quantum and classical parts) arises from interaction between the subsystems:

DE_k = e²
2
å
\Sb A ¹ B

A Î R,B Î M\endSbZ_A^RZ_B^MR_AB^-1+

+ á Y_k^appr| W(q,E)+V^M(q)+W(q,E)PQR(q,w)QPW(q,E)| Y_k^appr ñ .

(0)

To perform the averaging the explicit form of the operators of the Coulomb and resonance interactions is needed. The operator of the Coulomb interaction is taken to conserve the number of electrons in the subsystems and thus can be written as:

V^c(q) = å
\Sb pp^¢ Î R mm^¢ Î M\endSb
å
st
(p_s p_s ^¢||m_t m_t ^¢)p_s ⁺m_t ⁺m_t ^¢p_s^¢,
(0)
where

(p_s p_s ^¢||m_t m_t^¢) = (pp^¢ | mm^¢)-d_st(pm^¢ | mp^¢)
(0)
and the indici pp^¢ and mm^¢ refer to the spatial parts of one-electron states in the R- and M-subsystems, 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 ZDO approximation results in its turn in further simplification due to the restrictions: p = p^¢ and m = m^¢.

The operator of resonance interaction describes the one-electron transfers from the M-subsystem to the R-subsystem or vice versa:

V^r(q) = å
\Sb p Î R m Î M\endSbv_pm(q)
å
s
(p_s ⁺m_s +m_s ⁺p_s ).
(0)
It is convenient to introduce one-electron density matrices for each subsystem. The elements of this matrix for the R-subsystem are:

P_k^{spp^¢} = á F_k^R|p_s ⁺p_s ^¢| F_k^R ñ .
(0)
The APSLG structure of the M-subsystem assumes that the one-electron densities differ from zero only for orbitals belonging to the same geminal. Also, the one-electron density matrix of the M-subsystem is spin-independent [,]. Therefore the elements of the one-electron density matrix for the M-subsystem are:

P_m^{tt^¢} = á F₀₀^M| t_ms⁺t_ms^¢| F₀₀^M ñ = á0| g_mt_ms⁺t_ms^¢g_m⁺|0 ñ ;

P_m^rr = u_m²+w_m²;P_m^ll = v_m²+w_m²;P_m^rl = P_m^lr = (u_m+v_m)w_m.

(0)
It is also convenient to introduce reduced Coulomb integrals:

Y_pp^¢mm^¢ =
å
t
(p_sp_s ^¢||m_t m_t ^¢) = 2(pp^¢ | mm^¢)-(pm^¢ | mp^¢).
(0)

The Coulomb electron interaction between the subsystems contributes to the energy:

á Y_k^appr| V^c(q)| Y_k^appr ñ =
å
pp^¢ Î R

å
t_mt_m^¢ Î M
Y_{pp^¢i_mj_m}P_m^{tt^¢}
å
s
P_k^{spp^¢}.
(0)
If the ZDO approximation is used and orbitals are attached to different atoms (A ¹ B), this contribution to the energy can be rewritten in the form:

2 å
\Sb A Î R B Î M\endSb g_AB æ
è
å
t_m Î B
P_m^tt ö
ø æ
è
å
p Î A

å
s
P_k^spp ö
ø .
(0)
In the case of frontier atoms the integral g_AB must be replaced by Y_{ppt_mt_m} which can not be factorized from the sums over t_m and p , respectively. The operator of attraction of electrons in the R-subsystem to the cores of the M-subsystem V^M can be written as

V^M(q) = -e²\stackunderB^¢ å
Z_B^¢^M
| r-R_B^¢|
= -
å
pp^¢ Î R

å
B Î M
V_B^{pp^¢}Z_B^M
å
s
p_s⁺p_s ^¢,

V_B^{pp^¢} = -e² ó
õ d³r y_p^*(r)y_p^¢(r)
| r-R_B^¢|
.

(0)
The averaging of this expression over the approximate wave function gives:

á Y_k^appr| V^M(q)| Y_k^appr ñ = -
å
pp^¢ Î R

å
B Î M
V_B^{pp^¢}Z_B^M
å
s
P_k^{spp^¢}.
(0)
In the case of the ZDO approximation the p = p^¢ condition applies. Analogously, the contribution to the energy from the attraction of electrons in the M-subsystem to the cores of the R-subsystem is:

á Y_k^appr| V^R(q)| Y_k^appr ñ = -2
å
t_mt_m^¢ Î M

å
A Î R
V_A^{t_mt_m^¢}Z_A^RP_m^{tt^¢}.
(0)
It can be noted that the contributions Eqs. (23), (25), and (26) can be combined with the repulsion of the cores of the R- and M-subsystems yielding the total electrostatic contribution to the renormalization of the PES. Each of these four contributions is large but the total electrostatic correction is relatively small:

DE_k^el-st = e²
2
å
\Sb A ¹ B

A Î R,B Î M\endSbZ_A^RZ_B^MR_AB^-1+

+ á Y_k^appr| V^c(q)+V^M(q)+V^R(q)|Y_k^appr ñ .

(0)
If the ZDO approximation is used the electrostatic contribution from different atoms of the R- and M-subsystems acquires the familiar form

å
A Î R

å
B Î M
g_ABQ_A^RkQ_B^M,
(0)
where the effective atomic charges are defined by:

Q_A^Rk =
å
p Î R

å
s
P_k^spp-Z_A^R,

Q_B^M = 2
å
t_m Î B
P_m^tt-Z_B^M.

(0)
The contribution to the PES from the operator V^rr(q,E) corresponds to the second order of the perturbation theory with respect to the intersystem one-electron transfers. To estimate it we present the corresponding resolvent operator in the approximate form:

R(E) =
å
i Î ImQ
| i ñ á i|
E-E_i
,
(0)
where | i ñ are the states with one electron transferred from the M-subsystem to the R-subsystem or vice versa. Each state | i ñ is an antisymmetrized product of the ionized states | m ñ and | r ñ of the M- and R-subsystems respectively. The energy denominators in the resolvent expression Eq. (30) are expressed through the ionization potentials (IP) and the electron affinities (EA) of the respective subsystems:

E_i-E = ì
í
î

I_m -A_r -g_mr

I_r -A_m -g_rm
.
(0)
The approximate form of ionized states of the M-subsystem as it appears from the APSLG approximation is the following. The ionized states |m ñ are assumed to be obtained by removing or adding an electron from or to bonding and antibonding orbitals of the geminals of the M-subsystem. These orbitals can be written as []:

ì
í
î

b_ms = x_ml_ms+y_mr_ms

a_ms = -y_ml_ms+x_mr_ms
; (x_m²+y_m² = 1),
(0)
where x_m and y_m are chosen to represent the m-th geminal in the form:

g_m⁺ = U_mb_ma⁺b_mb⁺+V_ma_ma⁺a_mb⁺; (U_m²+V_m² = 1).
(0)
The coefficient sets (U_m,V_m;x_m,y_m) and (u_m,w_m,v_m) are uniquely related:

ì
ï
í
ï
î

u_m = U_my_m²+V_mx_m²

v_m = U_mx_m²+V_my_m²

w_m = ( U_m-V_m) x_my_m

(0)
and the APSLG-structure of the wave function can be expressed both in the basis of hybrid orbitals and in the basis of bond orbitals. The approximate wave functions of the ionized states of the M-subsystem have the form:

b_ms⁺
Õ
l ¹ m
g_l⁺ | 0ñ, a_ms⁺
Õ
l ¹ m
g_l⁺ | 0ñ,

b_ms⁺r_m-s⁺l_m-s⁺
Õ
l ¹ m
g_l⁺ | 0ñ, a_ms⁺r_m-s⁺l_m-s⁺
Õ
l ¹ m
g_l⁺ | 0ñ.

(0)
The corresponding IP's and the EA's entering the Eq. (31) are estimated in Ref. []. It is convenient to introduce some vacuum averages for the M-subsystem:

q_if^m = á 0| g_mi_maf_mb⁺r_ma⁺l_ma⁺| 0 ñ ,

h_if^m = á 0| g_mi_ma⁺f_mb⁺|0 ñ .

(0)
These averages can be readily evaluated:

q_rb^m = y_mw_m+x_mv_m, q_ra^m = x_mw_m-y_mv_m;

q_lb^m = -x_mw_m-y_mu_m, q_la^m = y_mw_m-x_mu_m;

h_rb^m = y_mu_m+x_mw_m, h_ra^m = x_mu_m-y_mw_m;

h_lb^m = y_mw_m+x_mv_m, h_la^m = x_mw_m-y_mv_m.

(0)

Using these notations we obtain the contribution to the energy from the operator V^rr(q,E):

DE_k^rr = á Y_k^appr| V^rr(q,E)| Y_k^appr ñ = -
å
pp^¢ Î R

å
m Î M

å
ij Î { r,l}

å
f Î {a,b}
v_{pi_m}v_{p^¢j_m}×

×
å
s
(
å
r Î ImO_R(N_R-1)
á F_k^R| p_s ⁺| r ñ á r| p_s ^¢| F_k^R ñq_if^mq_jf^m
I_kr-A_mf-g_{krf_m}
+

+
å
r Î ImO_R(N_R+1)
á F_k^R|p_s | r ñ á r| p_s^¢+| F_k^R ñ h_if^mh_jf^m
I_mf-A_kr-g_{f_mkr}
).

(0)
This expression can be simplified if the IP's and EA's are replaced by their average values:

DE_k^rr » -
å
pp^¢ Î R

å
m Î M

å
ij Î { r,l}

å
f Î {a,b}
v_{pi_m}v_{p^¢j_m}×

×
å
s
æ
ç
è P_k^{spp^¢}q_if^mq_jf^m
I_k-A_mf-g_{kf_m}
+ ( 1-P_k^{sp^¢p}) h_if^mh_jf^m
I_mf-A_k-g_{f_mk}
ö
÷
ø .

(0)
The contribution from the second projection can be presented as a sum:

áY_k^appr | W(q,E)PQR(q,w)QPW(q,E) | Y_k^apprñ =

= áY_k^appr | V^c(q)P QR(q,w)QPV^c(q) | Y_k^apprñ+

+áY_k^appr | V^R(q)PQR(q,w)QPV^R(q) | Y_k^apprñ+

+2áY_k^appr | V^c(q)P QR(q,w)QPV^R(q) | Y_k^apprñ+

+2áY_k^appr | V^c(q)P QR(q,w)QPV^rr(q,E) | Y_k^apprñ+

+2áY_k^appr | V^rr(q,E)P QR(q,w)QPV^R(q) | Y_k^apprñ+

+áY_k^appr | V^rr(q,E)PQR(q,w)QPV^rr(q,E) | Y_k^apprñ.

(0)
The first three contributions of Eq. (40) are large enough but the perturbation expansion is valid because these contributions are grouped into the relatively small correction. At the same time, for convenience we consider these contributions separately. The resolvent can be approximated by the sum over the states:

Â(q,w) = å
\Sb rm ¹ 0\endSb | r,m ñ á r,m|
w-w_rm
.
(0)
We assume that the dependence of the resolvent on w is weak and the w values of interest are much smaller than the resolvent poles which are all higher than 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 pass to the limit w® 0 in the above resolvent. The first contribution to the expression Eq. (40) is the second order correction due to the Coulomb repulsion operator. To write down this contribution it is convenient to introduce dynamic polarization propagators [] for the R- and M-subsystems:

P_pp^¢qq^¢^Rk(w) =
å
r ¹ k
(e_r -e_k) á F_k^R|p⁺p^¢| r ñ á r|q⁺q^¢| F_k^R ñ
(e_r -e_k)²-w²
,

P_mm^¢nn^¢^M(w) =
å
m ¹ 0
e_m á F₀₀^M| m⁺m^¢| m ñ á m| n⁺n^¢| F₀₀^M ñ
e_m²-w²
,

(0)
where e_r and e_m are the energies of excitations in the R- and M-subsystems, respectively. The excitation energies and matrix elements necessary for calculating the polarization propagator of the M-subsystem are given in Ref. []. Using those notations we obtain:

áY_k^appr | V^c(q)P QR(q,0)QPV^c(q) | Y_k^apprñ =

=
å
pp^¢qq^¢ Î R

å
i_mi_n^¢j_nj_m^¢ Î M
{ å
\Sb st

s^¢t^¢\endSb (p_s p_s ^¢||i_ms^¢i_ns^¢^¢)(q_tq_t ^¢||j_nt^¢j_mt^¢^¢)×

× é
ê
ë P_k^{spp^¢}P_k^{tqq^¢}P_{i_ms^¢i_ns^¢^¢j_nt^¢j_mt^¢^¢}^M(0)+ 2
p
ó
õ P_{p_sp_s ^¢q_t q_t ^¢}^Rk(iu)P_{i_ms^¢i_ns^¢^¢j_nt^¢j_mt^¢^¢}^M(iu)du ù
ú
û .

(0)
The expression using integrals is not useful in practice. At the same time we can approximate it taking into account the essentially higher energies of excitations in the M-subsystem in comparison with those of the R-subsystem:

2
p
ó
õ P_{p_s p_s ^¢q_t q_t ^¢}^Rk(iu)P_{i_ms^¢i_ns^¢^¢j_nt^¢j_mt^¢^¢}^M(iu)du »

» P_{i_ms^¢i_ns^¢^¢j_nt^¢j_mt^¢^¢}^M(0)
å
k^¢ ¹ k
á F_k^R| p_s ⁺p_s ^¢| F_k^¢^R ñ á F_k^¢^R|q_t ⁺q_t ^¢| F_k^R ñ ,

(0)
i.e. it is a sum of products of the static polarization propagator of M-subsystem and transition electron densities of R-subsystem. The contribution Eq. (44) mainly corresponds to the dispersion interaction between the subsystems DE_k^disp.

The contribution to the energy of the second order in operators V^R(q) does not change the relative energies of the states in the R-subsystem and, therefore, is k-independent, where k is the number of electronic state in the R-subsystem. The operator V^R(q) can be written in the form:

V^R = -
å
mm^¢ Î M

å
s
m_s⁺m_s ^¢
å
A Î R
V_A^{mm^¢}Z_A^R

» - å
\Sb B Î M\endSb
å
m Î BÇM

å
s
m_s ⁺m_s å
\Sb A Î R

A ¹ B\endSb g_ABZ_A^R.

(0)
The contribution to the energy is:

DE_k^RR = áY_k^appr | V^R(q)PQR(q,0)QPV^R(q) | Y_k^apprñ =

=
å
i_mi_n^¢j_nj_m^¢ Î M

å
AA^¢ Î R
V_A^{i_mi_n^¢}V_A^¢^{j_nj_m^¢}Z_A^RZ_A^¢^R
å
st
P_{i_msi_ns^¢j_ntj_mt^¢}^M(0).

(0)
Let us denote the atom at the right end of the m-th bond as R_m and that at the left end as L_m. In the case if AA^¢ Ï {R_m,L_m} the contribution Eq. (46) can be approximately rewritten in the physically transparent form in terms of the bond polarizability tensors [^(a)]^(m)(w) []:

-
å
m Î M
å
\Sb AA^¢ Î R AA^¢ Ï { R_m,L_m} \endSb Z_A^RZ_A^¢^R(ÑV(R_A^¢^(m)) ê
ê ^
a

^(m)(0) ê
ê ÑV(R_A^(m))),
(0)
where V(R_A^(m)) is the potential induced by a unit charge placed on the atom A at the center of the m-th bond. It is worth to mention that the bond polarizabilities [^(a)] ^(k)(0) are tabulated, for example, in Ref. [].

The cross term between the V^c(q) and V^R(q) equals to:

DE_k^cR = áY_k^appr | V^c(q)P QR(q,0)QPV^R(Q) | Y_k^apprñ =

= -
å
ss^¢t

å
pp^¢ Î R
P_k^{spp^¢}
å
A Î R
Z_A^R
å
i_mi_n^¢j_nj_m^¢ Î M
(p_s p_s^¢||i_ms^¢i_ns^¢^¢)V_A^{j_nj_m^¢}P_{i_ms^¢i_ns^¢^¢j_ntj_mt^¢}^M(0).

(0)
The contributions to the energy DE_k^cR, DE_k^RR and the first contribution of Eq. (43) can be grouped in the total correction DE_k^coul, which corresponds to the second order interaction between charges in the R-subsystem due to polarization in the M-subsystem. If the ZDO approximation is used this contribution can be written as:

DE_k^coul =
å
m Î M

å
ii^¢ Î { r,l}

å
AA^¢ Î R
g_{AI_m}g_{A^¢I_m^¢}Q_A^RQ_A^¢^R
å
st
P_{i_msi_ms^¢i_mt^¢i_mt}^M(0),
(0)
where I_m denotes the atom of the m-th bond corresponding to the i-th HO of the bond (R_m or L_m).

The next two contributions of Eq. (40), including the operator V^rr(q,E), correspond to the third order in the perturbation expansion. Therefore, they may be only roughly estimated. The main assumptions are the ZDO approximation and that the energies of an electron transfers between the subsystems are replaced by average ones depending only on the number of bond taking part in the electron transfer. The cross term between the V^c(q) and V^rr(q,E) is:

DE_k^crr = áY_k^appr | V^c(q)P QR(q,0)QPV^rr(q,E) | Y_k^apprñ »

»
å
pqq^¢ Î R

å
i_mj_mj_m^¢ Î M
v_{qj_m}v_{q^¢j_m^¢}
å
ss^¢t
(p_s p_s ||i_ms^¢i_ms^¢)×

×{(I_m-A_k-g_mk)^-1[P_k^spp(1-P_k^{tq^¢q})P_{i_ms^¢i_ms^¢j_mtj_mt^¢}^M(0)-

- 2
p
ó
õ P_{p_s p_s q_t ^¢q_t}^Rk(iu)P_{i_ms^¢i_ms^¢j_mtj_mt^¢}^M(iu)du]-

-(I_k-A_m-g_km)^-1[P_k^sppP_k^{tqq^¢}P_{i_ms^¢i_ms^¢ [(j)\tilde]_mt[(j)\tilde]_mt^¢}^M(0)+

+ 2
p
ó
õ P_{p_s p_s q_t q_t ^¢}^Rk(iu)P_{i_ms^¢i_ms^¢[(j)\tilde]_mt [(j)\tilde]_mt^¢}^M(iu)du]},

(0)
where [(j)\tilde] = r for j = l and [(j)\tilde] = l for j = r. The integral expressions can be replaced by sums using Eq. (44).

Analogously, the cross term between V^R(q) and V^rr(q,E) equals to:

DE_k^Rrr = áY_k^appr | V^R(q)PQR(q,0)QPV^rr(q,E) | Y_k^apprñ »

» -
å
st

å
pp^¢ Î R

å
A Î R
Z_A^R
å
j_mj_m^¢ Î M
v_{pj_m}v_{p^¢j_m^¢}×

×{(I_m-A_k-g_mk)^-1(1-P_k^{tp^¢p})×

×( g_{AR_m}P_{r_msr_msj_mtj_mt^¢}^M(0)+g_{AL_m}P_{l_msl_msj_mtj_mt^¢}^M(0)) -

-(-1)^{d_jj^¢}(I_k-A_m-g_km)^-1P_k^{tpp^¢}×

×( g_{AR_m}P_{r_msr_ms[(j)\tilde]_mt[(j)\tilde]_mt^¢}^M(0)+g_{AL_m}P_{l_msl_ms[(j)\tilde]_mt[(j)\tilde]_mt^¢}^M(0)) }.

(0)
The last contribution to the energy of Eq. (40) is of the fourth order with respect to the perturbation operator (counting four small operators of electron transfer) and may be totally neglected.

Discussion

Large molecular systems are very important from the practical point of view. Polymers and many biological molecules contain hundreds and thousands atoms. Modern ab initio methods do not allow to calculate the electronic structure of such systems because the computational time for these methods grows as N⁴¸N⁷, where N is number of one-electron basis functions. Even in the case of the standard semiempirical methods based on the Hartree-Fock approximation this growth is proportional to N³. It explains the especial role of the hybrid QM/MM schemes in the computational chemistry. There is a large diversity of the hybrid QM/MM schemes proposed in the literature due to the form of junction between quantum and classical subsystems. It leads to divergence of results obtained with different junctions. The results of the calculations seems to be ambiguous due to uncontrolled influence of the junction on the calculated PES of the molecular system.

We try to justify the form of the junction between subsystems. The complete neglect of the junction, i.e. the use of the bare (''ab initio'') Hamiltonian for the R-subsystem and standard molecular mechanical procedure without renormalization leads to the energy E_k0^R+E₀^M instead of the correct value E_k. Taking into account the interaction between the system (junction) results in a correction to the PES, which can be represented as a sum of contributions:

DE_k » DE_k^el-st+DE_k^rr+DE_k^disp+DE_k^coul+DE_k^crr+DE_k^Rrr,
(0)
where the main contribution DE_k^el-st is of the first order with respect to the perturbation operator and originates from the electrostatic interaction between the subsystems; the next four contributions are of the second order and follow from the coupling of one-electron transfers between subsystems, of Coulomb electron-electron interactions between them, of interactions between electrons of the M-subsystem and cores of the R-subsystem; and of interaction of the Coulomb electron-electron interaction with the interaction of the electrons of the M-subsystem and cores of the R-subsystem; the second term can be regarded as dispersion interaction between subsystems. Moreover, the contributions DE_k^disp and DE_k^coul can be considered as the second order interaction between the electronic polarization in the M-subsystem and the polarized R-subsystem. The last two contributions to Eq. (52) correspond to the third order in the interaction between the subsystems originating from the coupling two projection procedures. Physically it corresponds to interaction between two one-electron transfers and Coulomb interaction between electrons of the M-subsystem with electrons and cores of the R-subsystem.

One of us (A.M.T.) acknowledges financial support from the Haldor Topsø e A/S.

References

[]: CECAM-NSF Meeting on QC/MM methods, Int J Quantum Chem 1996, 60, No 6, Special Issue.
[]: Aqvist J.; Warshel A. Chem Rev 1993, 93, 2523.
[]: Bakowics D.; Thiel W.; J Comput Chem 1996, 17, 87.
[]: Carmer C.S.; Weiner B.; Frenklach M. J Chem Phys 1993, 99, 1356.
[]: Field M.J.; Bach P.A.; Karplus M. J Comput Chem 1990, 11, 300.
[]: Field M.J. In Computer Simulation of Biomolecular Systems: Theoretical and Experimental Applications, van Gunsteren W.F.; Weiner P.K.; Wilkinson A.J., eds. ESCOM: Leiden, 1993.
[]: Gao J.L. In Reviews in Computational Chemistry Lipkowitz K.B.; Boyd D.B., eds. VCH Publishers: New York, 1996.
[]: Gao J.L. J Phys Chem 1992, 96, 537.
[]: Gao J.L.; Xia X.F. Science 1992, 258, 631.
[]: Luzhkov V.; Warshel A. J Comput Chem 1992, 13, 199.
[]: Maseras F.; Morokuma K. J Comput Chem 1995, 16, 1170.
[]: Stanton R.V.; Hartsough D.S.; Merz K.M., Jr. J Comput Chem 1995, 16, 113.
[]: Thompson M.A. J Phys Chem 1995, 99, 4794.
[]: Tuñon I.; Martins-Costa T.C.; Millot C.; Ruiz-Lopez M.F.; Rivail J.-L. J Comput Chem 1996, 17, 19.
[]: Truong T.N.; Stefanovich E.V. Chem Phys Lett 1996, 256, 348.
[]: Morokuma K.; Froese R.D.J.; Humbel S.; Svensson M.; Matsubara T.; Sieber S. J Phys Chem 1996, 100, 2573.
[]: Bernardi F.; Olivucci M.; Robb M.A. J Amer Chem Soc 1992, 114, 1606.
[]: Tchougréeff A.L. Khim Fiz 1997, 16, 62 (in Russian); Chem Phys Reps 1997, 16, 1035 (in English).
[]: Tchougréeff A.L. Phys Chem Chem Phys 1999, 1, 1051.
[]: Tokmachev A.M.; Tchougréeff A.L.; Misurkin I.A. J Mol Struct (Theochem) 2000, 506, 17.
[]: Löwdin P.-O. J Math Phys 1962, 3, 969.
[]: Surján P.R. The Concept of the Chemical Bond, In Theoretical Models of Chemical Bonding; Part 2 Z.B. Maksi\'c, ed. Springer: Heidelberg, 1989.
[]: Wu W.; McWeeny R. J Chem Phys 1994, 101, 4826.
[]: Tokmachev A.M., Tchougréeff A.L. Zhurn Fiz Khim 1999, 73, 347 (in Russian); Russ J Phys Chem 1999, 73, 320 (in English).
[]: Tokmachev A.M.; Tchougréeff A.L. J Comput Chem (Accepted).
[]: Tokmachev A.M.; Tchougréeff A.L. (In preparation).
[]: Weinbaum S. J Chem Phys 1933, 3, 547.
[]: Tokmachev A.M.; Tchougréeff A.L.; Misurkin I.A. Int J Quantum Chem (Submitted).
[]: McWeeny R. Methods of Molecular Quantum Mechanics; 2nd Edition; Academic Press: London, 1992.
[]: Le Fevre R.J.W. Molecular Refractivity and Polarizability, In Advances in Physical Organic Chemistry, Vol. 3, Gold V. ed. Academic Press: London and New York, 1965.

File translated from T_EX by T_TH, version 2.67.
On 4 Oct 2000, 15:35.

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

Potential Energy Surfaces in Hybrid QM/MM Methods

Abstract

Introduction.

PES of composite quantum-classical system.

Discussion

References

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