Density Functional Theory of Molecular Solids: Local versus Periodic

Rubicelia Vargas, Jorge Garza, Richard A. Friesner, Harry Stern, Benjamin P. Hay, and David A. Dixon .... Wen-Ge Han , Emadeddin Tajkhorshid , Sándor...
1 downloads 0 Views 557KB Size
Download PDF
3950

J. Phys. Chem. 1996, 100, 3950-3958

Density Functional Theory of Molecular Solids: Local versus Periodic Effects in the Two-Dimensional Infinite Hydrogen-Bonded Sheet of Formamide Sa´ ndor Suhai Molecular Biophysics Department, German Cancer Research Center, Im Neuenheimer Feld 280, D-69120 Heidelberg, Germany ReceiVed: September 8, 1995; In Final Form: NoVember 22, 1995X

The performance of an ab initio computational scheme for molecular crystals based on density functional theory (DFT) was investigated by computing several structural and energetic properties of hydrogen bonded infinite chains and of two-dimensional infinite periodic networks of formamide. The applied DFT potentials covered a wide range in quality, starting with a simple local exchange (X) without correlation (C), and then gradually introducing C and gradient corrections for both X and C. At the same time, five atomic basis sets of systematically increasing size, in the range of DZ to TZ(2df,2pd) were used to construct the Bloch-type crystal orbitals, to optimize the structures, and to extrapolate different physical quantities to the limit of a hypothetical infinite basis set. Infinite lattice sums were computed by the multipole expansion technique, and basis set superposition errors were (partly) eliminated by the counterpoise method. To be able to assess the accuracy of the theoretical models, the formamide monomer and two different dimers were also investigated using the same methods. Detailed comparisons were made for all models also with recent results obtained by using different orders of many-body perturbation theory. Structural optimizations for the dimers and the infinite crystal demonstrated the importance of gradient terms both for exchange and correlation. For the most successful DFT functional, containing the Becke exchange and the Lee-Yang-Parr correlation term, the lengths of the hydrogen bonds, RHB, were reduced by 0.16-0.19 Å (depending on the basis set) as compared with dimers, due to cooperative interactions in the crystalline environment. The binding energies were increased typically by 60-70 %. The theoretical model explained why the RHB values for open-chain dimers become shorter in the crystal than those obtained for the cyclic ones (as opposed to free dimers), correctly predicted changes of bong lengths in going from the monomer to the crystal, and provided NsH‚‚‚O bond distances and lattice constants reasonably close to experiments.

I. Introduction Hydrogen bonded crystals provide a large group of interesting materials to test the accuracy of theoretical models in predicting structural and energetic properties of molecular solids. The strength of intermolecular forces that act among their building blocks, especially with regard to nonadditive interactions, is also a critical element in our understanding of the nature of binding in several important biological macromolecules. Formamide provides one example for such crystals. Its characteristic structural element is the NsH‚‚‚OdC hydrogen bond that plays a crucial role not only in several other molecular solids but also forms the basic structural feature in nucleic acids and proteins. In the solid state, formamide is composed of puckered layers (Figure 1), each of which can be described either as being built from alternating (zig-zag) chains or from cyclic dimers of the HCONH2 unit.1 The hydrogen bonds within the zig-zag chains should be stronger than within the cyclic dimers, as reflected by the NsH‚‚‚O bond length of 2.88 Å in the former case as compared to 2.93 Å for the latter one. The results reported in this paper for the formamide crystal form a part of the long-term research goal to develop first principles quantum mechanical methods for the computer-aided prediction of the geometrical and electronic structure of polymers and molecular crystals.2-10 To be able to describe the molecular details of the interactions at a given lattice site, but also to properly take into account the collective response of the system, such methods have to combine theoretical and computational procedures of both molecular and solid state X

Abstract published in AdVance ACS Abstracts, February 1, 1996.

0022-3654/96/20100-3950$12.00/0

Figure 1. Structure of the unit cell in the infinite two-dimensional model of the formamide crystal.

physics. Density functional theory (DFT), that not only had a long history in solid state physics11-14 but recently also successfully penetrated the molecular field,15-18 provides a most reasonable ansatz for this goal. Different versions of DFT were also applied to polymers.19-25 The conventional method to calculate the electronic properties of molecular crystals is to first perform a Hartree-Fock (HF) calculation26,27 followed by a subsequent computation of electron correlation effects using different orders of many-body perturbation theory (MBPT).2,7,9 In the case of a careful Coulomb integral handling and with the use of proper numerical schemes in performing the infinite © 1996 American Chemical Society

Density Functional Theory of Molecular Solids

J. Phys. Chem., Vol. 100, No. 10, 1996 3951

lattice sums,5 the first step scales essentially as O(N2-N2.5) for a polymer, where N is the number of atomic basis functions in the (translationally invariant) unit cell. The computational expense of MBPT procedures, on the other hand, starts formally with O(N5) in second order and increases rapidly for higher orders. Though the systematic use of localized Wannier orbitals and of advanced integral transformation and selection procedures reduces this burden, such calculations are clearly too demanding to become routine within the next years. For the investigation of physical and chemical properties that do not explicitly depend on the wave function itself, DFT may provide an alternative way to estimate the correlation effects at a relatively modest cost, formally O(N3), which may be reduced by an efficient implementation.28 Therefore, DFT polymer studies including electron correlation will scale only again as O(N2-N2.5), similar to the HF case. Hence, the quality of DFT methods as applied to molecular crystals deserves attention and careful comparison, both with experimentally observable quantities and with the results of MBPT calculations using extended atomic basis sets. Though the main object of our present study is the twodimensional, infinite periodic lattice of formamide molecules (Figure 1), we will also compare several properties of this system with those of the monomer and of two dimers of formamide. On the one hand, such comparisons will demonstrate how certain properties of the solid differ from the corresponding dimer interactions. On the other hand, the monomer and the dimers will serve as benchmarks characterizing the computational schemes. Before turning to the infinite solid, we will perform systematic structure optimizations and interaction energy calculations for them and compare the results with experiments and with the most sophisticated theoretical investigations. This procedure should enable us to further evaluate the accuracy of the DFT results obtained for the infinite system. In particular, we want (i) to understand the structural and energetic differences in the binding of two formamide molecules within isolated dimers and embedded into a crystalline environment, respectively, (ii) to investigate the role of electron correlation effects in those interactions, and (iii) to identify the exchangecorrelation potential(s) that provide(s) the best description of the structural and energetic properties in comparison with highaccuracy MBPT results and with experiments, respectively.

