ARTICLE pubs.acs.org/JPCA
Topography of Scalar Fields: Molecular Clusters and π-Conjugated Systems Sachin D. Yeole†,‡ and Shridhar R. Gadre*,† † ‡
Department of Chemistry, Indian Institute of Technology Kanpur, Kanpur 208016, India Department of Chemistry, University of Pune, Pune 411007, India ABSTRACT: The pioneering works due to Bader and co-workers have generated widespread interest in the study of the topography of molecular scalar fields, the first step of which is the identification and characterization of the corresponding critical points (CPs). The topography of a molecular system becomes successively richer in going from the bare nuclear potential (BNP) to the molecular electrostatic potential (MESP) through the molecular electron density (MED). The present work clearly demonstrates, through the study of some π-conjugated test molecules as well as molecular clusters, that the CPs could be economically located by following this path within ab initio level theory. Further, the topography mapping of large molecules, especially at a higher level of theory, is known to be a demanding task. However, it is rendered possible by following the above sequential mapping assisted by a divide-and-conquer-type method termed as the molecular tailoring approach (MTA). This is demonstrated with the topography mapping of β-carotene and benzene nonamer at MP2 and a (H2O)32 cluster at the HF level of theory, which are rather challenging problems with contemporary off-the-shelf computer hardware.
1. INTRODUCTION Study of molecular scalar fields is an active research area founded and extensively developed by Bader and co-workers.1 It is aimed at providing an understanding of molecular structure and reactivity. The scalar fields explored so far in the chemical literature include the bare nuclear potential (BNP), molecular electron density (MED), and molecular electrostatic potential (MESP), to name a few. Furthermore, some of these scalar fields are experimentally amenable and thereby provide a clear-cut bridge between theory and experiment. The real-life systems, such as molecules of biological interest, being large in size, involve a huge amount of numerical data to be examined for the study of a scalar field. This is where topography mapping is useful because it provides condensed information extracted from the vast numerical data of the scalar field. Topography mapping involves identification and characterization of the scalar field through isosurfaces, contours, gradient paths, and, most importantly, the critical points (CPs).13 The CP of a three-dimensional (3D) scalar field described by a function f(x1,x2,x3) is defined as the point P at which all three first-order partial derivatives vanish, that is, for i = 1, 2, 3 ∂f ðx1 , x2 , x3 Þ=∂xi ¼ 0 ð1Þ p
A CP is further characterized by evaluation of eigenvalues of the Hessian matrix Hij ¼ ð∂2 f =∂xi ∂xj Þ ð2Þ p
r 2011 American Chemical Society
at the CP, namely, P. A nondegenerate CP (where none of the eigenvalues are 0) is characterized by the rank (R), the number of nonzero eigenvalues, and the signature (S), the algebraic sum of the signs of eigenvalues of the Hessian matrix.13 Thus, for a 3D scalar field, a CP for which all of the eigenvalues of the Hessian matrix are positive is denoted as (3,+3), a minimum. Likewise, there can be three more types of nondegenerate CPs, namely, (3,3), a maximum, and two saddles, (3,+1) and (3,1).The CPs are of vital importance for the atoms-in-molecules (AIM) theory of Bader.1 BNP, Vbn(r), a 3D scalar function, is the nuclear contribution to MESP, also identified as the external potential for a molecular system within the HohenbergKohn theorem,4 given as Vbn ðrÞ ¼
∑α rαα Z
ð3Þ
MED is another scalar function of chemical interest that is of great conceptual value, forming the basis of the density functional theory (DFT), and can be determined from the corresponding wave function as Z jψðx, x 2 , ... , xN Þj2 d3 r2 , ... d3 rN ð4Þ FðrÞ ¼ N
∑σ
Here, x denotes jointly the spatial (r) and spin (σ) coordinates of an electron and N is the number of electrons. Parr and Berk5 Special Issue: Richard F. W. Bader Festschrift Received: April 27, 2011 Revised: August 21, 2011 Published: September 21, 2011 12769
dx.doi.org/10.1021/jp2038976 | J. Phys. Chem. A 2011, 115, 12769–12779
The Journal of Physical Chemistry A
ARTICLE
Figure 1. Molecular cluster test systems, namely, (CO2)5, (C6H6)4, (CH3OH)3, (HCl)4, and (H2O)9 clusters, along with their BNP CPs, the black dots denoting the (3,+1) CPs and the gray dots the (3,1) CPs. See text for details.
termed the BNP as a “harbinger” of the MED. Earlier, Parr et al.6 noticed that MED turns out to be a local function of BNP, within a simple local density functional theory. The scalar field of MED has been popularized over the past three decades by Bader and co-workers1,7 through their fundamental and exhaustive investigations. In their pioneering work, they have successfully exploited the topography of MED and the details of the Laplacian of MED, namely, r2F(r), for a crisp and compact description of the molecular structure. Their studies have revealed that a chemical bond is represented by a (3,1) CP, a ring in molecular framework shows a (3,+1) signature, and a cage gives rise to a (3,+3) CP in MED. Parallel to this development, Gadre et al. have investigated another 3D scalar field, namely, MESP,2,3 which is endowed with rich topographical features. MESP at a point is a sum of the corresponding BNP and negative potential engendered by the continuous electron charge density and is represented as Z N Zα Fðr0 Þ 3 0 dr ð5Þ V ðrÞ ¼ jr r0 j A ¼ 1 jr rα j
∑
MED as well as MESP for the ground states of neutral atoms is positive throughout.8 However, for molecular systems, MESP possesses some negative values in the reactive valence region, a unique feature that is not exhibited individually by the electronic or bare nuclear contributions. MESP exhibits negative-valued minima and saddles in the region of electron concentration, for example, for lone pairs and π-bonds. MESP cannot exhibit nonnuclear maxima,9 while MED for a few cases is known to show such a feature.10 Several years ago, Parr et al.6 brought out the similarity between contour diagrams of BNP and MED for the B2 molecule. This relationship was later investigated topologically by Tal et al.11 along the C2v dissociative path of the water molecule. Politzer and Zilles12 reported that similarity between BNP and MED does not hold for molecules such as cubane and oxirane.
Gadre and Bendale13 also examined these similarities, and they found that the contours of these two scalar fields appear to be strikingly similar over the ring regions for cyclopropane, oxirane, and so forth. Recently, Popelier et al.14 have reported an extensive study of the homeomorphism between MED and BNP for a large number of molecules. According to them, many molecules (but not all) show geometrical homeomorphism between BNP and MED topography. They also reported that cage minima, namely, (3,+3) CPs, cannot exist in a BNP distribution. The number of nondegenerate CPs in the topography of any 3D scalar field can be checked out through the PoincareHopf (PH) relationship nþ3 nþ1 þ n1 n3 ¼ χ
ð6Þ
Here, n+3 is the number of (3,+3) minima, n+1 is the number of (3,+1) saddles, and so on, and χ is the Euler characteristic of the scalar field. The value of χ is 1 for the scalar fields such as BNP and MED, which possess the same sign throughout space. For the scalar field of MESP, Libeouf et al. and Roy et al.15 proposed the PH relation as nþ3 nþ1 þ n1 n3 ¼ n nþ ¼ χ
ð7Þ
Here, n and n+ are the number of asymptotic maxima and minima, respectively. The PH relationship provides a necessary condition (although not sufficient) for the completeness of topography of a scalar field. For example, if a CP each of type (3,+3) and (3,+1) is missed, the PH relation will still hold true. The AIM analysis developed by Bader and co-workers starts with locating and characterizing the MED CPs. It also involves many further steps, such as tracing gradient paths, mapping out atomic basins, and so forth. Several packages are available for AIM analysis of molecules, such as the AIMAll suite of programs,16a MORPHY,16b XD2006,16c TOPXD,17a and BUBBLE,17b to name a few. All of these packages are being used to locate electron density critical points and properties, including ellipticities, energy densities, 12770
dx.doi.org/10.1021/jp2038976 |J. Phys. Chem. A 2011, 115, 12769–12779
The Journal of Physical Chemistry A
ARTICLE
Figure 2. The π-conjugated test systems, namely, heptacene, [18]annulene, C38H22, β-carotene, and C32H24S, along with their BNP CPs, the black dots denoting the (3,+1) CPs and the gray dots the (3,1) CPs. See text for details.
stress-tensor data, bond path lengths, bond path angles, and visualization of these properties. The XD2006 package also does the AIM analysis for the electron densities obtained from experimental results. Locating and characterizing all of the CPs of MED is an integral and important part of the AIM analysis. The studies on complete topography mapping of large molecules are rather scarce in the literature. In view of this, the present work is limited to mapping the complete topography of BNP, MED, and MESP of large molecules, without calculating the latter two functions on a grid. Carrying out ab initio calculations on large molecules is a difficult task, in particular, at higher levels of theory such as MP2. These methods require huge computational power and are timeconsuming. To overcome this difficulty, a state-of-the-art software called the molecular tailoring approach (MTA), is being developed in our laboratory since 1994.18 In the literature, some more active attempts based on a divide-and-conquer-type strategy have been reported. These include the DC method of Yang,19a GEBF approach from Li et al.,19b FMO from the group of Fedorov,19c and another DC method of the Nakai group,19d to name a few. In the present work, some representative molecular clusters and π-conjugated systems, which are of contemporary chemical interest and also exhibit rich topography, are chosen as test systems. A sequential scheme of topography mapping is presented, which begins by mapping the BNP topography, followed by MED. The rich MESP topography is subsequently explored, starting with the CPs of these two scalar fields. The last section presents the use of MTA for efficiently
capturing the topography of large molecular systems and assemblies. Ab initio calculations on such large systems are not routinely possible by the use of the conventional black box quantum chemical codes.
2. METHODOLOGY AND COMPUTATIONAL DETAILS For the present work, some weakly bound molecular clusters,20,21 namely, (CO2)5 , (C 6 H6)4 , (CH3 OH)3, (HCl)4 , and (H 2O)9 clusters, and few π-conjugated systems,22 namely, heptacene, annulene, C38H22, and β-carotene, are chosen as test systems. Further, the methodology for topography mapping of large molecules and clusters with the use of MTA is also tested out. The starting geometries of the test systems are optimized at the HF/6-31+G(d) level of theory. However, the first-order density matrix at the DFT/MP2 level with any other basis set could also be employed. These test systems are displayed23 in Figures 1 and 2, respectively. As a first step, the BNP topography is mapped for all of the test systems. For this purpose, BNP is calculated on a grid using eq 1. The present work does not exploit the symmetry of the systems because the method is intended for use with large molecular systems, which are generally devoid of symmetry. From the function values of BNP, guess points for CPs are generated. These guess points are then optimized using subroutine STEPIT.24 STEPIT is a direct-search, function-based optimizer for a function of many variables. It finds the local minima of a function by the method of successive relaxation. In the present context, the number of variables is three, namely, x, y, and z coordinates, which are continuously varied over a prescribed 12771
dx.doi.org/10.1021/jp2038976 |J. Phys. Chem. A 2011, 115, 12769–12779
The Journal of Physical Chemistry A
ARTICLE
Table 1. Critical Points of BNP for the Test Casesa n+3
n+1
n1
n3
total
(CO2)5
0
0
14
15
29
(C6H6)4
0
5
52
48
105
system
Table 2. Critical Points of MED (at the HF/6-31+G(d) Level) for the Test Casesa system
n+3
n+1
n1
n3
total
(CO2)5
0
2
16
15
33
5 0
16 1
58 18
48 18
127 37 17
(CH3OH)3
0
1
18
18
37
(HCl)4
0
1
8
8
17
(C6H6)4 (CH3OH)3
(H2O)9
0
2
28
27
57
(HCl)4
0
1
8
8
109
(H2O)9
1
6
31
27
65
heptacene
0
7
54
48
109
Heptacene
0
7
54
48
[18]annulene
0
1
36
36
73
C38H22 β-carotene
0 0
9 2
68 97
60 96
137 195
[18]annulene
0
4
39
36
79
C38H22
0
15
74
60
149
C32H24S
0
7
63
57
127
β-carotene
0
8
103
96
207
C32H24S
0
7
63
57
127
a
n+3, n+1, n1, and n3 are the numbers of different types of critical points. See text and Figures 1 and 2 for structures and details.
interval, starting from a guess point. The function to be minimized is chosen to be (∂f/∂x)2 + (∂f/∂x)2 + (∂f/∂x)2, where f may stand for BNP, MED, or MESP. The cutoff value for the minimum of the above function is set at 1016, ensuring that none of the individual Cartesian gradients numerically exceed 108. The steps in the procedure for mapping BNP topography are outlined below. Procedure I: Topography Mapping of BNP 1. Calculate the BNP of a molecule on a grid (typically 251 251 251). 2. Generate guess points from the grid and optimize them using STEPIT. 3. Characterize the CPs obtained by calculating the Hessian and check whether the PH relation is satisfied. 4. If the PH relation is not satisfied, recalculate the BNP with another grid and follow steps 2 and 3. Because the BNP calculation on a grid is computationally inexpensive, the function calculation with different grid sizes is affordable. The connection between the BNP and MED topography has been pointed out earlier.5,6,14 Keeping this in mind, in the next step, the MED topography is mapped by making use of the BNP topography. The MED topography of all of the test systems is obtained by using the density matrices (DM) at the HF/ 6-31+G(d) level.25 The studies pioneered by Bader et al.1,7 and several other works2629 reveal that the MED CPs lie within or quite close to the nuclear framework. Also, there exists a CP between a pair of nuclei that are bonded, that is, the (3,1) bond CP and a ring of nuclei that has a (3,+1) saddle. If the system contains a large network of CPs, mapping the complete MED topography of the system turns out to be a formidable task. Reports of such studies are indeed rather scarce in the literature. In view of this, the guess points for MED CPs for larger systems were generated by taking a few points on the line joining every pair of nuclei that falls within a 5 au distance. These points and BNP CPs are given as a guess to INDPROP30 for calculating MED CPs. The stepwise procedure for MED topography mapping for small/medium-sized molecules is outlined below. The STEPIT subroutine used in INDPROP optimizes the function over a given range of (x,y,z). There is a possibility that a guess point may not converge but gets closer to the CP. Such nonconverged guess points are again fed to INDPROP for optimization; the process is termed recycling in what follows. For the sake of brevity, the following terms are used in Procedures IIV discussed below.
a
Notations are identical to those used in Table 1. See text and Figures 3 and 4 for structures and details.
Search by Connecting CPs: Guess points for MED CPs are generated by taking a few points on a line joining the pairs of (3,+3) as well as those of (3,+1) CPs lying within a cutoff distance. These guess points are fed to INDPROP for optimization and characterization. The nonconverged points are recycled to ensure consistency in the number of CPs. Void Sphere Search: CPs that have no other CP lying within a sphere of a given cutoff radius are identified. Guess points are generated around the former at the corners of a cube of appropriate edge length and submitted for optimization. The nonconverged points are recycled to ensure consistency in the number of CPs. Procedure II: Topography Mapping of MED 1. Generate guess points for MED CPs by connecting all pairs of nuclei that lie within a cutoff distance. 2. Feed these points along with BNP CPs for optimization as a guess for MED CPs. Recycle the nonconverged points, to ensure consistency in the number of CPs. 3. Perform the step of Search by Connecting CP. 4. Perform the step of Void Sphere Search. 5. Check for the PH relation. If PH is not satisfied, repeat steps 3 and 4, with altered parameters. As noted above, the MESP topography is rich with CPs for lone pairs and π-bonds,31 which are conspicuous by their absence in the MED distribution. The earlier studies on MESP topography2,3,15 have shown that many of the negative-valued MESP CPs lie outside of the nuclear framework. For capturing MESP topography, the guess points for MESP CPs are hence generated at the corners and face centers of a cube around each nucleus at a typical 3 au distance. This distance is chosen due to the earlier observations of Gadre and Bhadane32 that the negative-valued CPs typically lie within 1.53 au from the nuclei of first-row atoms. These guess points along with MED CPs are given as input to INDPROP for calculating MESP CPs at the HF/6-31+G(d) level of theory. The detailed methodology for mapping the MESP topography of small/medium-sized molecules is as follows. Procedure III: Topography Mapping of MESP 1. Generate guess points for MESP CPs by setting up a cube around all the nuclei at 3 au distance. 2. Feed these points along with MED CPs for optimization. Recycle the nonconverged points to ensure consistency in the number of CPs. 3. Perform the step of Search by Connecting CP. 12772
dx.doi.org/10.1021/jp2038976 |J. Phys. Chem. A 2011, 115, 12769–12779
The Journal of Physical Chemistry A
ARTICLE
Figure 3. Critical points of the MED (at the HF/6-31+G(d) level) found for test systems, namely, (CO2)5, (C6H6)4, (CH3OH)3, (HCl)4, and (H2O)9 clusters. The black dots denote the (3,+1) CPs, the gray dots denote the (3,1) CPs, and the blue dots denote the (3,+3) CPs. See text for details.
4. Perform the step of Void Sphere Search. 5. Check for the PH relation. If PH is not satisfied, repeat steps 3 and 4, with altered parameters. It may be noted that the Euler characteristic χ does not bear a fixed value of 1 for MESP, but it is determined by examining asymptotic maxima and minima.15 After obtaining and characterizing the MESP CPs, the PH relation as specified in eq 7 is verified for all of the systems. Topography mapping of MED and MESP of large molecules is a formidable task. For this purpose, use of MTA is made in two different ways. First, MTA is employed for carrying out geometry optimization of large molecules, which cannot be readily treated by conventional black box packages, especially at high (correlated) levels of theory. Second, it is further used for evaluating the DM, which is subsequently used for mapping the topography of MED and MESP.33 The MTA method (details may be found in ref 18) is based on the set inclusion and exclusion principle, wherein one breaks the molecule into smaller, overlapping fragments. Ab initio calculations at the appropriate level are then effected for each of these fragments to get their DMs, which are suitably stitched together using eq 8 to form the synthesized DM for the whole molecular system. f i ∩ fj fi ∩ f j ... ∩ fk þ ... þ ð1Þk 1 Pab Pab ¼ Pfabi Pab
∑
∑
∑
ð8Þ
Here, Pab denotes the element (a,b) of the density matrix, Pfiab denoting the element (a,b) taken from the fragment i and so on. The MED and MESP are accurately computed from the DM that is assembled from the fragment DMs. The steps in topography mapping of MED and MESP of large molecules are outlined below. Procedure IV: Topography Mapping of MED for Large Molecules 1. Fragment the molecule into smaller fragments, keeping some overlap.
2. Map the BNP topography of individual main fragments by employing Procedure I. 3. Map the MED topography for each fragment by employing the fragment DMs and Procedure II described above. Check the PH relation for each fragment. 4. Calculate DM for the parent molecule from fragment DMs by employing MTA. 5. Patch all of the CPs from the fragments and feed them as a guess to the INDPROP package using the MTA-computed DM of the parent molecule. 6. Check for the PH relation. If not satisfied, carry out connecting the CP and Void Sphere Search. Procedure V: Topography Mapping of MESP for Large Molecules 1. Map the MESP topography for each fragment using Procedure III described above. Check the PH relation for each fragment. 2. Patch all of the CPs from the fragments and feed them as a guess to the INDPROP package using the MTA-computed DM of the whole molecule. 3. Check for the PH relation. If not satisfied, do Search by Connecting CP and Void Sphere Search. For the present work, β-carotene and (C6H6)9 at MP2/6-31 +G(d) and a (H2O)32 cluster at the HF/6-31+G(d) level of theory, respectively, are taken as test cases for the use of MTA in topography mapping.
3. TOPOGRAPHY OF TEST MOLECULES The above-described methodology is applied for mapping the BNP, MED, and MESP topography (cf. Figures 1 and 2) of small test systems at the HF/6-31+G(d) level of theory, and the results are discussed below. 3.1. Bare Nuclear Potential. The BNP topography is mapped for all of the test systems by the method described in Procedure I. 12773
dx.doi.org/10.1021/jp2038976 |J. Phys. Chem. A 2011, 115, 12769–12779
The Journal of Physical Chemistry A
ARTICLE
Figure 4. Critical points of the MED (at the HF/6-31+G(d) level) found for π-conjugated test systems, namely, heptacene, [18]annulene, C38H22, β-carotene, and C32H24S. The black dots denote the (3,+1) CPs, the gray dots the (3,1) CPs, and the blue dots the (3,+3) CPs. See text for details.
Table 3. Critical Points of MESP (at the HF/6-31+G(d) Level) for the Test Casesa n+3
n+1
n1
n3
(CO2)5
11
18
24
15
68
(C6H6)4
13
28
66
48
155
5 4
7 1
21 8
18 8
51 21
(H2O)9
18
30
39
27
104
heptacene
32
51
68
48
199
[18]annulene
18
22
41
36
117
C38H22
28
59
92
60
239
β-carotene
34
46
109
96
285
C32H24S
15
23
65
57
160
system
(CH3OH)3 (HCl)4
total
a
Notations are identical to those used in Table 1. See text and Figures 5 and 6 for structures and details.
The BNP CPs for all of the test systems are depicted in Figures 1 and 2, and the details of number of CPs of each kind are shown in Table 1. The search for nuclear maxima is not done because the nuclei are known to behave as pseudomaxima in the BNP distribution.1,7 First, let us summarize the results of the BNP topography of molecular clusters. For the case of (CO2)5, the total number of
BNP CPs is 24. Here, 4 (3,1) CPs are found in addition to the 10 bond CPs, indicating weak bonding present between the terminal O atom of a monomer and the C atom of another. For (C6H6)4, there appear five (3,+1) CPs, out of which four are centered in the individual benzene rings and one at the center of the cluster. The total numbers of BNP CPs for the case of (CH3OH)3, (HCl)4, and (H2O)9 clusters are 37, 17, and 28, respectively, with details as shown in Table 1. In π-conjugated systems (cf. Table 1), the first test case examined is from the acene family, namely, heptacene. For heptacene, 7 (3,+1) CPs are captured for the 7 rings along with 54 bond CPs. For C38H22, [18]annulene, C32H24S, and β-carotene, the number of CPs found is 68, 137, 127, and 195 respectively. Further details regarding their nature and number are reported in Table 1. It may be further noted that the PH relation in eq 6 holds true, that is, the value of χ is 1 for all of the cases. The CPU time required for the BNP topography, for all of the test systems, is less than 2 min on one core of a CPU with 8 GB of RAM. All of the CPU times reported in this work are reported in a similar fashion. 3.2. Molecular Electron Density. The MED topography of all of the test systems is calculated using INDPROP by employing the DM computed at the HF/6-31+G(d) level of theory employing Procedure II described in the above section. The numbers of 12774
dx.doi.org/10.1021/jp2038976 |J. Phys. Chem. A 2011, 115, 12769–12779
The Journal of Physical Chemistry A
ARTICLE
Figure 5. Critical points of the MESP (at the HF/6-31+G(d) level) found for test systems, namely, (CO2)5, (C6H6)4, (CH3OH)3, (HCl)4, and (H2O)9 clusters. The black dots denote the (3,+1) CPs, the gray dots the (3,1) CPs, and the blue dots the (3,+3) CPs. See text for details.
different types of MED CPs for all of the systems are reported in Table 2 and displayed in Figures 3 and 4. The nuclei of a molecule behave as a (3,3) maxima1,7 in the MED distribution; therefore, their search is not needed. For (CH3OH)3 as well as (HCl)4 clusters, the total number and nature of all of the CPs is identical to that in BNP. Thus, these two systems, namely, (CH3OH)3 and (HCl)4 clusters, show homeomorphism between the MED and BNP distribution as noted by Popelier14 for many of the systems studied by him. However, for the test cases of (CO2)5, (C6H6)4, and (H2O)9 clusters, there appear to be more CPs in MED as compared to those in the BNP topography. For the cases of heptacene and C32H24S, the number of MED CPs obtained is identical to the corresponding ones for BNP, as may be seen from Table 1. For [18]annulene, β-carotene, and C38H22, the number of MED CPs is more than the BNP CPs (cf. Table 1).
All of the MED CPs are captured by this method, including (3,+3) minima, which are conspicuous by their absence in the BNP distribution. For (CO2)5 and (CH3OH)3 clusters, the time required is less than 5 min, while for (C6H6)4, (H2O)9, heptacene, C38H22, and C32H24S, the time required is 49, 6, 19, 39, and 42 min, respectively. The number of MED CPs displayed in Table 2 also shows that the PH relation holds well for all of the systems, that is, the value of χ is ensured to be 1 for all of the cases. 3.3. Molecular Electrostatic Potential. As reported in the literature,2,3 MESP is endowed with rich topographical features. Thus, it is a daunting task to map the complete topography of MESP, especially for medium and large molecular systems. The MESP topography of all of the systems is mapped starting from the MED topography by employing Procedure III, as described in section 2. In MESP also, nuclear positions are 12775
dx.doi.org/10.1021/jp2038976 |J. Phys. Chem. A 2011, 115, 12769–12779
The Journal of Physical Chemistry A
ARTICLE
Figure 6. Critical points of the MESP (at the HF/6-31+G(d) level) found for π-conjugated test systems, namely, heptacene, [18]annulene, C38H22, β-carotene, and C32H24S. The black dots denote the (3,+1) CPs, the gray dots the (3,1) CPs, and the blue dots the (3,+3) CPs. See text for details.
Table 4. Critical Points of Three Scalar Fields for the Large Test Systems Employing the DM Generated with the Help of MTAa system
n+3
n+1
n1
n3
total
β-carotene
0
2
(C6H6)9
0
9
97
96
195
116
108
(H2O)32
0
0
95
233
96
191
β-carotene
0
8
(C6H6)9
11
38
103
96
207
134
108
(H2O)32
22
89
291
162
96
369
β-carotene (C6H6)9
34 31
46 64
109 142
96 108
285 345
(H2O)32
64
142
172
96
474
Bare Nuclear Potential
Molecular Electron Density
Molecular Electrostatic Potential
a
Notations are identical to those used in Table 1. See text and Figure 7 for structures and details.
pseudomaxima;2,3 hence, a search for nuclear maxima is not done. The number and types of MESP CPs are reported in Table 3 and are also displayed in Figures 5 and 6.
For the (CO2)5 cluster, a total of 68 MESP CPs are found. There are 11 (3,+3) minima, 18 (3,+1), and 24 (3,1) saddles. The MESP topography of (C6H6)4, (HCl)4, and (H2O)9 clusters is seen to be comprised of 107, 21, and 104 CPs in all, the details of which are given in Table 3. For the heptacene, 32 (3,+3) minima, 51 (3,+1), and 68 (3,1) saddle CPs are found in the MESP distribution. For [18]annulene, there are 18 (3,+3) minima, 22 (3,+1), and 41 (3,1) saddles found. In the MESP topography of C38H22, β-carotene, and C32H24S, the total number of CPs turns out to be 149, 187, and 102, respectively, as per the details given in Table 3. In this way, the MESP topography for all of the test systems is built by using the MED topography as a guess, supplemented by the cube search as described in Procedure III. All of the CPs are captured by this method, including (3,+3) minima for lone pairs of atoms and π-bonds. For finding the values of asymptotic maxima (n) and minima (n+), MESP for a given test system is generated on a large sphere. This sphere is then visualized in the UNIVIS package, and the values of n, n+, and hence χ are obtained. The correctness of the MESP topography is checked by finding asymptotic minima and maxima and then applying a PH relation in eq 7 to all of the systems. The relation holds true for the MESP topography of all test systems. The time required for the MESP topography calculation of (CO2)5, (C6H6)4, and (H2O)9 clusters is approximately 18, 365, and 26 min, respectively. For conjugated systems, heptacene, [18]annulene, 12776
dx.doi.org/10.1021/jp2038976 |J. Phys. Chem. A 2011, 115, 12769–12779
The Journal of Physical Chemistry A
ARTICLE
Figure 7. Critical points of (A) the MED and (B) the MESP found for β-carotene (left) and the (C6H6)9 cluster (right) at the MP2/6-31+G(d) level of theory. The black dots denote the (3,+1) CPs, the gray dots the (3,1) CPs, and the blue dots the (3,+3) CPs. See text for details.
Figure 8. Critical points of (A) the MED and (B) the MESP found for the (H2O)32 cluster at the HF/6-31+G(d) level of theory. The black dots denote the (3,+1) CPs, the gray dots the (3,1) CPs, and the blue dots the (3,+3) CPs. See text for details.
C38H22, and C32H24S, the average time required is somewhat larger, with the individual CPU times being 202, 84, 350, and 359 min, respectively. It may be noted that the CPU time required for function and gradient evaluation on a large grid is generally much larger, a step that is avoided in the present methodology.
4. TOPOGRAPHY OF LARGE MOLECULES Here, results of MTA applied for calculating the topography of three large systems viz. β-carotene and (C6H6)9 at MP2/6-31 +G(d) and (H2O)32 at HF/6-31+G(d) level of theory are presented. As stated in Section 2 above, use of MTA is made for calculating the DM of the large systems by applying eq 8.
For demonstrating the mapping of the topography of large molecules, β-carotene optimized at the HF level is chosen as a test systems at the MP2/6-31+G(d) level of theory, employing the geometry optimized at the HF/6-31+G(d) level. The MTAoptimized geometries of cluster systems (C6H6)9 and (H2O)32 at the MP2/6-31+G(d) and HF/6-31+G(d) levels of theory, respectively, are used for the topography mapping. The MED and MESP topography of these large systems is mapped by employing Procedures IV and V, respectively. as described in section 2 above. The nature and numbers of MED and MESP CPs obtained for these molecules are reported in Table 4 and are displayed in Figures 7 and 8. For the MED topography, the PH relation in eq 6 holds true, that is, the Euler characteristic is 12777
dx.doi.org/10.1021/jp2038976 |J. Phys. Chem. A 2011, 115, 12769–12779
The Journal of Physical Chemistry A χ = 1, and for MESP, the PH relation in eq 7 is found to be satisfied for all three systems, with an appropriate value of χ. As can be seen from Table 4, the β-carotene molecule and the cluster systems (C6H6)9 and (H2O)32 exhibit a large number of MED CPs, and the corresponding MESP is even richer in the number of CPs for all three systems. For the test case of β-carotene at the HF/6-31+G(d) level (cf. section 3.3), the five lowest (3,+3) MESP minima lie between 0.0332 and 0.0322 au. Although there are small changes in the minimum values at the MP2/6-31+G(d) level of theory, the CP coordinates are almost identical to those at the HF level. This is in agreement with the earlier work reporting the effect of the basis set and level of theory on the MESP topography.34
5. CONCLUDING REMARKS The path-breaking works due to Bader and co-workers1 have generated widespread interest in the topography of molecular scalar fields. In the literature, there are two main approaches for mapping the topography of scalar fields such as MED and MESP. The first of these methods is based on a detailed grid-based search, wherein the function values are enumerated at a large number of grid points, followed by a numerical interpolation for locating the CPs. The other one is to follow gradient paths starting from the set of CPs that have already been located. Both of these alternatives are computationally expensive for medium as well as large molecules, especially when treated at higher levels of theory. As mentioned earlier, the work of Popelier et al.14 reports the topography mapping of BNP and MED for many systems. However, reports on the mapping of all CPs of MED and MESP are conspicuous by their absence from the literature studies. The present work has demonstrated, through the study of chosen test molecules, that the critical points could be economically located, starting with the scalar field of BNP, followed by MED and MESP. After a thorough search, a check on the number of CPs of each type is made in terms of the PH relation. An earlier study has revealed that the nature and number of CPs in MESP is independent of the level of theory and basis set if a sufficiently large basis (incorporating diffuse/polarization functions) is employed.34 The building up of the topography of a molecule from the scalar field of the BNP to MED and then to MESP seems a natural choice because the topography successively becomes richer along this path. Thus, the present approach, enabled by MTA, paves the way for mapping the topography of large molecules. For example, the cluster systems of (CO2)5 and (C6H6)4 are endowed with 29 and 105 BNP CPs, respectively. These numbers grow to 33 and 127 for the MED distribution. For the scalar field of MESP, the respective numbers are even larger, namely, 68 and 155, respectively. This trend is observed for conjugated systems also, divulging the sparse nature of BNP topography. The MED topography is enriched further with the addition of cage minima as well as some other CPs. This makes the MED topography, in general, richer than the BNP topography. The feature of MESP to attain both negative as well as positive values makes its topography even more exuberant in terms of the number of CPs. For large molecular systems, mapping of the MED and MESP topography is indeed a formidable task because even the geometry optimization of such systems cannot be readily done, especially at correlated levels of theory, due to high-order scaling of the methods. The benzene nonamer at the MP2/6-31+G(d) level of theory offers an example of success of MTA, not only for enabling
ARTICLE
the geometry optimization but also for obtaining the DM of this system for calculating MED and MESP distributions. It is thus hoped that the present fragment-based MTA study paves the way for mapping the MED and MESP topography of large systems, including biomolecules and molecular aggregates. Our earlier work has shown that molecules recognize each other at large intermolecular separations in terms of the topography of MESP.35 Thus, the present approach can be used for studying molecular recognition in large chemical and biological systems through the topography of MESP.
’ AUTHOR INFORMATION Corresponding Author
*E-mail:
[email protected].
’ ACKNOWLEDGMENT Authors thank the Center for Development of Advanced Computing (C-DAC), Pune, for providing financial and computational support. S.R.G. is thankful to the Department of Science and Technology (DST) for the award of the J. C. Bose National Fellowship. ’ REFERENCES (1) Bader, R. F. W. Atoms in Molecules: A Quantum Theory; Clarendon: Oxford, U.K., 1990. (2) Gadre, S. R. In Computational Chemistry: Reviews of Current Trends; Lesczynski, J., Ed.; World Scientific: Singapore, 2000; Vol. 4. (3) Gadre, S. R.; Shirsat, R. N. In Electrostatics of Atoms and Molecules; Universities Press: Hyderabad, India, 2000. (4) Hohenberg, P.; Kohn, W. Phys. Rev. 1964, 136, B864–B871. (5) Parr, R. G.; Berk, A. In Chemical Applications of Atomic and Molecular Electrostatic Potentials; Politzer, P., Truhlar, D. G., Eds.; Plenum: New York, 1981; pp 5162. (6) Parr, R. G.; Gadre, S. R.; Bartolotti, L. J. Proc. Natl. Acad. Sci. U.S. A. 1979, 76, 2522–2526. (7) (a) Bader, R. F. W.; Nguyen-Dang, T. T.; Tal, Y. J. Chem. Phys. 1979, 70, 4316–4329. (b) Malcolm, N. O. J.; Popelier, P. L. A. J. Comput. Chem. 2003, 24, 437–442. (c) Popelier, P. L. A. Coord. Chem. Rev. 2000, 197, 169–189. (8) Weinstein, H.; Politzer, P.; Srebrenik, S. Theor. Chim. Acta 1975, 38, 159–163. (9) Pathak, R. K.; Gadre, S. R. J. Chem. Phys. 1990, 93, 1770–1773. (10) (a) Cao, W. L.; Gatti, C; McDougall, P. J.; Bader, R. F. W. Chem. Phys. Lett. 1987, 141, 380–385. (b) Cioslowski, J. J. Phys. Chem. 1990, 94, 5496–5498. (c) Wiberg, K. B.; Bader, R. F. W.; Lau, C. D. H. J. Am. Chem. Soc. 1987, 109, 985–1001. (11) Tal, Y.; Bader, R. F. W.; Erkku, J. Phys. Rev. A 1980, 21, 1–11. (12) Politzer, P.; Zilles, B. Croat. Chem. Acta 1984, 57, 1055–1064. (13) Gadre, S. R.; Bendale, R. D. Chem. Phys. Lett. 1986, 130, 515–519. (14) Popelier, P. L. A.; Bremond, E. A. G. Int. J. Quantum Chem. 2009, 109, 2542–2553. (15) (a) Leboeuf, M.; K€oster, A. M.; Jug, K.; Salahub, D. R. J. Chem. Phys. 1999, 111, 4893–4905. (b) Roy, D. K.; Balanarayan, P.; Gadre, S R. J. Chem. Phys. 2008, 129, 174103–174108. (16) (a) Keith, T. A.; Gristmill, T. K. AIMAll, version 11.06.19; Software: Overland Park, KS, 2011; aim.tkgristmill.com. (b) Popelier, P. L. A. Comput. Phys. Commun. 1995, 93, 212–240.(c) Volkov, A.; Macchi, P.; Farrugia, L. J.; Gatti, C.; Mallinson, P.; Richter, T.; Koritsanszky, T. XD2006 Program; 2006. (17) (a) Volkov, A.; Gatti, C.; Abramov, Yu.; Coppens, P. TOPXD Program. Acta Cryst. A. 2000, 56, 252.(b) Krug, P. Program BUBBLE; McMaster University: Hamilton, Ontario, 1990. 12778
dx.doi.org/10.1021/jp2038976 |J. Phys. Chem. A 2011, 115, 12769–12779
The Journal of Physical Chemistry A (18) Ganesh, V.; Dongare, R. K.; Balanarayan, P.; Gadre, S. R. J. Chem. Phys. 2006, 125, 104109–104118. (19) (a) Yang, W. Phys. Rev. A 1991, 44, 7823–7826. (b) Yang, Z.; Hua, S.; Hua, W.; Li, S. J. Phys. Chem. A 2010, 114, 9253–9261. (c) Fedorov, D G; Ishida, T; Uebayasi, M; Kitaura, K J. Phys. Chem. A 2007, 111, 2722–2732. (d) Akama, T; Kobayashi, M; Nakai, H J. Comput. Chem. 2007, 28, 2003–2012. (e) Bettens, R. P. A.; Lee, A. M. J. Phys. Chem. A 2006, 110, 8777–8785. (f) Collins, M. A. J. Chem. Phys. 2007, 127, 24104–10. (20) (a) Rahalkar, A. P.; Katouda, M.; Gadre, S. R.; Nagase, S. J. Comput. Chem. 2010, 31, 2405–2418. (b) Mahadevi, A. S.; Rahalkar, A. P.; Gadre, S. R.; Sastry, G. N. J. Chem. Phys. 2010, 133, 164308–164319. (21) Yeole, S. D.; Gadre, S. R. J. Chem. Phys. 2011, 134, 084111–084119. (22) (a) Yeole, S. D.; Gadre, S. R. J. Chem. Phys. 2010, 132, 094102–094109. (b) Nguyen, H. T.; Truong, T. N. Chem. Phys. Lett. 2010, 499, 263–267. (23) UNIVIS-2000: a comprehensive visualization package; see: Limaye, A. C.; Gadre, S. R. Curr. Sci. 2001, 80, 1296–1301. (b) The package MeTA Studio: Ganesh, V. J. Comput. Chem. 2009, 30, 661–672. (24) Chandler, J. P. The subroutine STEPIT; distributed by QCPE, University of Indiana: Bloomington, IN, 1967. (25) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A., GAUSSIAN03, revision B.04; Gaussian Inc.: Pittsburgh, PA, 2003. (26) (a) Bader, R. F. W.; Slee, T. S.; Cremer, D.; Kraka, E. J. Am. Chem. Soc. 1983, 105, 5061–5068. (b) Cremer, D.; Kraka, E.; See, T. S.; Bader, R. F. W.; Lau, C. D. H.; Nguyen-Dang, T. T.; MacDougall, P. J. J. Am. Chem. Soc. 1983, 105, 5069–5075. (c) Cremer, D.; Kraka, E. Croat. Chem. Acta 1984, 57, 1265–1287. (27) (a) Matta, C. F.; Hernandez-Trujillo, J.; Tang, T.-H.; Bader, R. F. W. Chem.—Eur. J. 2003, 9, 1940–1951. (b) Espinosa, E.; Molins, E.; Lecomte, C. Chem. Phys. Lett. 1998, 285, 170–173. (c) Baranov, A. I.; Kohout, M. J. Phys. Chem. Solids 2010, 71, 1350–56. (28) (a) Cremer, D; Kraka, E. J. Am. Chem. Soc. 1985, 107, 3800–3810. (b) Gatti, C Z. Kristallogr. 2005, 220, 399–427. (c) Espinosa, E.; Molins, E. J. Chem. Phys. 2000, 113, 5686–5694. (29) (a) Matta, C. F.; Castillo, N.; Boyd, R. J. J. Phys. Chem. B 2006, 110, 563–578. (b) Bader, R. F. W.; Matta, C. F. Int. J. Quantum Chem. 2001, 85, 592–607.(c) Matta, C. F., Boyd, R. J., Eds. The Quantum Theory of Atoms in Molecules: From Solid State to DNA and Drug Design; Wiley-VCH: Weinheim, Germany, 2007. (30) INDPROP, the molecular property calculation package developed at Theoretical Chemistry Group, Department of Chemistry, University of Pune, Pune, India. See: Bapat, S. V.; Shirsat, R. N.; Gadre, S. R. Chem. Phys. Lett. 1992, 200, 373–378. (31) (a) Gadre, S. R.; Pathak, R. K. Proc.—Indian Acad. Sci., Chem. Sci. 1990, 102, 189–192. (b) Gadre, S. R.; Kulkarni, S. A.; Shrivastava, I. H. J. Chem. Phys. 1992, 96, 5253–5260. (c) Kulkarni, S. A.; Gadre, S. R. J. Mol. Struct.: THEOCHEM 1996, 361, 83–91. (d) Luque, F. J.; Orozco, M.; Bhadane, P. K.; Gadre, S. R. J. Phys. Chem. 1993, 97, 9380–9384. (32) Gadre, S. R.; Bhadane, P. K. J. Chem. Phys. 1997, 107, 5625–5626. (33) (a) Gadre, S. R.; Shirsat, R. N.; Limaye, A. C. J. Phys. Chem. 1994, 98, 9165–9169. (b) Babu, K; Gadre, S. R. J. Comput. Chem. 2003,
ARTICLE
24, 484–495. (c) Babu, K; Ganesh, V.; Gadre, S. R.; Ghermani, N. E. Theor. Chem. Acc. 2004, 111, 255–263. (34) Gadre, S. R.; Kulkarni, S. A.; Suresh, C. H.; Shrivastava, I. H. Chem. Phys. Lett. 1995, 239, 273–281. (35) Roy, D. K.; Balanarayan, P.; Gadre, S R. J. Chem. Sci. 2009, 121, 815–821.
12779
dx.doi.org/10.1021/jp2038976 |J. Phys. Chem. A 2011, 115, 12769–12779
