- Email: [email protected]
Bond orbital framework for rapid calculation of environmental effects on molecular potential surfaces
15 April 1983
CHEWCAL PHYSXS LETTERS
Volume 96, number 4
BOND ORBITAL FRAMEWORK FOR RAPID CALCULATION
OF ENVIRONMENTAL
EFFECTS
ON MOLECULAR POTENTIAL SURFACES G5bor N&Y-SZABd CHlh’olN
P?zannaceurical
and Peter R. SUFUh and Chemical
Works, H-1325
Budapest.
P-0. B. 1 ZO, Hungag?
Received 8 December 1982; in final form 8 February 1983
A procedure is outlined to calculate rapidly potential mrfaces of very large (bio)molewles Using strictly localized orbitals as a basis set, the molecule is divided into a central part, treated at the SCF level. and the environment where the strictly localized character of the orbitals is maintained. Conformation of the active serine sidechain in a-chymotrypsin is discussed.
1. Introduction Molecules containing more than fifty non-hydrogen atoms are too large to be treated by conventional quantum-chemical methods. Procedures to handle such systems by dividing them into a central part, subject
to a quantum-chemical calculation, and its environment mimicked by a set of point charges or multipole moments, are known in the literature [l-3] _ These methods suffer from two drawbacks_ First, the model of the environment may be crude and, second, isolating the central part is problematic_ In this paper we propose a new method which does not possess these drawbacks. Division into central and environmental parts becomes straightforward if a basis set of strictly localized bonding and antibonding molecular orbitals is used. At present, we formulate the theory within the zero-differential overlap approximation and we apply the CNDO/2 parameterization.
tion exclusively from SLBOs is a rough approximation but it works reasonably in describing, for example, electrostatic potential maps [8-lo]_ If bend polarities are determined fully variationally, net atomic charges can be calculated [ 1 l] _This approach corresponds to the zeroth-order PCILO method [7] _However, conformational properties are incorrectly predicted since tails, completing SLBOs to exact localized MOs obtain-
2. Method Let us construct in-phase and out-of-phase combinations of atomic hybrid orbitals (HYOs) which are orthogonalized to other hybrids centered on the same atom. The combinations are called strictly localized bond orbitals (SLBOs) and their use in quantumchemical methods has been proposed by several authors [4-l l] _Construction of the total molecular wavefunc0 009-2614/83/0000-0000/S
03.00 0 1983 North-Holland
able by a linear transformation
from the solutions of
the SCF equations, play a decisive role in this respect
[12-l 6]_ Accordingly, in order to describe conformational properties, delocalization effects are to be accounted for. Recent analyses [14,15] show that different kinds of delocalization effects have different significance in governing conformations_ In an orthogonal basis, such as the ZDO one, barriers to internal rotation arise primarily due to through-space (direct) delocalization interactions (TSI) [13-151. Through-bond iuteractions (TBI) can lead to significant delocalizations but the corresponding “tails” are rather insensitive to the relative bond orientations [IS] _The orientation-dependent tails are those due to TSI but they vanish rapidly with increasing distance_ The characteristic radius of action of TSI is ==3-500 pm [15,16]. From this it follows that partition into subsystem and environment is possible. To treat the subsystem as precisely as possible we apply the conventional SCF procedure for it. In order to make partitioning straightforward we use a basis set of bonding and antibonding SLBOs. Then the 499
subsystem is chosen as a closed-shell ensemble of SLBOs. directly affected by the chemical rearrange-
the environment block of the density matrix as given in eq. (2), eq. (4) involves computational work which is proportional onIy to the first power of the number of bonds in the environment, nE_ However, due to the optimization of the bond polarities by solving 2 X 2 factorized secular equations [11,16] the computationaI work in our program became proportional to u;_ Since the proportionality constant is small, rapid calculations for relatively large systems become possible_
ment .
We decompose iollows:
15 April 1983
CHEMICAL PHYSICSLETTERS
Volume 96. number 4
the fockian for the subsystem
as
3. Application to ar-chymotrypsin
where i.,.
k and I refer to SLBOs. while C and E denote the centrid part and environment, respectively. 11 is the core hannltonian. Partition of the molecule nit0 C zind E. combined with the orthogonality of SLBOs. allows us to write the following for the densit> matri\ P :
I’/$, = 16k,
it-k. I E E and X-.I occupied.
=0
if k. I E E and k or I virtual.
=0
ifkEEandIECor IEEandkEC.
(2)
whereas f’,-, (i. j E C) is treated as a full block in the densit> matrix of the whole system. Introducing the ZDO condrtion and eq. (2) into eq. (1) we obtain * i-$=-i+’
•r-x-.,~ c c [‘Al l(iiW)
\v nh an effective hamiltonian
-f
(ikliO1.
(3)
To analyze the practical applicability of our method we studied the rotation of the active serine side chain around the C?Cfl bond in cu-chymotrypsin [ 16]_ Our molecular model is depicted in fig. 1. thrre-dimensional coordinates were taken from the Protein Data Bank [ 1S ] as submitted by Tulinsky and co-workers [ 19]_ We partitioned the whole model into five parts. Region I consists of the rotating group and directly attached atoms. Region II is the o-type lone pair (a single HYO) at NE2 while region III includes the imidazole moiety but two relatively far-lying bonds. At last, region IV is the NH bond of the protein backbone while the remaining fifth part is the rest of the whole model. Fig. 2 depicts the rotational potential curves as obtained by our method cutting out three different subsystems from the model. Conventional full SCF results are also given for reference. If we treat region I
b
defined as
WC Applied the fockian of eq. (3) in the CNDO/2 p.rrameteri,zrtion [ 171 throughout this work. The total energy of the system WJS obtained as the exact especration value of the ltaniiltonian using the special density matrix of eq. (2). Due to the simple form of
.-
_p===c ; CH3
bonding-anrlbondmg
b&s integrals of Qpe trr*lkk). {ii’]kk*) do sunive the ZDO approGmation (the st.u refers IO dntibondinp orbit&).
x Note
500
rlmr
in
the
Fig. 1. Molecular model of the active site of a-chymotrypsin [ 18.191. Lone-pair SLBOsare indicated by double dots. For definition of rqions I. II. III and IV. see text.
.
CHEMICAL
Volume 96, number 4
an IBM 303 1 computer_ This is reduced to 2 min if regions I, II, III and IV are cut out of the model and to 1 mm if only regions I and II are treated. More can be spared if the whole molecule is larger. Our computer program is far from optimal and work to improve it is m progress.
ml eneqyTy+IkJlmoll
:
15 April 1983
PHYSICS LETTERS
;
Acknowledgement Discussions with Drs. J_ Angyin and I. Mayer (Budapest) are gratefully acknowledged.
-180
-120
-60
0
60
120
Fig. 2. Torsional potential curves for the CW6 rotation in the model of fig. 1 as obtained by various approximations. p= 0” corresponds to tbe K-ray geometry [ 18.19]_ The curdes are shifted to have absolute minima at the zero level. Absolute energy values (in kJ/mol) for 9 = -90” are given in parentheses. (a) 1 + II (-602.547); (b) I + II+ III (-603.691); (c) I + II + III + IV (-603.596): (d) SW reference (-606.675). only as the central part, irrational results are obtained, therefore we do not even give the corresponding torsional potential curve in fig_ 2. This failure is due to the formation of the 0rNE2 hydrogen bond involving considerable charge transfer between His-57 and Ser195. However, if region II, the single lone pair at NE2 is included into the subsystem and a qualitatively reasonable torsional potential curve, indicated by b in fig. 2, is obtained_ Location of minima and maxima become practically identical with those obtained from the reference SCF calculation if region III, the imidazole ring is also treated in the central part. At last, the agreement with the SCF curve is almost quantitative, the maximal deviation being 8 kJ/mol, if region IV is also considered. Summing up, we may state that the simplified method, presented in this paper, is applicable to the conformational analysis of large molecules. Some chemical intuition is needed to select the central region properly. This may be confumed in problematic cases by comparing approximate results to exact SCF calculations on suitable models as done above. The computer time is considerably reduced: a standard SCF CNDO/2 calculation [17] for the whole model in fig. 1 for one point on the potential curve costs 15 rnin CPU using
References [ 1J D.M. (1976)
Hayes and P.A. KoIlman. 1. Am. Chem. Sot.
98
7811.
111 0. Tapia, F_ Sussman and E. Poulain. J. Theoret. Biol. 7 1(1978) 49. [3] P-T_ Van Duijnen, B-T_ Thole and W.G.J. HOI, Biophys Chem. 9 (1979) 273. [4 ] G-G. HaII. Proc. Roy_ Sot. A205 (1951) 54 1. [5] C. Sandorfy. Canad. J. Chum. 33 (1955) 1337. [6] G. Del Re. J. Chem. Sot. (1958) 4031. [ 71 J-P. MaIrieu, in: hiodern theoretical chemistry, Vol. 7. ed. G.A. Se@ (Plenum. New York. 1977) p_ 69. [ 8 ] R. Bonaccorsi. E. Scrocco and J. Tomasi. J. Am. Chem. Sot. 98 (1976) 4049. [9] G. Ndray-Szab6. Intern. J. Quantum Chem. 16 (1979) 265. [ lo] G. Ndray-Srabb, A. Grofcsik, K. K&a, hi. Kubinyi and A. Martin, J. Comput. Chem. 2 (1981) 58. [ 11) P.R. Surjrin. hl. Rh&z and I. hlayer. J. Chem. Sot Faraday II 77 (1981) 1129. [ 121 W. England and h1.S. Gordon. J. Am. Chem. Sot. 93 (1971) 4649; 94 (1972)4818_ 1131 T-K. Brunck and F. Weinbold. J. Am. Chem. Sot. 101 (1979) 1700. [ 141 P-R_ Surjdn and I. Mayer, Theoret. Chim. Acta 59 (1981) 603. [ 1.51 P-R. SurjBn. I. Mayer and Y. Kertesz. J. Chem. Phys 77 (1982) 2454. [ 161 P-R Surjin. G. Ndray-Szab6 and 1. Mayer. Intern. J. Quantum
Chem. 12 (1962)
929.
[ 17) J-A. Pople and D-L. Beveridge, Approximate
molecular orbital theory (hlcGraw-Hi& New York, 1970). [ 181 F-C. Bernstein, T-F. Koetzie, G.J.B. Williams, E.F. Mayer Jr-, h1.D. Brice, J.R. Rogers, 0. Kemrard, T. Shimanouchi and hf. Tasumi. J. Mol. BioL 112 (1977) 535. [ 191 A. Tulinsky, N-V_ Mani. C.N. Morimoto and R.L. Vandlen, Acta Cryst. B29 (1973) 1309_
501
CHEWCAL PHYSXS LETTERS
Volume 96, number 4
BOND ORBITAL FRAMEWORK FOR RAPID CALCULATION
OF ENVIRONMENTAL
EFFECTS
ON MOLECULAR POTENTIAL SURFACES G5bor N&Y-SZABd CHlh’olN
P?zannaceurical
and Peter R. SUFUh and Chemical
Works, H-1325
Budapest.
P-0. B. 1 ZO, Hungag?
Received 8 December 1982; in final form 8 February 1983
A procedure is outlined to calculate rapidly potential mrfaces of very large (bio)molewles Using strictly localized orbitals as a basis set, the molecule is divided into a central part, treated at the SCF level. and the environment where the strictly localized character of the orbitals is maintained. Conformation of the active serine sidechain in a-chymotrypsin is discussed.
1. Introduction Molecules containing more than fifty non-hydrogen atoms are too large to be treated by conventional quantum-chemical methods. Procedures to handle such systems by dividing them into a central part, subject
to a quantum-chemical calculation, and its environment mimicked by a set of point charges or multipole moments, are known in the literature [l-3] _ These methods suffer from two drawbacks_ First, the model of the environment may be crude and, second, isolating the central part is problematic_ In this paper we propose a new method which does not possess these drawbacks. Division into central and environmental parts becomes straightforward if a basis set of strictly localized bonding and antibonding molecular orbitals is used. At present, we formulate the theory within the zero-differential overlap approximation and we apply the CNDO/2 parameterization.
tion exclusively from SLBOs is a rough approximation but it works reasonably in describing, for example, electrostatic potential maps [8-lo]_ If bend polarities are determined fully variationally, net atomic charges can be calculated [ 1 l] _This approach corresponds to the zeroth-order PCILO method [7] _However, conformational properties are incorrectly predicted since tails, completing SLBOs to exact localized MOs obtain-
2. Method Let us construct in-phase and out-of-phase combinations of atomic hybrid orbitals (HYOs) which are orthogonalized to other hybrids centered on the same atom. The combinations are called strictly localized bond orbitals (SLBOs) and their use in quantumchemical methods has been proposed by several authors [4-l l] _Construction of the total molecular wavefunc0 009-2614/83/0000-0000/S
03.00 0 1983 North-Holland
able by a linear transformation
from the solutions of
the SCF equations, play a decisive role in this respect
[12-l 6]_ Accordingly, in order to describe conformational properties, delocalization effects are to be accounted for. Recent analyses [14,15] show that different kinds of delocalization effects have different significance in governing conformations_ In an orthogonal basis, such as the ZDO one, barriers to internal rotation arise primarily due to through-space (direct) delocalization interactions (TSI) [13-151. Through-bond iuteractions (TBI) can lead to significant delocalizations but the corresponding “tails” are rather insensitive to the relative bond orientations [IS] _The orientation-dependent tails are those due to TSI but they vanish rapidly with increasing distance_ The characteristic radius of action of TSI is ==3-500 pm [15,16]. From this it follows that partition into subsystem and environment is possible. To treat the subsystem as precisely as possible we apply the conventional SCF procedure for it. In order to make partitioning straightforward we use a basis set of bonding and antibonding SLBOs. Then the 499
subsystem is chosen as a closed-shell ensemble of SLBOs. directly affected by the chemical rearrange-
the environment block of the density matrix as given in eq. (2), eq. (4) involves computational work which is proportional onIy to the first power of the number of bonds in the environment, nE_ However, due to the optimization of the bond polarities by solving 2 X 2 factorized secular equations [11,16] the computationaI work in our program became proportional to u;_ Since the proportionality constant is small, rapid calculations for relatively large systems become possible_
ment .
We decompose iollows:
15 April 1983
CHEMICAL PHYSICSLETTERS
Volume 96. number 4
the fockian for the subsystem
as
3. Application to ar-chymotrypsin
where i.,.
k and I refer to SLBOs. while C and E denote the centrid part and environment, respectively. 11 is the core hannltonian. Partition of the molecule nit0 C zind E. combined with the orthogonality of SLBOs. allows us to write the following for the densit> matri\ P :
I’/$, = 16k,
it-k. I E E and X-.I occupied.
=0
if k. I E E and k or I virtual.
=0
ifkEEandIECor IEEandkEC.
(2)
whereas f’,-, (i. j E C) is treated as a full block in the densit> matrix of the whole system. Introducing the ZDO condrtion and eq. (2) into eq. (1) we obtain * i-$=-i+’
•r-x-.,~ c c [‘Al l(iiW)
\v nh an effective hamiltonian
-f
(ikliO1.
(3)
To analyze the practical applicability of our method we studied the rotation of the active serine side chain around the C?Cfl bond in cu-chymotrypsin [ 16]_ Our molecular model is depicted in fig. 1. thrre-dimensional coordinates were taken from the Protein Data Bank [ 1S ] as submitted by Tulinsky and co-workers [ 19]_ We partitioned the whole model into five parts. Region I consists of the rotating group and directly attached atoms. Region II is the o-type lone pair (a single HYO) at NE2 while region III includes the imidazole moiety but two relatively far-lying bonds. At last, region IV is the NH bond of the protein backbone while the remaining fifth part is the rest of the whole model. Fig. 2 depicts the rotational potential curves as obtained by our method cutting out three different subsystems from the model. Conventional full SCF results are also given for reference. If we treat region I
b
defined as
WC Applied the fockian of eq. (3) in the CNDO/2 p.rrameteri,zrtion [ 171 throughout this work. The total energy of the system WJS obtained as the exact especration value of the ltaniiltonian using the special density matrix of eq. (2). Due to the simple form of
.-
_p===c ; CH3
bonding-anrlbondmg
b&s integrals of Qpe trr*lkk). {ii’]kk*) do sunive the ZDO approGmation (the st.u refers IO dntibondinp orbit&).
x Note
500
rlmr
in
the
Fig. 1. Molecular model of the active site of a-chymotrypsin [ 18.191. Lone-pair SLBOsare indicated by double dots. For definition of rqions I. II. III and IV. see text.
.
CHEMICAL
Volume 96, number 4
an IBM 303 1 computer_ This is reduced to 2 min if regions I, II, III and IV are cut out of the model and to 1 mm if only regions I and II are treated. More can be spared if the whole molecule is larger. Our computer program is far from optimal and work to improve it is m progress.
ml eneqyTy+IkJlmoll
:
15 April 1983
PHYSICS LETTERS
;
Acknowledgement Discussions with Drs. J_ Angyin and I. Mayer (Budapest) are gratefully acknowledged.
-180
-120
-60
0
60
120
Fig. 2. Torsional potential curves for the CW6 rotation in the model of fig. 1 as obtained by various approximations. p= 0” corresponds to tbe K-ray geometry [ 18.19]_ The curdes are shifted to have absolute minima at the zero level. Absolute energy values (in kJ/mol) for 9 = -90” are given in parentheses. (a) 1 + II (-602.547); (b) I + II+ III (-603.691); (c) I + II + III + IV (-603.596): (d) SW reference (-606.675). only as the central part, irrational results are obtained, therefore we do not even give the corresponding torsional potential curve in fig_ 2. This failure is due to the formation of the 0rNE2 hydrogen bond involving considerable charge transfer between His-57 and Ser195. However, if region II, the single lone pair at NE2 is included into the subsystem and a qualitatively reasonable torsional potential curve, indicated by b in fig. 2, is obtained_ Location of minima and maxima become practically identical with those obtained from the reference SCF calculation if region III, the imidazole ring is also treated in the central part. At last, the agreement with the SCF curve is almost quantitative, the maximal deviation being 8 kJ/mol, if region IV is also considered. Summing up, we may state that the simplified method, presented in this paper, is applicable to the conformational analysis of large molecules. Some chemical intuition is needed to select the central region properly. This may be confumed in problematic cases by comparing approximate results to exact SCF calculations on suitable models as done above. The computer time is considerably reduced: a standard SCF CNDO/2 calculation [17] for the whole model in fig. 1 for one point on the potential curve costs 15 rnin CPU using
References [ 1J D.M. (1976)
Hayes and P.A. KoIlman. 1. Am. Chem. Sot.
98
7811.
111 0. Tapia, F_ Sussman and E. Poulain. J. Theoret. Biol. 7 1(1978) 49. [3] P-T_ Van Duijnen, B-T_ Thole and W.G.J. HOI, Biophys Chem. 9 (1979) 273. [4 ] G-G. HaII. Proc. Roy_ Sot. A205 (1951) 54 1. [5] C. Sandorfy. Canad. J. Chum. 33 (1955) 1337. [6] G. Del Re. J. Chem. Sot. (1958) 4031. [ 71 J-P. MaIrieu, in: hiodern theoretical chemistry, Vol. 7. ed. G.A. Se@ (Plenum. New York. 1977) p_ 69. [ 8 ] R. Bonaccorsi. E. Scrocco and J. Tomasi. J. Am. Chem. Sot. 98 (1976) 4049. [9] G. Ndray-Szab6. Intern. J. Quantum Chem. 16 (1979) 265. [ lo] G. Ndray-Srabb, A. Grofcsik, K. K&a, hi. Kubinyi and A. Martin, J. Comput. Chem. 2 (1981) 58. [ 11) P.R. Surjrin. hl. Rh&z and I. hlayer. J. Chem. Sot Faraday II 77 (1981) 1129. [ 121 W. England and h1.S. Gordon. J. Am. Chem. Sot. 93 (1971) 4649; 94 (1972)4818_ 1131 T-K. Brunck and F. Weinbold. J. Am. Chem. Sot. 101 (1979) 1700. [ 141 P-R_ Surjdn and I. Mayer, Theoret. Chim. Acta 59 (1981) 603. [ 1.51 P-R. SurjBn. I. Mayer and Y. Kertesz. J. Chem. Phys 77 (1982) 2454. [ 161 P-R Surjin. G. Ndray-Szab6 and 1. Mayer. Intern. J. Quantum
Chem. 12 (1962)
929.
[ 17) J-A. Pople and D-L. Beveridge, Approximate
molecular orbital theory (hlcGraw-Hi& New York, 1970). [ 181 F-C. Bernstein, T-F. Koetzie, G.J.B. Williams, E.F. Mayer Jr-, h1.D. Brice, J.R. Rogers, 0. Kemrard, T. Shimanouchi and hf. Tasumi. J. Mol. BioL 112 (1977) 535. [ 191 A. Tulinsky, N-V_ Mani. C.N. Morimoto and R.L. Vandlen, Acta Cryst. B29 (1973) 1309_
501