The ab initio method to investigate the electronic structure of infinite polymers using HF and various orders of MBPT was described in detail.2,7,9 Since in extending the HF crystal orbital procedure to DFT we have to perform the same steps as those used when introducing DFT to molecular systems,16 we only summarize here the basic expressions of the LCAO Bloch crystal orbital (CO) ansatz that substitutes the MO picture in the case of extended (infinite) systems with periodic boundary conditions. In fact, our computational procedure is a generalization of the methods worked out by Pople’s group over the past decades38-41 for molecules, to the case of infinite systems using additional symmetry operations and different numerical procedures of solid state theory for the manipulations in the Brillouin zone (BZ). The computer programs developed in our laboratory for polymer electronic structure calculations have also made intensive use of several techniques and modules of the Gaussian molecular packages starting from Gaussian 76 up to Gaussian 92/DFT.41 Within the traditional HF approach, in the case of the spinrestricted (RHF) theory the Fock operator of the twodimensional crystal consists of kinetic energy, nuclear contributions, Coulomb, and exchange terms in the form

II. LCAO Density Functional Approach for Molecular Solids

r l), themselves will be obtained as selfThe Bloch orbitals, φBnk (b consistent solutions of the HF equations of the crystal:25,26

In their local spin density (LSD) form, DFT methods simply replace the HF exchange functional with an exchange-correlation (XC) functional of the local electron spin densities FR and Fβ.13-16 Gradient corrected methods include also the spin density gradients ∇FR and ∇FR to account for local fluctuations in the electron density. The Hohenberg-Kohn theorem11 ensures that there exists a unique functional which exactly determines the ground state energy. Since the exact form of the functional is unknown, several approximate forms of the XC potential have been proposed.29-37 In the most widely used formulation of DFT, the Kohn-Sham (KS) theory,12 an LCAOtype basis set expansion of one-electron orbitals is obtained via a self-consistent procedure analogous to that of conventional HF theory. This allows correlation to be included at the SCF level at some nominal extra cost. An alternative to computing the KS density is to obtain the electron density via the oneelectron orbitals resulting from a HF calculation and substituting it into the LSD expression. This hybrid method involves performing a single numerical integration after solving the HF problem and may be even more cost effective than the KS procedure.

1 Fˆ (b r l) ) Tˆ + Vˆ + Jˆ + K ˆ ) - ∆br l 2 )Nc/2 +Nc/2 nA F(b r l, b r l) ZA 3 + d b r ∫ ∑ ∑ ∑ r - BR - BR m |b rl - b r m| j)-Nc/2 h)-Nc/2 A)1 |b l jh A F(b r l,b r m)

∫d3br m |br - br l

r m) (1) Pˆ (b r l,b

m|

where (Nc + 1)2 is the number of elementary cells in the square lattice, nA is the number of atoms per cell, and ZA is the core charge of the atom A at position B RA. The first-order density bm), will be constructed from (doubly) occupied matrix, F(r bl,r Bloch spin orbitals by numerical integration over the first BZ as occ. {BZ}

r m) ) ∑ ∑ φBnk (b r l)[φBnk (b r m)]* F(b r l, b n

(2)

B k

r l) ) Bnk φBnk (b r l) Fˆ (b r l) φBnk (b

(3)

As mentioned before, in a first approximation we will model the hydrogen-bonded network of formamide molecules in the crystal by a two-dimensional infinite periodic array of HCONH2 units and neglect interactions between the sheets. Our computational unit cell will thus consist of four formula units. For segregated (one-dimensional) chains, due to the presence of a screw symmetry (translation along the chain axis plus a rotation around it by 180°), the unit cell will be even reduced to one formula unit. The substantially decreased size of the computational problem will allow, consequently, much more accurate calculations in this case. According to the translational symmetry, the one-electron Bloch functions can be identified with the quasi-momentum B k and the band index n. The values of the reciprocal space vectors B k(kx,ky) of the two-dimensional rectangular lattice will be restricted during the computations for the irreducible part of the first BZ. The Bloch functions will be expressed as linear combinations of symmetry adapted Bloch basis orbitals:

3952 J. Phys. Chem., Vol. 100, No. 10, 1996

Suhai

ν

k φBnk (b) r ) ∑cBan ψBak (b) r

(4)

a)1

r ψabr (b),

per elemenwhere ν is the number of basis functions, tary cell of the crystal, which themselves will be linear combinations of contracted Gaussian-type atomic orbitals (CGTOs): +Nc/2

r ) (Nc + 1)-1 ψBak (b)

+Nc/2

exp{ik BB Rjh}χjh r ∑ ∑ a (b) j)-N /2 h)-N /2 c

(5)

c

r ) χa(r b- B Ra - B Rjh) will be centered The atomic orbital χjh a (b) Rjh. The N-electron in the cell with the index pair [j,h], at B Ra + B wave function ΦHF will be written as a Slater determinant built from the doubly filled symmetry adapted Bloch functions in the form

