The mechanistic approach to the energy of molecules based on the combination of classical force fields is a useful tool for analysis of the potential energy surfaces for large molecular systems. It is now called molecular mechanics (MM) and exists in many incarnations: CHARMM [], OPLS [], MM3 [], AMBER [] to mention a few. The particular schemes differ in the number and structure of the force fields as well as in the sets of parameters used. It is important to stress the purely empirical character of standard MM constructs. It allows to make these schemes computationally feasible but it also precludes understanding and detailed analysis of physical interactions in the molecule. In practice, the question whether to include or not of some ''weak'' force fields like ''off-diagonal'' interbond interactions is solved on the ''school-wise'' basis [].
The low computational cost of the MM treatment of molecules (especially with respect to high quality quantum mechanical (QM) computations) has led to an idea [] of constructing hybrid QM/MM procedures where small part of the molecule is calculated by using precise QM method while the rest is described by relatively cheap MM procedure. During past decades the progress in this field was enormous which is clearly shown in reviews Gao,Sherwood,PTCP. At the same time the question about the interface between subsystems and the account of their interactions is not completely solved []. In practice different recipes based on the intuitive models are used and their success in description of that or another system is considered as an ultimate criterion of their value. In a series of articles [,,] we proposed another way of thinking: the junction between subsystems should be obtained as a result of sequential derivation using explicit separation of electronic variables. This approach based on the effective Hamiltonian technique is quite general.
The essential features of the approach [,,] are (i)
unique assignment of one-electron states (which are chosen to be hybrid
orbitals) to the subsystems and (ii) explicit using of underlying wave
function for the inert (MM) subsystem. The key moment of the required
underlying QM procedure is its local character recovering the common
concepts of chemical bonds and lone pairs. The antisymmetrized product of
strictly localized geminals (APSLG) wave function is used here as an
underlying function for the MM. Its semiempirical implementation, we rely
upon, has been described in Refs. [,]. Here we give
only a brief account of this method. The APSLG wave function is constructed
as the antisymmetrized product of two-electron wave functions (geminals):
| (0) |
| (0) |
| (0) |
| (0) |
| (0) |
| (0) |
The use of the APSLG wave function as an underlying one for the MM subsystem is based not only on the internal similarity between the APSLG and the MM energy expressions but also on the possibility of classical representation of the potential energy surface (PES) of an organic molecule on the basis of the APSLG approach. The key feature of the semiempirical implementation of the APSLG wave function, we rely upon, differing it from others proposed in the field [] and previously used to analyse additive systematics (say in Ref. []) is the variational treatment of the strictly localized HOs which are determined on equal footing with other electronic structure parameters (ESPs), namely, the parameters characterizing the bonds (amplitudes of the geminal components), i.e. by minimizing the molecular total energy. Also, the important step here is the possibility to use the guess values of the ESPs and to construct the non-iterative (MM-like) procedure based on the APSLG wave function []. The accepted way to treat the ESPs serves as a basis for classification of the approximate schemes []: FA or TA for the fixed or tuned geminal amplitudes respectively and FO or TO for the fixed or tuned hybrid orbitals. Essential progress was achieved in Ref. [] where the linear response approximation was used for analysis of interrelations between molecular geometry and hybridization for the sp3 carbon atom. This construct may be qualified as a deductive molecular mechanics since each of its components has a transparent counterpart in the underlying QM description and the approximations and simplifications used can be uniquely characterized and formulated. The analysis has shown that the APSLG approach can be a basis for introducing the transferable energy contributions thus settling the QM foundation for the ''force fields'' known in the MM realm. Moreover, it allows to find corrections to the transferable QM prototypes of the MM force fields and in some cases to substantiate the presence of certain ''weak'' force fields. It is important to stress that the values of the elasticity constants obtained in the framework of the deductive molecular mechanics and those accepted in the standard force fields are very close. At the same time the purpose of the deductive molecular mechanics is not merely a substantiation of the standard MM force fields but also the detailed study of hybridization patterns in different situations as emphasized in Ref. [].
In the present paper we extend the treatment proposed in Ref. [] to organic molecules containing sp3 nitrogen. The next Section results in a deductive MM description of the nitrogen-containing compounds with particular stress on the source of the pyramidalization potential. Then, the implications of our treatment for development of consistent hybrid QM/MM (quantum mechanics/molecular mechanics) methods for analysis of PESs of chemical transformations of large molecular systems are considered.
For a century two fundamental facts shaped the area of stereochemistry: the terahedral carbon introduced by van t'Hoff and Le Bel and the pyramidal nitrogen. Despite that long development history no common viewpoint upon source of the angular dependence of the energy has not been developed within MM and stereochemistry itself. On the other hand, even very simple quantum chemical methods reproduce the observed features. Also, in quite general terms it is clear that the form of the coordination polyhedron is controlled by relation between the bonding (two-center) energy which favors population of excited states of an atom under consideration and the excitation energies themselves which tend to keep an atom in its ground state. Recently these notions again have been brought to discussion in Ref. []. The Gillespie explanatory scheme [] must be mentioned in this context. According to it the angular dependence of energy appears due to Coulomb repulsion between electron pairs (valence shell electron pair repulsion - VSEPR). In the literature there exists quite a number of attempts to reconcile this qualitative and intellectually very attractive picture with the results of the quantum chemical calculations, which are reviewed in Ref. []. These attempts, however, turned out to be discouraging. It has been found that the intraatomic energy terms responsible for the molecular shape cannot be identified with the interpair Coulomb interactions. At the same time the stereochemistry (form of the coordination polyhedron) and related to it the hybridization pattern can be fairly reproduced with no relation to any Coulomb repulsion at all, as in, say, the extended Hückel method []. This finding applies both to the carbon stereochemistry described previously and to the nitrogen one. An attempt to derive sequentially the stereochemistry rules with the forms of hybrids obtained semiquantitatively and a posteriori was undertaken in Ref. []. This simplified analysis based on the s-/p- ratio in HOs has led to qualitative description of stereochemistry with driving force differing from that claimed within the VSEPR approach. It was qualitatively shown in Ref. [] that the known stereochemistry rules in fact rely on the relative strength of electron attraction to the core on the s- and p-orbitals, respectively. In what follows we analyze different contributions to the energy of the APSLG description of the model ammonia molecule and demonstrate that the same energy terms are responsible for determining the observed stereochemistry and specifically are the source of the pyramidalization potential of the nitrogen atom.
As it has been shown in Ref. [] the form of the HOs within the FA class of the approximations leading to the MM-like description of the four-coordinated carbon atom is ultimately defined by the two-center resonance interactions. The SO(4) group structure of the hybrid manifold significantly restricts the ability of the one-electron states residing on the atom to adjust themselves to the variations of geometry arrangement of the surrounding atoms (groups). For some variations (hybridization incompatible) it is not possible to construct a system of HOs which would even approximately fit the direction of the internuclear axes. This rigidity, however, has nothing to do with the Coulomb repulsion of electrons populating the HOs. Moreover, in the case of the FA approximation the sum of the intraatomic Coulomb interaction terms does not depend on the hybridization at all. No significant bond-bond interaction occurs either to maintain the observed forms of the hybridization patterns (tetrahedra) and thus of the molecules themselves, rather than the topology of the hybridization manifold which assures the latter. The situation with other organogenous atoms seems to differ significantly from this picture. In the case of nitrogen and oxygen atoms even in the FA class of approximations the one-center energy is strongly hybridization dependent. But in contradiction with the VSEPR conjecture this happens due to the one-electron terms describing the core attraction of electrons in the lone pair and sensitive to the relative weights of the s- and p-AOs in the corresponding HO. The source of this sensitivity is of course strong difference between the core attraction in the s- and p-subshells with large preference towards purely s lone pair (in a line with qualitative reasoning of Ref. Mayer). In the free atom (purely one-center energy) this immediately resulted in no hybridization for nitrogen and oxygen and in a 90° valence angles predicted by older theories [] for water and ammonia with subsequent need to explain the observed form of these molecules with the valence angles only slightly smaller than the tetrahedral one and both exceeding 100°. Curiously enough that the authors of the VSEPR model seem to overlook this result well known for decades and did not consider it as a starting point and incidentally the limiting case of the electron pair repulsion and started their theory from a scratch.
If we reside for a while in the FA domain we have to admit that the only
source of the observed stereochemistry can be found in the interplay between
one-center hybridization dependent terms and the resonance energy, which was
clear yet to Coulson []. For nitrogen (in ammonia) one can
easily write the part of the electronic energy which depends on the
hybridization of the nitrogen atom:
| (0) |
It is instructive to parameterize the hybridization manifold. This problem
was adderssed in Refs. [,] where the parameterization of the SO(4) group by Jacobi angles was used. This six-parametric group can be
parameterized by the angles wsxA, wsyA, wszA, wyzA, wxzA, and wxyA
corresponding to rotations in two-dimensional subspaces of the whole
4-dimensional space spanned by one s- and three p-functions. The first
triple of (pseudorotation) angles ([(w)\vec] bA) defines a
hybridization tetrahedron - a shape uniquely defining relative weights of
the s- and p-functions in the hybrids and by this also the interhybrid
angles, while the second triple of (quasirotation) angles is responsible for
rotation of the hybridization tetrahedron as a whole (The prefix quasi
refers to the fact that no physical body rotates under its action - only
the hybridization tetrahedron defined right above). It is known that the SO(4) group is a direct product of two SO(3) subgroups. Using an
intermediate parameterization of the whole SO(4) manifold by a pair of
quaternions in Ref. [] allowed to determine the first order
correction to the HOs due to small quasi- and pseudorotations d[(w)\vec] l and d[(w)\vec] b:
| (0) |
For the symmetry reasons the resonance energy of the C3v ammonia
molecule is a function of only one pseudorotation angle wsz and
of the pyramidalization angle d which equals to zero for the planar
ammonia molecule. The resonance energy as a function of hybridization and
geometry acquires the form (assuming equivalence of all bonds):
| (0) |
The hybridization-dependent part of the one-center energy of the nitrogen
atom is:
| (0) |
Minimization of the sum of resonance Eq. (9) and one-center
Eq. (10) energies on wsz allows to determine the
optimal value of this pseudorotation angle as a function of the angle d. This minimization is not possible analytically. So, we try to find
realistic estimates for the dependence wsz(d). First of
all, we can determine the value of the pyramidalization angle d
optimal for a given value of pseudorotation angle wsz. It can be
easily done by taking the derivative with respect to d of the
resonance energy only since the lone pair contribution does not depend on d explicitly. It should be noted that the relatively small specific
correction to the core-core interaction adopted in semiempirical schemes is
pyramidalization angle-dependent. Also, small polarization of the N-H bonds
within the TA approximation framework modifies the energy expression. We
neglect the effect of these contributions. It allows to obtain a simple
relation:
| (0) |
| (0) |
The second derivatives matrix Ñ[(w)\vec] 2E which is necessary
for construction of linear response relations between the geometry
variations and the can be readily calculated for the sp3 nitrogen atom
in the vicinity of the equilibrium values of the quasi- and pseudorotation
angles. It can be shown that the second order correction to the resonance
energy is an average of the 6×6 matrix over vector parts of the
intermadiate quaternions p and q (see for details Ref.
TchSO4). It has the form:
| (0) |
| (0) |
| (0) |
| (0) |
It is interesting to find how the equilibrium value of pseudorotation angle wsz is reproduced in the linear response approximation:
| (0) |
| (0) |
One of the purposes of the present work is to elucidate the precise form of the junction between the QM and MM regions in the hybrid QM/MM methods. The latter is more or less covered by the polarization effects when it goes about the junctions across space, i.e. when no chemical (covalent) bond binds the regions to be treated by the MM and QM, respectively. An opposite situation arises when there are chemical bonds between the QM and MM regions like in the case of enzymatic reactions or even reactions of large organic molecules where only a small fraction of bonds breaking or forming is considered with use of QM. In these cases it is now accepted that the interregion boundary must not cut bonds in halves, but is assumed either to pass ''through nuclei'' leaving some of the bonds completely in the MM region and keeping others completely in the QM one or is replaced by a kind of the ''grey zone'' for which neither fully MM nor fully QM treatment applies. Theoretical substantiation of namely a ''through nuclei'' choice of the interregion border has been given in Ref. []. However, in either of approaches there remains a problem of how to select the HOs residing on the frontier atoms which do not belong to the QM region completely and how the taking part of an atom in the QM affects its MM parameters. The list of recipes existing in the field is quite impressive and includes ''capping'' hydrogen atoms, ''link'' atoms, etc. In all these cases the question arises how to set the parameters of these atoms (particularly if advanced ab initio schemes are employed to treat the QM part). Also it is clear that the HOs included in the QM region but residing on the frontier atoms must adjust at least their orientation to the variation of the positions of the MM region atoms attached to the frontier atom. This obvious fact did not escape the attention of the workers in the field (see Ref. [] and references therein). However, the solutions of the problem proposed so far in the literature largely reduce to ad hoc recipes of adjusting the orientation of the QM residing HOs to the geometry variations in the MM part. These recipes are as well introduced without too much explanation or substantiation. An obvious part completely missed by this approach is the renormalization of the MM part of the system due to variations of the ESPs on the frontier atom caused by its participation in the QM subsystem.
The adjustment of the QM residing HOs must follow the variational principle for the energy when the geometry (there is nothing else) of the MM part changes. Even in such a setting the energy parameters at the QM residing HOs vary following the variations of the geometry changes in the MM part since the latter may result also in variation of the form of the HOs (i.e. of the relative weights of the s- and p-contributions to the HOs). Our previous approach [,] describes only the energy adjustment with pretty unclear prospects of including the variation of the form of the HOs in the QM region. Here we by contrast consider the modification of the form of the HOs and the renormalization of the QM and MM parameters corresponding to the frontier atom. More precisely, we consider a special case when the frontier atom serving as the QM/MM junction is the sp3 nitrogen atom supplying its lone pair to the QM subsystem. Such a setting seems to be quite natural since the basicity or the nucleophilicity functions of the nitrogen atom are both due to interaction of its lone pair with acceptor orbitals. This interaction is naturally to be treated by some kind of QM technique while leaving the rest of the nitrogen neigbours in the MM region.
The one- and two-electron density matrix elements for the QM residing HO are
evaluated by the corresponding (QM) procedure thus invoking the TA type of
description. Those for the MM region HOs of the atom under consideration are
fixed at their invariant values according to the FA setting. The
nonvanishing intersystem two-electron densities reduce to the products of
the corresponding one-electron density matrix elements. With these
assumptions we get the corrected hybridization-dependent one-center energy
for the frontier nitrogen atom:
| (0) |
| (0) |
The derivatives of the energy correction with respect to the angles [(w)\vec] b, [(w)\vec] l yield additional quasi- and pseudotorques
acting upon the hybridization tetrahedron:
| (0) |
In the QM part of the system the variation of the bond orders can take place
as well. In variance with the pure APSLG picture [,]
accepted in the present work as the underlying QM method for the MM part of
the system the atoms in the QM part of the combined system may have
off-diagonal elements of the one-electron density matrix between arbitrary
one-electron states ascribed to the QM subsystem. The latter are obviously
the bond orders for the QM part of the system. The corresponding
contribution to the energy reads:
| (0) |
| (0) |
| (0) |
| (0) |
The total pseudo- and quasitorques then become:
| (0) |
On the other hand any deformation in the MM system results in the variation
of the pseudo- and quasirotation angles. The shifts of the positions of the
MM neighbours of frontier atom result in quasi- and pseudotorques acting
upon the hybridization tetrahedron. In its turn this produces variations of
both one-center parameters correspoding to the QM residing HO and of the
resonance paprameters for the QM residing HO and all other orbitals in the
QM region. The variation of the one-center matrix elements of the
Hamiltonian corresponding to the QM HO is:
| (0) |
| (0) |
The renormalization of the QM resonanse intergrals to which the HO at hand
is involved is somewhat more complex. It nevertheless uses the same (DCF)
representation of the resonance integrals as previously:
| (0) |
In the present paper we extended the QM substatiation of MM of aliphatic hydrocabons proposed in Ref. [] to compounds containing aliphatic nitrogen. We analyzed the energy expression for ammines and revealed sources of pyramidality of the nitrogen atom. That allowed to construct approximate (interpolation) formulae for the pyramidalization potential for the sp3 nitrogen atom departing from the APSLG description of electronic structure of ammonia molecule and related (locally) C3v tertial ammines. From another point of view it leads to the special treatment of frontier atoms responsible for the junction between the QM and MM parts of the system under the hybrid QM/MM treatment. The effect of deformation in the MM-treated part upon the parameters of the effective QM Hamiltonian and the reverse effect of the QM driven variation of the elements of the one- and two-electron density matrices upon the MM treated part of the combined system have been considered and explicit formulae and numerical estimates of these effects have been presented.
Acknowledgments
This work has been performed with financial support of RAS through the grant # 6-120 dispatched by the Young Researchers' Commission of RAS. ALT gratefully acknowledges valuable discussions with Prof. A.A. Levin, Dr. I.V. Pletnev. The authors are grateful to Prof. I. Mayer for sending a reprint of his work Ref. [].