r i) R(σi) φBnk (b r i + 1) β(σi+1)...] (6) ΦHF ) (N!)-1/2det[...φBnk (b For the DFT formulation of the same problem, the Fock operator takes the form analogous to eq 1 with the difference that the HF exchange will be substituted with the appropriate exchangecorrelation term:

Fˆ DFT ) Tˆ + Vˆ + Jˆ + Fˆ XC

(7)

and the crystal HF equations will go over to the corresponding Kohn-Sham equations.12 The self-consistent and gradient procedures for solids are again in complete analogy to the molecular case, and the appropriate moduls of the G92/DFT program package41 were utilized, after extension for translational and helical symmetry operations, for their solution. Several computational aspects of this procedure specific to infinite systems (handling of infinite lattice sums, numerical accuracy of Brillouin zone manipulations, etc.) will be reported elsewhere.42 The functionals used in our computations consists of separate exchange and correlation parts, respectively. For the exchange part, we will use either the free-electron gas functional proposed by Slater (S)29 or the gradient-corrected functionals designed by Becke, B 35 and B3,43 respectively, to correctly reproduce the exact asymptotic behavior of the exchange energy density for a finite many-electron system. Either the correlation part will be ignored (leading to the HFS and HFB theories, respectively) or it will be treated by the local spin density theory as parametrized by Vosko, Wilk, and Nusair (VWN)31 or by the gradient-corrected functional of Lee, Yang, and Parr (LYP)36 as transformed by Mielich et al.37 While the VWN functional reproduces exact uniform electron gas results, the LYP functional was designed (initially by Colle and Salvetti30) to calculate the exact nonrelativistic energy of the He atom from the HF density. Both correlation functionals can be combined with exchange terms providing the potential schemes SVWN, SLYP, BVWN, BLYP, and B3LYP, respectively, as proposed by Johnson et al.44 The LYP functional can also be combined with the exact exchange leading to the HFLYP method. As it is well-known, the intermolecular interaction energy of the formamide dimer, ∆E, may sensitively depend on the extension of the atomic basis set applied and it contains, especially in the case of basis sets of limited size, a considerable error called basis set superposition error, BSSE. In the crystalline environment, we may expect this error to further increase partly as a consequence of shorter hydrogen bond distances in the solid (as compared to the dimer due to a cooperative attraction in the extended system) and partly because

TABLE 1: List of the Applied Sets of Uncontracted and Contracted Atomic Basis Functions, Respectively uncontracted basis components Na

symbol of basis set DZ DZ(d,p) DZ(2d,2p) TZ(2d,2p) TZ(2df,2pd)

(9s5p/4s) (9s5p1d/4s1p) (9s5p2d/4s2p) (11s6p2d1f/6s2p) (11s6p2d1f/6s2p1d)

84 111 138 153 189

contracted basis components Na

refb

[4s2p/2s] 36 55 [4s2p1d/2s1p] 63 55,41 [4s2p2d/2s2p] 90 55,41 [5s3p2d1f/4s2p] 102 56,57 [5s3p2d1f/4s2p1d] 138 56,57

a

The total number of uncontracted or contracted atomic basis functions, respectively, resulting by using the given basis set for the formamide monomer. b The first reference refers to orbitals in the basis set which are occupied in the ground state of the corresponding atom; the second one, to the polarization functions used.

further neighboring molecules provide their basis functions for use of the monomer in the reference cell. To at least partially correct for BSSE, we will be using the counterpoise correction (CPC) scheme as proposed by Jansen and Rose45 and by Boys and Bernardi.46 If EA(A) denotes the energy of the water monomer A calculated with its own (monomer) basis functions at the geometry optimized for the isolated monomer, and EA(Ab) indicates the use of basis A plus the ghost basis on monomer B according to the geometry optimized for the dimer, then the interaction energy of the dimer containing BSSE will be given by

∆E(BSSE) ) E(A,B) - EA(A) - EB(B)

(8)

where E(A,B) stand for the computed total energy of the dimer. In a similar manner, the CPC interaction energy will be obtained as

∆E(CPC) ) E(A,B) - EA(Ab) - EB(aB)

(9)

It is not obvious, however, that the relaxation of the monomers upon formation of the dimer will be negligible for this interaction. Correction terms complementing the Boys-Bernardi scheme to take into account this effect have been proposed by Mayer and Surjan.47 Denoting the energy of molecule A computed at the relaxed (dimeric) geometry using the monomeric basis set by ErA(A), the interaction energy corrected for monomer relaxation can be obtained in the form

∆Er(CPC) ) E(A,B) - EA(Ab) - EB(aB) + [EAr (A) EA(A)] + [EBr (B) - EB(B)] (10) This will be our working expression to correct the BSSE both for the dimer and the polymer in the next sections. III. Monomer and Dimer of Formamide The structure and properties of the formamide molecule were the subject of numerous experimental and theoretical studies. Here, we can cite only a few recent papers48-54 which contain references to previous works. To calibrate our theoretical models, we investigated the structure of the monomer using atomic basis sets of systematically increasing size as defined in Table 1. The first part of the basis sets refers to orbitals occupied in the ground state of the corresponding atoms. The exponents and contraction coefficients of these orbitals have been taken from the double-ζ (DZ) and triple (TZ) schemes of Dunning.55,56 The exponents of the polarization functions used to complement the isotropic part originate partly from the library of the Gaussian program packages41 and partly from the correlation optimized sets of Dunning.57 The structure of the monomer was optimized by assuming planarity. The nonpla-

Density Functional Theory of Molecular Solids

J. Phys. Chem., Vol. 100, No. 10, 1996 3953

TABLE 2: Molecular Total Energy (hartrees), Bond Lengths (Å), and Dipole Moment (D) of the Formamide Molecule Computed with Several DFT Potentials Using the TZ(2df,2pd) Atomic Basis Set Defined in Table 1a method

Etot.

RC,O

RC,N

RN,H

HFS HFB SVWN SLYP BVWN BLYP B3LYP HFLYP HF MP2 MP4 exptl (ED)b exptl (MW)c

-167.045 974 -169.108 617 -169.097 938 -167.870 036 -171.163 637 -169.931 880 -169.973 282 -169.837 338 -169.007 586 -169.660 554 -169.704 516

1.2227 1.2344 1.2091 1.2108 1.2204 1.2220 1.2095 1.1770 1.1865 1.2133 1.2179 1.2110 1.2190

1.3653 1.3905 1.3459 1.3454 1.3696 1.3690 1.3572 1.3283 1.3460 1.3550 1.3600 1.3670 1.3520

1.0329 1.0243 1.0168 1.0232 1.0086 1.0148 1.0064 0.9841 0.9885 0.9998 1.0014 1.0210 1.0020

µ 3.8875 3.8134 3.9652 3.9461 3.8920 3.9006 3.9645 4.1779 4.1353 3.8808 -

a For comparison, HF and MBPT results from ref 10 are also shown together with two experimental data sets. b Reference 59. c Reference 60.

Figure 3. RC,N bond length in the formamide molecule obtained using different methods with the five atomic basis sets defined in Table 1.

Figure 2. RC,O bond length in the formamide molecule obtained using different methods with the five atomic basis sets defined in Table 1.

narity of formamide has been discussed in detail previously.49,58 Though MP2 calculations led to a slight energetic preference of a pyramidal amino group, the energy difference between the C1 and Cs structures is only a few microhartrees and it will not play a role in our present investigations. The molecular total energy, the most important bond lengths, and the dipole moment of the monomer obtained with different methods using the TZ(2df,2pd) basis set are collected in Table 2 together with two experimental data sets. Apart from the HFB and HFLYP methods, the overall agreement of the DFT bond lengths with the MP results and with experiments is surprisingly good. This observation is also in agreement with the DFT results obtained recently by Florian and Johnson for the formamide monomer.52 Obviously, we shall have to look at more difficult problems (like the description of the intermolecular H-bond in the dimers) to be able to better differentiate between the various DFT potentials. It should be mentioned, however, that the quality of the above results sensitively depends on the extension of the basis set. For two potentials, BLYP and B3LYP, which will play a major role for the solid, the strong dependence of two bond lengths on the number of contracted Gaussian-type orbitals, NCGTO, per monomer will be documented in Figures 2 and 3, respectively, together with the corresponding MP4 values. It is pleasing that both DFT methods excellently mimic the MP4 curves, though as expected, the DFT total energy itself is less sensitive to NCGTO

Figure 4. Total energy of the formamide molecule obtained using different methods with the five atomic basis sets defined in Table 1.

(Figure 4). On the other hand, it is interesting to observe that the extrapolation of the MP4 energy (which is still far from being convergent) to the hypothetical infinite basis set limit would provide -169.845 hartrees, in very good agreement with the corresponding values of -169.833 and -169.872 hartrees, obtained for BLYP and B3LYP extrapolations, respectively. With subtraction of the previously obtained HF limit of -169.016 hartrees,10 the MP4 estimate of the correlation energy of the formamide molecule is Ecorr. ) -0.829 hartree. Measured on this value, we can see from Table 2 that the VWN potential with Ecorr. ) -2.052 hartrees seriously overestimates correlation (in both of the SVWN and BVWN combinations), while LYP gives with Ecorr. ) -0.824 hartree a surprisingly good account of it. Also the formamide dimer was the subject of a sustained interest over the past 2 decades, and it provides, therefore, a rich and reliable further basis for the evaluation of our theoretical models concerning their performance in describing the relevant intermolecular interactions. Specifically, we would like to see how accurately the basis sets used above for the monomer, can predict the geometrical and energetic parameters of hydrogen

3954 J. Phys. Chem., Vol. 100, No. 10, 1996

Suhai TABLE 4: H-Bond Lengths, RNsH‚‚‚O, and Binding Energies per H-Bond, ∆Er(CPC), According to Equation 10, of Formamide Dimers Computed Using the BLYP and B3LYP Potentials and the Atomic Basis Sets Defined in Table 1a zig-zag dimer

Figure 5. Structure of the zig-zag and cyclic dimers of formamide.

TABLE 3: H-Bond Lengths, RNsH‚‚‚O, and Binding Energies per H-Bond, ∆Er(CPC), According to Equation 10, in Formamide Dimers Computed Using Different DFT Potentials and the DZ(d,p) Basis seta zig-zag dimer 4

cyclic dimer

method

RNsH‚‚‚O

∆E (CPC)

RNsH‚‚‚O

∆Er(CPC)

HFS HFB SVWN SLYP BVWN BLYP B3LYP HFLYP HF MP2 MP4

2.8663 3.2483 2.7830 2.7384 3.1057 3.0051 2.9819 2.9237 3.1002 2.9985 3.0107

-8.12 -3.72 -9.49 -11.14 -4.04 -5.51 -5.62 -6.96 -4.91 -5.46 -5.51

2.7224 3.0616 2.6469 2.6117 2.9438 2.8508 2.8443 2.8287 2.9923 2.8802 2.8949

-10.82 -4.79 -12.37 -14.11 -5.46 -7.15 -7.26 -8.38 -6.20 -6.87 -7.27

a For comparison, the HF and MBPT results from ref 10 are also included (distances in angstroms; energies in kilocalories per mole).

bonding using various DFT schemes and how far the BSSE plays a role in these calculations. Two formamide molecules may have, of course, a number of relative configurations. For our later investigations, two configurations will be of special interest, the zig-zag and cyclic dimers, respectively, as depicted schematically in Figure 5. The structure, energetics, and force fields of the cyclic dimer have recently investigated by Florian and Johnson54 using MP2 and DFT methods, respectively. To get a first impression of the capability of the proposed exchangecorrelation potentials in predicting the lengths of the H-bonds in these systems, we used first the medium-sized DZ(d,p) basis set to determine the various dimer geometries. Though Table 3 containes only the NsH‚‚‚O distances for the two types of formamide dimers, all structures have been fully optimized for each basis set. The hydrogen bond will be taken as linear for the first model, but it will be bent for the cyclic one. Some qualitative tendencies will be immediately apparent by analyzing these results. At first, for a given theoretical level, the H-bond lengths within the cyclic dimer are always shorter than those obtained for the zig-zag structure. This is in contradiction to the previously cited experimental values of the NsH‚‚‚O distances in the formamide crystal and should be traced back to fundamental differences in the corresponding crystal field effects. We shall return to this point in the next section. At second, taking the MP2 and MP4 distances as a measure, we can see that the methods using local exchange (HFS, SVWN, and SLYP) seriously underestimate RNsH‚‚‚O. The opposite is true for the methods incorporating the gradient-corrected Becke exchange alone (HFB) or together with a local correlation functional (BVWN). Only the functional combinations BLYP and B3LYP seem to be able to reasonably reproduce the MBPT values. They deserve, therefore, a more

cyclic dimer

method/basis set

RNsH‚‚‚O

∆Er(CPC)

RNsH‚‚‚O

∆Er(CPC)

BLYP/DZ BLYP/DZ(d,p) BLYP/DZ(2d,2p) BLYP/TZ(2d,2p) BLYP/TZ(2df,2pd)

2.9682 3.0051 2.9762 3.0402 3.0450

-6.63 -5.51 -5.40 -4.45 -4.39

2.8420 2.8508 2.8272 2.9028 2.9094

-8.46 -7.15 -7.01 -5.58 -5.48

B3LYP/DZ B3LYP/DZ(d,p) B3LYP/DZ(2d,2p) B3LYP/TZ(2d,2p) B3LYP/TZ(2df,2pd)

2.9358 2.9819 2.9556 3.0084 3.0104

-7.71 -5.75 -5.61 -4.94 -4.88

2.7947 2.8443 2.8210 2.8878 2.8858

-8.78 -7.21 -7.06 -5.95 -5.87

HF/DZ HF/DZ(d,p) HF/DZ(2d,2p) HF/TZ(2d,2p) HF/TZ(2df,2pd)

3.0093 3.1002 3.1015 3.1464 3.1261

-5.89 -4.91 -4.82 -4.60 -4.56

2.8904 2.9923 2.9892 3.0054 3.0115

-7.18 -6.20 -5.71 -5.36 -5.21

MP2/DZ MP2/DZ(d,p) MP2/DZ(2d,2p) MP2/TZ(2d,2p) MP2/TZ(2df,2pd

2.9918 2.9985 2.9792 2.9856 3.0012

-6.46 -5.46 -5.18 -5.12 -5.08

2.8902 2.8802 2.8492 2.8506 2.8588

-7.90 -6.87 -6.44 -6.38 -6.52

MP4/DZ MP4/DZ(d,p)

2.9963 3.0107

-6.69 -5.51

2.9158 2.8949

-8.08 -7.27

a For comparison, the corresponding HF and MBPT results from ref 10 are also included (distances in angstroms; energies in kilocalories per mole).

TABLE 5: Different Components of the Binding Energy (kcal/mol) per H-Bond in the Cyclic Formamide Dimer According to Equations 8-10 Computed Using the BLYP Potential and the Atomic Basis Sets Defined in Table 1 method/basis set

∆E(BSSE)

∆E(CPC)

∆Er(CPC)

BLYP/DZ BLYP/DZ(d,p) BLYP/DZ(2d,2p) BLYP/TZ(2d,2p) BLYP/TZ(2df,2pd)

-9.98 -8.43 -8.67 -7.19 -6.97

-9.37 -7.93 -7.89 -6.52 -6.39

-8.46 -7.15 -7.01 -5.58 -5.48

detailed basis set study. Its results are presented in Table 4 together with the corresponding HF and MBPT values obtained earlier.10 In view of our numbers of highest quality, the MP2/TZ(2df,2pd) values in Table 4, both DFT methods perform excellently, the B3LYP results for RNsH‚‚‚O being a little bit closer to the MP2 ones for both structures. Basis set effects seem to play a larger role for DFT than for MP2, especially the DZ to TZ transition, which expands the H-bonds by 0.05-0.06 Å for the zig-zag structure and by 0.07-0.08 Å for the cyclic one, while such effects seem to be absent in MP2. The DFT binding energies closely follow the MP2 trends and also predict the increased binding in the cyclic structure. We may conclude from these results that at least a TZ(2d,2p)-type basis set will have to be used for the infinite systems to be able to reliably predict structural details at the resolution of a few hundredths of an angstrom. In addition, to be able to assess the role of the various correction terms in the DFT interaction energies (eqs 8-10), we collected them for different basis sets in the case of the cyclic dimer (Table 5). Finally, the formamide dimers provide us an adequate model system to investigate the role of the local and nonlocal (gradient correction) parts of the DFT correlation potentials in describing intermolecular hydrogen bonds. Taking the HFLYP functional

Density Functional Theory of Molecular Solids

J. Phys. Chem., Vol. 100, No. 10, 1996 3955

Figure 7. Effect of local, l-LYP, and nonlocal (gradient correction) parts, nl-LYP, of the DFT correlation potential in determining the binding energy in the cyclic (cyc) formamide dimer.

Figure 6. Effect of local, l-LYP, and nonlocal (gradient correction) parts, nl-LYP, of the DFT correlation potential in determining the length of the hydrogen bond in the zig-zag (z-z) and cyclic (cyc) structures of the formamide dimer.

as a typical example, Figure 6 demonstrates the effect of gradually switching on first the local part (l-LYP) and then the nonlocal part (nl-LYP), respectively, of the LYP functional. The contraction of the RN-H‚‚‚O distance nearly linearly depends on the percentage of both parts of the correlation term, and their separate contribution to ∆R is approximately equal. The same linear dependence can be observed for the binding energy of the dimer (∆E) but, interestingly, the total contribution from the gradient correction is twice as large as the local one (Figure 7). Similar conclusions hold for other potential combinations not shown here. We would like to also mention that our results concerning both the structure and the energetics of the cyclic formamide dimer are in excellent quantitative agreement with the corresponding ones obtained by Hrouda et al.53 and by Florian and Johnson54 apart from small differences due to basis set effects. IV. Infinite One-Dimensional Chain of Formamide Molecules Concerning the cooperative mechanisms of bonding in the formamide crystal, we first present results for the onedimensional model system, a chain of formamide molecules built as an infinite array according to the relative geometry of the zig-zag dimers shown in Figure 5. This chain is in fact a onedimensional infinite segment of the formamide crystal cut out along the [010] direction in Figure 1. Table 6 demonstrates that, independently of the basis set applied, cooperative interactions within the chain substantially reduce the lengths of the hydrogen bonds for both DFT potentials. For comparison, the corresponding HF and MP2 results10 are included again. A similar contraction of the hydrogen bonded network was observed in clusters of formamide molecules61-63 and in acetic acid oligomers.64 Its development may be demonstrated by MO results obtained for oligomers using the BLYP/DZ(d,p) method. The (F)n clusters (n ) 2, 3, 4, 5) were built according to the geometry of the zig-zag chain. From these results, Figure 8 shows how the NsH‚‚‚O distances in the innermost H-bonds of the clusters will be reduced with increasing n. Parallel to this effect, other structural parameters significantly change as

well. Two aspects of these model studies deserve to be emphasized. Firstly, the convergence of the geometrical and energetic parameters is rather slow. On the average, even the five-membered clusters recover only 60-70% of the changes observed between the dimer and the infinite chains. Secondly, such calculations become very expensive with increasing cluster size (and basis set). The basic difference between these clustertype model calculations and the crystal orbital studies is that while in the case of the periodic one-dimensional polymers the size of the computational problem is proportional to only one formamide molecule (from which the whole system can be generated by a screw symmetry operation consisting of a translation along the polymer axis and of a rotation around it by 180°), for the clusters the size of the problem steadily increases with n. Though calculations with systematically increasing cluster size will be able to predict certain properties for the bulk system, at higher theoretical levels their application becomes computationally very demanding. Furthermore, beyond this economical point of view, there are physical properties for which cluster boundary effects (e.g., accumulating surface charges at the ends of the cluster) exhibit an inherent problem. From the quantitative point of view, in going from the dimers (Table 4) to the polymers (Table 6) within the DFT theory, the length of the NsH‚‚‚C bond is shrinking for all basis sets by 0.16-0.19 Å. The same tendency with 0.10-0.12 Å was obtained at the MP2 level.10 Furthermore, we may observe an additional cohesion energy per hydrogen bond of about 6070% due to cooperative effects for both DFT methods. On the other hand, the RNsH‚‚‚O values provided by the DFT/BLYP method turn out to be closer to the corresponding MP2 predictions for larger basis sets than the B3LYP results. Since the optimization of the two-dimensional structure is very expensive due to its increased unit cell, we will use the former potential in the subsequent calculations. V. Infinite Two-Dimensional Network in the Formamide Crystal According to the X-ray diffraction investigations of Ladell and Post,1 formamide crystallizes in a monoclinic unit cell with a ) 3.69 ( 0.01 Å, b ) 9.18 ( 0.025 Å, c ) 6.87 ( 0.02 Å, and β ) 98 ( 0.25°. The space group is P21/n with 4 molecules per unit cell. Through a two-dimensional network of hydrogen bonds, the HCONH2 units associate to sheets as schematically depicted in Figure 1. These sheets are parallel to the (101) plane of the crystal. The basic structural elements of these sheets

3956 J. Phys. Chem., Vol. 100, No. 10, 1996

Suhai

TABLE 6: H-Bond Lengths, RNsH‚‚‚O, and Binding Energies per H-Bond, ∆Er(CPC), in Infinite One-Dimensional Zig-Zag Chains of Formamide Computed Using Two DFT Potentials with Atomic Basis Sets Defined in Table 1a BLYP

B3LYP

basis set

RHsH‚‚‚O

∆Ee(CPC)

DZ DZ(d,p) DZ(2d,2p) TZ(2d,2p) TZ(2df,2pd)

2.7711 2.7981 2.7821 2.8462 2.8516

-10.54 -8.96 -8.77 -7.78 -7.66

a

HF

RNsH‚‚‚O

∆Er(CPC)

2.7116 2.7796 2.7605 2.8146 2.8209

-11.76 -9.39 -9.24 -8.55 -8.42

MP2

RNsH‚‚‚O

∆Er(CPC)

RNsH‚‚‚O

∆Er(CPC)

2.8714 2.9558 2.9504 2.9916 2.9644

-10.10 -7.77 -7.72 -7.28 -7.44

2.8664 2.8756 2.8678 2.8970 2.8784

-10.80 -8.65 -8.50 -8.02 -8.14

For comparison, the corresponding HF and MP2 results from ref 10 are also included (distance in angstroms; energies in kilocalories per mole).

TABLE 7: NsH‚‚‚O Distances in the Two-Dimensional Formamide Crystal (Figure 1) Computed Using the DFT/ BLYP Method with Three Different Basis Sets Defined in Table 1a H-bonds along the [010] axis

H-bonds along the [1h01] axis

basis set

BLYP

HF

MP2

BLYP

HF

MP2

DZ DZ(d,p) TZ(2d,2p)

2.7580 2.7864 2.8357

2.8594 2.9411

2.8466 2.8579

2.8863 2.8992 2.9431

2.9396 3.0377

2.9388 2.9230

a For comparison, the corresponding HF and MP2 results (Å) from ref 10 are also included.

TABLE 8: Binding Energy per Formamide Molecule in the Two-Dimensional Infinite Formamide Crystal (Figure 1) Computed Using the DFT/BLYP Method with Three Different Basis Sets Defined in Table 1a Figure 8. Variation of the hydrogen bond length, RNsH‚‚‚O, in the oligomers and in the infinite polymer (after the break) of formamide: Zig-zag geometry, DZ(d,p) atomic basis set, DFT/BLYP theory.

are coplanar bimolecular associations (cyclic dimers) whose tilts relative to one another leads to a puckered structure of the sheets. Within each sheet hydrogen bonds of two types exist: one type, 2.935 Å long in ref 1, links monomers together to form the cyclic dimers; the other type, 2.880 Å long, links these dimers together to form the sheets. The structure may be alternatively described as consisting of one-dimensional chains of formamide molecules cross-linked by perpendicular hydrogen bonds to form twodimensional sheets. Adjacent sheets are held together by van der Waals forces. The closest approach between atoms in different sheets is 3.39 Å, i.e., between a carbon atom of one layer and a nitrogen atom of another. The structure of the infinite two-dimensional model of formamide is reminiscent of the parallel β-pleated sheets in proteins whose electronic structure was previously investigated at the HF and MP2 levels, respectively.65,66 Several technical details of crystal orbital calculations for such systems were discussed in ref 66. A common feature of the protein and formamide studies is the extensive use of Wannier functions (WFs) to perform the lattice sums in computing different perturbation theoretical terms. The technique to produce symmetry adapted real WFs for hydrogen bonded molecular crystals and to ensure the electrostatically balanced, symmetry preserving truncation of the infinite lattice sums is discussed in more detail in ref 66, together with the numerical problems related to the Brillouin zone integrations and the treatment of distant electrostatic interactions by multipole expansions. Besides the structural and cohesive properties, our method also provides the electronic band structure of the crystal, a useful quantity for investigating further properties. The study of the above protein model,66 along with numerous other polymer crystals2-6 showed, however, that the single particle band structures obtained by solving equation 3 will have to be further corrected for correlation effects to give meaningful band gaps, excitation energies, etc. In the case of the parallel β-pleated sheets, for example, the HF band gaps were reduced by more

Basis set

BLYP

HF

MP2

DZ DZ(d,p) TZ(2d,2p)

-14.55 -11.32 -10.26

-14.23 -9.66

-15.29 -10.63

a For comparison, the corresponding HF and MP2 results (kcal/mol) from ref 10 are also included.

than 50% due to correlation. Therefore, in the present analysis we will restrict ourselves to the optimized structures and total (binding) energies and leave the presentation of detailed electronic band structures to a future publication. Table 7 presents the optimized hydrogen bond lengths along the [010] and [1h01] directions, respectively. In the case of the two-dimensional formamide models with four molecules in the unit cell, computational resources limited our structure optimizations for the three basis sets up to TZ(2d,2p). At the MP2 level, the same limitations allowed us to include at most the DZ(d,p) basis set.10 The comparative calculations for the onedimensional chain cited in Table 6 indicate, however, that while the MP2/DZ(d,p) results are very close to the MP2 values obtained for the larger basis set, this does not seem to hold for DFT, for which only the TZ(2d,2p) basis set provides satisfactory saturation. The [010] direction corresponds to the hydrogen bonded zig-zag chains in the crystal (Figure 1), and by comparing the first part of Table 7 with the one-dimensional values in Table 6, we can see that the lateral interactions within the two-dimensional sheet slightly strengthen the hydrogen bonds along this chain and decrease RNsH‚‚‚O by 0.011 Å for BLYP/TZ(2d,2p). Similar effects were observed for both at HF and MP2, respectively.10 On the other hand, the H-bonds in the perpendicular direction, along the [1h01] axis, will be weakened as compared with the isolated cyclic dimer: the BLYP/TZ(2d,2p) value of RNsH‚‚‚O ) 2.9028 Å (Table 4) will be elongated to 2.9431 Å in the infinite system (Table 7). Though our predicted DFT H-bond length of 2.8357 Å along the [010] direction is somewhat shorter than the best MP2 value of 2.8579 Å,10 it can still be reasonably compared with the corresponding experimental values, 2.880 Å (X-ray at 223 K1),

Density Functional Theory of Molecular Solids

J. Phys. Chem., Vol. 100, No. 10, 1996 3957

TABLE 9: Comparison of the Structural Parameters of the Formamide Molecule and of the Crystal, Respectively, Obtained with Different Experimental and Theoretical Methodsi method monomer EDb MWc MO/BLYP/TZ(2d,2p)d MO/MP4/TZ(2df,2pd)e crystal X-ray at 223 Kf X-ray at 90 Kg Neutron diffr. at 7 Kh CO/BLYP/TZ(2d,2p)d CO/MP2/DZ(d,p)e

RC,O

RC,N

RN,Ha

1.211 1.219 1.2226 1.2179

1.367 1.352 1.3698 1.3600

1.021 1.002 1.0140 1.0014

1.255 1.239 1.239(7) 1.2553 1.2356

1.300 1.326 1.294(5) 1.3548 1.3312

1.010 1.048(6) 1.0353 1.0216

b

c

9.180 9.041 8.9512(5) 9.0236 8.9942

6.870 6.994 6.9741(4) 6.9278 6.9580

a Proton involved in the hydrogen bond. b Reference 59. c Reference 60. d Present paper. e Reference 10. f Reference 1. g Reference 67. h Reference 68. i The molecular and crystal orbital (MO and CO) results refer to calculations performed for the isolated monomer and for the infinite twodimensional crystal, respectively. b and c are the lattice parameters according to Figure 1. All quantities in angstrom.

2.885 Å (X-ray at 90 K 67), and 2.853 Å (neutron diffraction at 7 K 68), respectively. In the perpendicular direction, along the (1h01) axis, the predicted value of 2.9431 Å favorably compares with 2.935, 2.948, and 2.958 Å, respectively, provided by the above three experiments (the corresponding MP2 value is 2.9230 Å 10). Table 8 contains the DFT/BLYP binding energies per formula unit obtained in the infinite crystal. They have been corrected again by the CPC formula using the basis sets of the first lattice neighbors in the strict sense.5 We may observe, comparing these values with the corresponding values obtained for the segregated zig-zag chains in Table 6, that interchain interactions further stabilize the crystal by about 30%. The formation of the hydrogen bonds influences the internal geometry of the constituents in the crystalline environment as well. In Table 9 we summarize the values of the three most important bond lengths of formamide, observed experimentally and predicted theroetically by different molecular orbital (MO) and crystal orbital (CO) methods. As discussed earlier, for the gas phase the MO/BLYP/TZ(2d,2p) values for the bond distances are in good agreement with the MP4 results and with measurements. Comparing, on the other hand, the single molecule electron diffraction and microwave values with those measured in the crystals, we can see that hydrogen bonding typically elongates the CdO and NsH bonds, and slightly reduces the CsN one. From the two groups of measurements, we may approximately infer ∆RC,O ) +0.028 to +0.036 Å, ∆RC,N ) -0.026 to -0.073 Å, and ∆RN,H ) +0.010 to +0.046 Å, respectively. The differences between the MO and CO results clearly show the same tendencies, providing ∆RC,O ) +0.0327 Å, ∆RC,N ) -0.0150 Å, and ∆RN,H ) +0.0213 Å, respectively, for the BLYP/TZ(2d,2p) procedure. The corresponding MP2/DZ(d,p) values are +0.0075, -0.0354, and +0.0161 Å, respectively. The theoretically predicted lattice vector components of b ) 9.024 Å and c ) 6.928 Å favor the values measured at low temperatures by neutron diffraction.68 Acknowledgment. The author is indebted to Gerd Raether and Bertram Scholz for their contribution in programming the DFT moduls of the POLYGAUSS system, to Tobias Reber for his continuous system advice, to Istva´n Mayer for useful suggestions concerning the BSSE problem, and to Bela Paizs and Eike Bergner for technical assistance in computing the BSSE correction terms. The use of the X - R code of John W. Mintmire in the early phase of this work and the generous support with computer time on the CONVEX and CRAY supercomputers and on IBM/SP2 and SGI Power Challenge clusters of DKFZ and of the University of Stuttgart are also gratefully acknowledged. This research has been partly supported by the Commission of the European Union (Grant No.

CT-93-006) and by the German Federal Ministry of Research and Technology (Grant No. 01 1B 303A). References and Notes (1) Ladell, J.; Post, B. Acta Crystallogr. 1954, 7, 559. (2) Suhai, S. Phys. ReV. B 1983, 27, 3506. (3) Suhai, S. Phys. ReV. B 1984, 29, 4570. (4) Suhai, S. J. Chem. Phys. 1986, 84, 5071. (5) Suhai, S. Collective Excitations in Biological Macromolecules Photoelectron and Exciton Spectra of Polyene, Polypeptide, and Polynucleotides. In Molecules in Physics, Chemistry, and Biology; Maruani, J., Ed.; Kluwer: Dordrecht, The Netherlands, 1989; Vol. IV, pp 133-194. (6) Suhai, S. Int. J. Quantum Chem. 1992, 42, 193. (7) Suhai, S. Int. J. Quantum Chem. Quantum Chem. Symp. 1993, 27, 131. (8) Suhai, S. J. Chem. Phys. 1994, 101, 9766. (9) Suhai, S. Phys. ReV. B 1994, 50, 14791. (10) Suhai, S. J. Chem. Phys. 1995, 103, 7030. (11) Hohenberg, P.; Kohn, W. Phys. ReV. B 1964, 136, 864. (12) Kohn, W.; Sham, L. J. Phys. ReV. A 1965, 140, 1133. (13) Senatore, G.; March, N. H. ReV. Mod. Phys. 1994, 66, 445. (14) Kryachko, E. S.; Ludena, E. V. Energy Density Functional Theory of Many Electron Systems; Kluwer: Dordrecht, The Netherlands, 1990. (15) Trickey, S. B., Ed. AdVances in Quantum Chemistry, Vol. 21; Academic Press: New York, 1990. (16) Parr, G. R.; Yang, W. Density Functional Theory of Atoms and Molecules; Oxford University Press: New York, 1989. (17) Labanowski, J. K., Andzehn, J., Eds.; Density Functional Methods in Chemistry; Springer-Verlag: New York, 1991. (18) Ziegler, T. Chem. ReV. 1991, 91, 651. (19) Mintmire, J. W.; White, C. T. Phys. ReV. Lett. 1983, 50, 101. (20) te Velde, G.; Berendsen, E. F. Phys. ReV. B 1991, 44, 7888. (21) von Boehm, J.; Kuivalainen, P.; Calais, J.-L. Phys. ReV. B. 1987, 35, 8177. (22) Springborn, M.; Andersen, O. K. J. Chem. Phys. 1987, 87, 7125. (23) Suhai, S. J. Phys. Chem. 1995, 99, 1172. (24) Suhai, S. Phys. ReV. B 1995, 51, 16553. (25) Suhai, S. Phys. ReV. B 1995, 52, 1674. (26) Del Re, G.; Ladik, J.; Biczo´, G. Phys. ReV. 1967, 155, 967. (27) Andre´, J. M.; Gouvameur, L.; Leroy, G. Int. J. Quantum Chem. 1967, 1, 427. (28) Gill, P. M. W.; Johnson, B. G.; Pople, J. A. Chem. Phys. Lett. 1994, 217, 65. (29) Slater, J. C. Quantum Theory of Molecules and Solids; McGrawHill: New York, 1974; Vol. 4. (30) Colle, R.; Salvetti, O. Theor. Chim. Acta 1975, 37, 329. (31) Vosko, S. H.; Wilk, L.; Nusair, M. Can. J. Phys. 1980, 58, 1200. (32) Perdew, J. P.; Zunger, A. Phys. ReV. B 1981, 23, 5048. (33) Perdew, J. P. Phys. ReV. B 1983, 33, 8822. (34) Perdew, J. P.; Wang, Y. Phys. ReV. B 1986, 33, 8800. (35) Becke, A. D. Phys. ReV. A 1988, 38, 3098. (36) Lee, C.; Yang, W.; Parr, R. G. Phys. ReV. B 1988, 37, 785. (37) Mielich, B.; Savin, A.; Stoll, H.; Preuss, H. Chem. Phys. Lett. 1989, 157, 200. (38) Pople, J. A.; Binkley, J. S.; Seeger. Int. J. Quantum Chem. Symp. 1976, 10, 1. (39) Krishnan, R.; Pople, J. A. Int. J. Quantum Chem. 1978, 14, 19. (40) Krishnan, R.; Pople, J. A.; Replogle, E. S.; Head-Gordon, M. J. Phys. Chem. 1990, 94, 5579. (41) Frisch, M. J.; Trucks, G. W.; Head-Gordon, M.; Gill, P. M. W.; Gong, M. W.; Foresman, J. B.; Johnson, B. G.; Schlegel, H. B.; Robb, M.

3958 J. Phys. Chem., Vol. 100, No. 10, 1996 A.; Replogle, E. S.; Gomperts; Andres, J. L.; Raghavachari, K.; Binkley, J. S.; Gonzalez, C.; Martin, R. L.; Fox, D. J.; Defrees, D. J.; Baker, J.; Stewart, J. J. P.; Pople, J. A. Gaussian 92; Gaussian: Pittsburgh, PA, 1992. (42) Raether, G.; Suhai, S. Chem. Phys., to be published. (43) Becke, A. D. J. Chem. Phys. 1993, 98, 5648. (44) Johnson, B. G.; Gill, P. M. W.; Pople, J. A. J. Chem. Phys. 1993, 98, 5612. (45) Jansen, H. B.; Rose, P. Chem. Phys. Lett. 1969, 3, 140. (46) Boys, S. F.; Bernardi, F. Mol. Phys. 1970, 19, 553. (47) Mayer, I.; Surja´n, P. R. Chem. Phys. Lett. 1992, 191, 497. (48) Kwiatkowski, J. S.; Leszczynski, J. J. Mol. Struct. 1992, 270, 67. (49) Kwiatlowski, J. S.; Lesczczynski, J. J. Mol. Struct. 1993, 297, 277. (50) Burton, N. A.; Chiu, S. S.-L.; Davidson, M. M.; Green, D. V. S.; Hilliar, I. H.; McDouall, J. J. W.; Vincent, M. A. J. Chem. Soc., Faraday Trans. 1993, 89, 2631. (51) Lisini, A.; Keane, M. P.; Lunell, S.; Correia, N.; Naves de Brito, A.; Svensson, S. Chem. Phys. 1993, 169, 379. (52) Florian, J.; Johnson, B. G. J. Phys. Chem. 1994, 98, 3681. (53) Hrouda, V.; Florian, J.; Polasek, M.; Hobza, P. J. Phys. Chem. 1994, 98, 4742. (54) Florian, J.; Johnson, B. G. J. Phys. Chem. 1995, 99, 5899.

Suhai (55) Dunning, T. H.; Hay, P. J. Modern Theoretical Chemistry; Plenum: New York, 1976, p 1. (56) Dunning, T. H. J. Chem. Phys. 1971, 55, 716. (57) Dunning, T. H. J. Chem. Phys. 1989, 90, 1007. (58) Wiberg, K. B.; Hadad, M.; LePage, T. J.; Breneman, C. M.; Frisch, M. J. J. J. Phys. Chem. 1992, 96, 671. (59) Kitano, M.; Kuchitsu, K. Bull. Chem. Soc. Jpn. 1974, 47, 67. (60) Hirota, E.; Sugisaki, R.; Nielsen, C.; Soerensen, G. J. Mol. Spectrosc. 1974, 49, 251. (61) Hinton, J. F.; Harpool, R. D. J. Am. Chem. Soc. 1977, 99, 349. (62) Mehler, E. L. J. Am. Chem. Soc. 1980, 102, 4051. (63) Shivaglal, M. C.; Singh, S. Int. J. Quantum Chem. 1992, 44, 679. (64) Turi, M.; Dannenberg, J. J. J. Am. Chem. Soc. 1994, 116, 8714. (65) Seel, M. Solid State Commun. 1990, 45, 1563. (66) Suhai, S. Int. J. Quantum Chem. 1991, 40, 559. (67) Stevens, E. D. Acta Crystallogr. 1978, B34, 544. (68) Torrie, B. H.; O’Donovan, C. O.; Powell, B. M. Mol. Phys. 1994, 82, 643.

JP9526399