Development of a Cyclic Periodic Wave Function Approach for the Study of Infinitely Periodic Solid-State Systems.

The ab initio cyclic periodic wave function (CPWF) approach is developed for the treatment of infinitely periodic systems. Using the full infinite Hamiltonian operator, as well as symmetrically identical basis set wave functions that preserve the translational symmetry of the electron density of the system, this approach can be applied at the Hartree-Fock level, or correlation can be directly included by the usual modes. In this approach, all many-body interactions are included, and no edge effects occur. Initial test calculations of the CPWF method at the ab initio Hartree-Fock level are performed on the chains of hydrogen fluoride molecules.

Introduction

In the past, many models and approaches have been developed for the determination of the bulk chemical properties in solid-state systems, including some approaches which are used to study infinitely periodic molecular solids.

The simplest and most popular models for studying infinitely periodic molecular solids are cluster models[1−5] and pair- potential, or N-body potential, methods.[6] The most common approaches to their studies are through the use of Bloch orbitals[7] and density function theory (DFT).[8,9]

Since the 1970s, cluster models have been applied in many kinds of solids.[1−5,10] The assumption in cluster models is that the electronic structure of the infinite crystal can be reasonably simulated by performing high-level calculations on a system consisting of a few representative molecules that are geometrically selected to imitate the overall crystal structure, as shown in Figure . By systematically increasing the number of molecules in the clusters, it is hoped that the properties of the clusters will smoothly converge to those of the crystal, if a large enough cluster is chosen. In a review of cluster models, Fink et al.[11] pointed out that cluster termination (or edge) effects are the major deficiency of the approach, when a finite cluster is used to simulate an extended or infinite system. Adachi[12] concluded that systematically increasing the cluster size in the solid will eventually yield a large enough cluster to provide an adequate model for the bulk crystal. However, this approach requires great computational effort to obtain a reasonable representation of the system studied and reduce the edge effects. Zunger proposed the cyclic cluster model (CCM) in 1970s,[13−15] which places the atoms themselves at cyclic periodic positions in space. Jug et al. successfully implemented the CCM at different levels of calculations to various crystalline systems.[16−20] Janetzko et al.[21] applied the CCM to different covalent periodic systems, treated by the Kohn–Sham auxiliary density function theory, which successfully retained the symmetry of the corresponding periodic systems.

Figure 1

Three different cluster models labeled with solid line―, dashed line---, and dotted line.... respectively, for an imaginary periodic system.

In pair-potential or N-body potential methods, an important assumption[6] is that the intermolecular energies of N molecules can be sufficiently well approximated by the sum of interactions between N(N – 1)/2 individual pairs. This approach works very well when the interaction energies of three-body or higher-body interactions are negligible compared to the two-body interactions.[22−24] Raynor[4] showed that large oscillations in error occurred in studies of molecular hydrogen solid at high pressure as higher many-body terms were added. In addition, studies on solid molecular hydrogen showed that anisotropies become more important either at short intermolecular distance,[25] or at high pressure.[22] Although those anisotropies make the pair-potential methods more difficult to use, pair potentials are still being applied to model some solid-state systems.[26−30] More recently, Bučko et al.[31,32] successfully used a pair interactions model with a DFT approach to treat molecular crystals, which can be applied at any level of approximations.

For the study of weak interactions in large molecular systems, it is necessary to simultaneously use quite large basis sets and reasonably high-order correlation corrections,[33−35] which is still very computationally demanding. To solve some of the problems involved with these systems, we have developed a very promising new ab initio method, called the cyclic periodic wave function (CPWF) method, which uses the complete infinite Hamiltonian for the infinite crystal and preserves the translational symmetry of the system in the wave function. In this paper, we will show the results of some initial test calculations of our method at the ab initio Hartree–Fock level for one-dimensional (1-D) systems. We will test and apply our CPWF approach at the AM1 and PM3 semiempirical levels of approximation to a series of two-dimensional (2-D) and three-dimensional (3-D) systems to demonstrate the viability and power of our approach, and the results will be presented in two subsequent papers.

Computational Methods

The novel ab initio cyclic periodic wave function (CPWF) approach uses the full infinite Hamiltonian operator for the crystal. First, the infinitely periodic system is divided into chemically equivalent repeat units, as demonstrated by Figure b. Figure a,b shows two different methods for choosing repeat units in an infinitely periodic line of molecules. Figure a places the cell boundaries in positions that truncate molecules, while Figure b places them between molecules. The cell boundaries in Figure b are such that the coefficients for the molecules at odd positions remain identical to those at even positions, no matter how many cell interactions are included in the calculation. The full infinite Hamiltonian operator for the system is then be written aswhere I and J denote individual repeat units, i and j denote electron labels, and A and B denote the atomic centers within repeat units I and J.

Figure 2

Two alternative unit cell definitions for a Bloch orbital model of an infinite zig-zag line of molecules. (a) Cell boundaries truncate molecules; (b) cell boundaries occur between molecules.

To solve the Schrödinger equation for the infinitely repeating system, we will assume that each repeat unit in the perfect crystal is electronically identical to all others, differing only perhaps by one or more symmetry transformations. Since each repeat unit therefore has a chemical environment that is identical to that of every other repeat unit in the system, the electron density distribution near repeat unit I must be kept identical to that near repeat unit J. Thus, we will begin by defining the ground-state Slater determinant of repeat unit molecular orbitals to have the following formwhere the molecular orbitals (MOs) associated with repeat unit I must have the identical mathematical form as those associated with repeat unit J—i.e., the MOs on repeat unit I are simple translations, and perhaps also rotations or reflections, of those on repeat unit J. Thus, if we can determine the “best” form for the MOs on one repeat unit, we will have determined the best form for all repeat units in the infinite system. The assumption is that a single Slater determinant ground-state wave function can be readily relaxed by incorporating Møller–Plesset perturbation theory corrections to recover the correlation energy.[34,35]

To create infinitely periodic molecular orbitals which preserve the translational symmetry in the system, we follow a similar prescription to that used by us in studies of infinitely periodic molecular and ionic solids,[4,36−40] but where the Hartree–Fock–Roothaan coefficients between neighboring repeat units are directly evaluated.

The basis set used in CPWF is constructed from a set of cyclic periodic wave functions, instead of traditional Bloch orbitals. The principal difference stems from using the individual repeat units themselves in the complex exponential expansion, rather than the unit cell, as is shown in the equationNote that eq is dependent upon the assignment of position indices (J, J, J). Each repeat unit is first assigned a unique set of three independent position indices (J, J, J), from which the geographic position for the center of each unit can be determined. In eq , ϕ( is an atomic orbital located on the monomer at position (j, j, j). Note that this atomic orbital may be a linear combination of degenerate atomic orbitals, with the particular linear combination being determined by the orientation of the monomer at position (j, j, j), relative to its orientation at the (0,0,0) position. For example, suppose that the third and fourth atomic orbitals in the basis set on the monomer at position (0,0,0) are a p and a p orbital, respectively, and further suppose that the monomer in the (1,0,0) position is related to that at the (0,0,0) position by a clockwise rotation of 90° about the z-axis. Then, in order for the third and fourth atomic orbitals at (1,0,0) be symmetrically identical to the third and fourth ones at (0,0,0), they should be defined as followsNext, the repeat length N in each direction must be selected. Choosing N = 4 ensures that all of the interactions for repeat units within ±1 in that index direction are treated exactly, (these repeat units are referred to as nearest neighbors); choosing N = 6 ensures that all of the interactions for repeat units within ±2 in that index direction are treated exactly (these repeat units are referred to as next-nearest neighbors); and so on. Equation then ensures that the cyclic periodic basis set function, ψ, is N-fold periodic in the x-index, N-fold periodic in the y-index, and N-fold periodic in the z-index.

Next, the periodic basis set is symmetrically orthogonalized to give a new orthonormal basis setwhere j references different atomic wave functions for the atoms inside cell (J, J, J), and k references atomic wave functions in cell (K, K, K). The elements of the inverse half-root overlap matrix, S–1/2, in eq are generated in the usual way,[36,41] i.e., from the eigenvalues and eigenvectors of the overlap matrix formed from the cyclic periodic wave functions defined in eq . The limits for J, J, and J are defined relative to K, K, and K, as followswith similar relations for the y- and z-summations.

Finally, it is these functions that form the orthonormalized basis set that is used to determine the Hartree–Fock–Roothaan coefficients for the molecular orbitals. Thus, the ith MO for the repeat unit at position (I, I, I) are given bywhere K, K, and K are defined relative to I, I, and I, as followswith similar relations for the y- and z-summations.

Because the full infinite Hamiltonian is used and the infinite translational symmetry is preserved in both the basis functions and the Hamiltonian, the infinite Fock matrix becomes block diagonal. Thus, we need only diagonalize the block representing the range of (K, K, K) = (0,0,0) to (N – 1, N – 1, N – 1), since the coefficients derived from it will be identical to those needed for all other blocks. Furthermore, even this block is further diagonalized to produce N × N × N individual blocks containing complex elements, as a consequence of the fact that all of the matrix elements between periodic basis functions at different index positions are zero:whereNow, redefining the second sets of summation indices as followswhere Δ varies from −((1/2)N – 1) to + ((1/2)N – 1), with similar ranges for Δ and Δ. Note that any matrix element involving an interaction that is further away than N/2 – 1 in any position index will be 0, we get the following form for eq Since f( ≡ f(0,0,0)→(Δ, and because each term in the first set of summations reduces to the following,with similar relations for the y- and z-summations. We find that the only nonzero terms in eq will be those involving center J = center I, leading to the following formwhere the f(0,0,0)→(Δ elements are the Fock matrix elements for the interaction of atomic orbital i on the monomer located at (0,0,0) with atomic orbital j on the monomer located at (Δ,Δ,Δ). Thus, the Hartree–Fock problem becomes a matter of finding the eigenfunctions of N × N × N individual matrices, each of which is an M × M matrix, where M is the number of atomic orbitals per monomer. However, the more efficient eigenvalue–eigenfunction routines are designed for matrices with real components only. By taking appropriate linear combinations of the cyclic periodic wave functions in eq , these matrices can be transformed into fully real ones. The linear combinations needed are illustrated below in the equationwith similar transformations for K and K. The largest block that can occur after these transformations becomes one that involves matrices with dimensions of 8 M × 8 M. However, for most systems, further symmetry relations aid in reducing the blocks still further, and the largest blocks that have occurred in any 3-D system that we have studied thus far had dimensions of 4 M. Thus, a Hartree–Fock calculation on the full infinite solid-state system is no more computationally intense than doing a few calculations on the isolated monomer, with a basis set 4 times that used for the lone monomer.

Because the translational symmetry of the infinitely periodic crystal has been preserved in the chosen form for the wave functions, the final total electron density at each repeat unit will be identical to that at every other repeat unit. Furthermore, the total energy of the system becomes partitioned into identical components from each repeat unit. Thus, each repeat unit contributes an identical energy to the total infinite system. In the familiar terms of standard closed-shell Hartree–Fock calculations, the energy of the Ith repeat unit is given bywhere H and F are evaluated using the wave functions in eq , and density matrix elements P are as followsand where the D coefficients are those for the orthogonalized components to the Hartree–Fock coefficients,The summation of J = “nearest”, or “next nearest” in eq , means summations relative to I. For example, if I = 000, then J nearest would involve all terms with J = −1 or +1, J = −1 or +1, and J = −1 or +1, etc. The EMad term given in eq is the Madelung correction term that arises as a consequence of applying our method for correcting the imbalanced attractive and repulsive contributions that arise when truncated summations occur.[40] Since we include all intercharge terms in the infinite system into this term, it also acts to recover some of the electrostatic interactions that are missing as a result of using truncated sums. Our method converges more quickly than does the more commonly used Delhalle method.[42,43] Specifically, our method subtracts the average of the first terms in the multipole expansions of the Coulomb and nuclear attraction integrals from the integrals themselves, as they are added into the H and F matrices. The infinite sum of all such multipole expansions is then added to both the Fock matrix elements and the overall energy. These terms have the following formswhere qA is the Mulliken charge on atom A and whereOur process has two beneficial effects—it prevents the occurrence of unbalanced contributions from attractive or repulsive terms, and it accounts at least partially for any long-range charge–charge interactions that were not directly included in the calculation. Thus, even though finite values are used for the repeat lengths, N, N, and N, which limits the number of intermolecular interactions that are included exactly into the matrix elements, the effects from neighbors further out are at least partially recovered. Thus, our method is far less sensitive to the sizes chosen for the repeat lengths, N, N, and N, than would otherwise be expected.

Finally, using cyclic periodic wave functions, it is a straightforward process to determine basis set superposition errors (BSSE) using the Boys–Bernardi counterpoise method.[44] Since for each choice of the repeat length parameters (N, N, N), only certain specific monomers are explicitly included in the matrix elements, the calculation of the counterpoise correction to the energy of the isolated monomer in each case should contain exactly the same interactions in their matrices. Thus, for example, if N is 4, then all atomic orbitals on positions with I = −1, 0, and 1 would be used in the basis set expansion relative to an isolated monomer at the position (0,0,0) (i.e., with I = I = I = 0), but integrals involving |ΔI| > 1 would be set equal to 0. These integrals would then be partially recovered, in exactly the same way as for the infinite system, by adding in the Madelung correction terms for these interactions. Thus, the BSSE-corrected interaction energy per monomer would be determined as followswhere EI is the energy per repeat unit for the infinite system, as given in eq , ECPI is the counterpoise correction to the isolated repeat unit’s energy, calculated as described above, and ΔEdefI is the deformation energy, i.e., the energy necessary to deform the isolated repeat unit from its equilibrium geometry to its new geometry in the infinite system.

Finally, it should be noted that although the basis set wave functions are cyclically periodic, the positions of the repeat units are not. Thus, the physical positions of each repeat unit will be translationally periodic, i.e., the molecules (and their associated atomic orbitals) are placed at their expected physical positions within the crystal.

Results and Discussion

To gain some understanding of how sensitive the interaction energies in hydrogen-bonded systems will be to the choice of the repeat length parameter, we decided to investigate two models involving the strongest hydrogen bond—that between hydrogen and fluorine in hydrogen fluoride (HF), and perform ab initio calculations on chains of hydrogen fluoride molecules. The first model is an infinite linear chain of hydrogen fluoride molecules shown in Figure a, and the second one is the single infinite HF catemer chain shown in Figure b. Our goal is to examine how different choices of the repeat length parameter, N, affect the predicted hydrogen-bond energy in a one-dimensional system.

Figure 3

Models of HF molecules. (a) Simple linear chain model of HF molecules; (b) single infinite HF catemer chain model.

It has been clearly demonstrated that the accuracy of hydrogen-bond interactions is strongly influenced by the quality of the basis set used and by the degree of correlation included.[4,39,45,46] Slater-type orbitals are known to have better long-range properties than Gaussian-type orbitals, and the trends in energy as a function of repeat length N were more important to us than the attainment of high accuracy. There are some very successful software packages that have been developed using DFT so far, such as the GAUSSIAN program package,[47] the Vienna ab initio package (VASP)[48−52] for the study of molecular solid systems,[53,54] and the most recent CRYSTAL[55,56] and CRYSCOR[57,58] programs. In our case, a software package is available that applies our cyclic periodic wave function method using Slater orbitals, so we decided to employ a simple minimum basis set of Slater orbitals for our test calculations.

Effect of Repeat Index Variation on the Linear Chains of HF Molecules

The first set of calculations was performed on a very simple model: an infinite linear chain of hydrogen fluoride molecules, as illustrated in Figure a. Our purpose in choosing this model was twofold: (i) to investigate the effect of different choices for the repeat length parameter, N, on the calculated hydrogen-bond energies, and (ii) to compare the efficiency of our approach when compared to cluster-type calculations. For these initial calculations, we kept the molecular H–F bond distance fixed at 1.80 bohr (0.950 Å), which was the optimal distance found for the isolated HF molecule.

First, we investigated the effect on the H-bond energy of varying the F···F intermolecular distance, using three different values for N:N = 4, which includes nearest-neighbor interactions exactly; N = 6, which also includes next-nearest neighbors; and N = 8, which includes next-next-nearest neighbors, as well. In addition to varying the values of N, we also investigated the importance of the Madelung correction term EMad in eq . Note that all of the curves which included this term provide a significant lowering of the overall energy per monomer, as shown in Figure . Since this term requires very little computation time, all further discussion will be for those curves which explicitly include the calculation of EMad. As shown in Figure , the results for N = 6 and 8 are nearly indistinguishable from one another, predicting nearly identical values for the optimal F···F intermolecular distance and for the H-bond energy. Although the N = 4 calculations predict noticeably weaker H-bond energies than the other two, they recover most of the intermolecular H-bond energy (about 80%). The N = 4 results also predict a slightly longer F···F intermolecular distance, but it differs by only about 0.1 Å from the other two. Thus, even for a system with interaction energies nearly as strong as covalent ones, the energy converges very rapidly as the length parameter, N, is increased.

Figure 4

Effect of N, and the Madelung correction terms on the calculated H-bond energies, ΔHHBond vs the intermolecular distance RF···F. All dashed lines are without Madelung correction, and all solid lines are with Madelung correction. Lines with diamond symbols are for N = 4, lines with circle symbols are for N = 6, and lines with square symbols are for N = 8, respectively.

Next, we used the same infinite linear chain model to compare the trends in the repeat parameter, N, to the trends found using a truncated chain approach. For the cluster calculations, we calculated the H-bond energy in two different ways. The first corresponds to the difference between the average energy per monomer in the truncated chain and the energy of a lone monomerwhere E is the total energy of the N-mer and E1 is the energy of the isolated molecule. The second method views an N-mer as being composed of (N – 1)-mer plus a monomer and one hydrogen bond. The average contribution per hydrogen bond then becomeswith E and E1 are defined as above. Figure shows the effect of cluster sizes on the predicted hydrogen-bond energy convergence for the cluster models. Also shown are the results of our cyclic periodic wave function method, with N = 4, 6, and 8. Since N = 4 corresponds to exact inclusion of nearest-neighbor interactions, its value is plotted as corresponding to the amount of computational efforts on a two-membered cluster. Similarly, our results for N = 6 and 8 are plotted as corresponding to three- and four-membered clusters, respectively. What is clear from Figure is that our cyclic periodic wave function method converges far more rapidly to the infinite limit than do the cluster models. Also, it is clear that ΔE2 converges more rapidly than ΔE1 for the cluster models. However, regardless of the method used to determine it, the convergence of the cluster models is extremely slow when compared to the convergence of the cyclic periodic wave function method.

Figure 5

Comparison of the computational efforts for the H-bond energy convergence vs the cluster length. For the truncated chain approaches, the x-axis represents the number of molecules in the truncated chain. For the cyclic periodic approach, the x-axis represents the number of nearest neighbors explicitly included in the matrix elements. The line with square symbols shows the values of ΔE1 for the truncated chain models, eq , the line with diamond symbols shows the values of ΔE2 for the truncated chain models, eq , and the line with circle symbols shows the calculated H-bond energies from the cyclic periodic model. The dashed line indicates the infinite limit for all methods.

The importance of rapid convergence becomes even more apparent when we consider the effect on computation time. Two factors strongly affect the amount of time needed to perform a Hartree–Fock calculation on a system: the number of unique Coulomb-type integrals that must be evaluated and the size of the Fock matrix that must be diagonalized. As a measure of the first factor, we have plotted both the ΔE2 data and the cyclic periodic wave function data as a function of the number of unique Coulomb integrals that had to be calculated. These plots are shown in Figure . Here, we can clearly see the benefits of using a model that preserves the symmetry of the infinite system. For the cluster models, a relatively large increase in computation time produces very little improvement in convergence. This is because in the cluster models, each molecule in the cluster occupies a different chemical environment. This asymmetry greatly increases the complexity of the calculation.

Figure 6

Comparison of the calculated energy error as a function of the number of unique Coulomb integrals that must be evaluated. The line with diamond symbols is for values of ΔE2 for the truncated chain models; the line with circle symbols is for the calculated H-bond energies from the cyclic periodic model. The dashed line shows the infinite limit for all methods.

It is easy to compare these approaches with regard to the second effect, Fock matrix size. For a cluster model, the Fock matrix will have a dimension of NM × NM, where M is the number of basis functions used on each monomer and N is the number of monomers in the cluster. Since the time needed to diagonalize matrices is on the order of n3 for an n × n matrix, these times will rise very rapidly with cluster size. In our approach, we will need to perform several diagonalizations of smaller Fock matrices. For a one-dimensional chain of identical molecules (linear or zig-zag), we will need to diagonalize the following: two sets of M × M matrices and [(1/2)N – 1] sets of 2 M × 2 M matrices. Thus, increasing the repeat length does not affect the sizes of the matrices that are diagonalized—it only affects how many such matrices need to be diagonalized. As a result, the time dependence for the cluster models will be proportional to (NM)3, but it will be proportional to only 16[2N – 3]M3 for the cyclic periodic wave function approach. This means that performing a Hartree–Fock calculation on an infinite system is only moderately more time consuming than performing the calculation on an isolated molecule.

Results for Fully Optimized Catemers at SCF and Møllet–Plesset Levels

Since our method is wave function-based, we are able to directly determine correlation corrections using Møllet–Plesset perturbation theory. These corrections are demonstrated using our second model system, which replicates more closely the geometry of HF in the solid state. In the solid state, HF crystallizes in staggered chains,[59] as illustrated in Figure b. Using our cyclic periodic wave function approach, we studied the effect of varying the repeat length parameter, N, on the H-bond energy in a single infinite HF catemer chain, calculated with and without correlation. We once again employed a minimum basis set of Slater orbitals for these calculations. Since the experimental FHF angle is nearly linear at 176°,[59] we assumed that this angle remained exactly linear and optimized the H–F and F···F bond lengths, as well as the FFF angles, at the Hartree–Fock level. We then used our best geometries while calculating the Møllet–Plesset perturbation theory corrections through MP4 for N = 4 and 6, and through MP3 for N = 8. Further corrections for N = 8 did not seem justified, since agreement was so close between the other N = 8 results with the N = 6 ones.

Our results for these studies are summarized in Table , where we report both the total energies for each calculation as well as the predicted hydrogen-bond energy at each level of calculation. It should be noted that although our bond angles are in excellent agreement with the experimental values of 116 and 120°,[59] our final bond lengths are in rather poor agreement. In the crystal, the F–H bond length is found to be approximately 0.95 Å and the H···F distances are found to be around 1.53 Å.[59] Our results appear to be converging on nearly equal distances of 1.10 Å for both interactions. This is most likely an artifact of using a very incomplete basis set. The shortening of the distances is probably due to the system attempting to bring all of the basis functions closer together so that the electronic structure around the intramolecular bonds can be better described. However, it is the trends as a function of N that are our real concern with this model, and for these, we see that all three sets predict very similar geometries. For the energies, we have tabulated in parentheses the % errors relative to the MP4 results for N = 6. Here, we see a very interesting trend as we add correlation. First, we see that the N = 6 and 8 results are nearly identical to one another with the N = 4 results predicting slightly weaker hydrogen bonds. Most interesting is the trend as higher-order perturbations are included. The MP2 level seems to overcorrect toward predicting more stable H-bonds relative to the Hartree–Fock results, with the latter actually showing smaller errors than the MP2 values. The MP3 level values are all less than 1 kcal different from the final MP4 values. Thus, the greatest errors are seen at the MP2 level. Clearly, care must be taken when applying the Møllet–Plesset perturbation theory corrections. For systems where relatively strong interactions occur between monomers, if any correlation corrections are to be made using the Møllet–Plesset perturbation theory, caution must be exhibited before accepting corrections calculated only at the MP2 level.

Table 1

Calculated Geometries and Energies of HF Catemer

	Nx = 1	Nx = 4	Nx = 6	Nx = 8
geometries
RFH (Å)	0.942	1.016	1.080	1.080
RH···F (Å)		1.185	1.080	1.085
∠F···F···F (degree)		110.4	115.3	114.0
energies (hartrees)
ESCF	–99.531767	–99.554606	–99.557580	–99.558156
εMP2	–0.018435	–0.023157	–0.021758	–0.021057
εMP3	–0.005930	–0.003219	–0.002376	–0.002205
εMP4 - singles	–0.001996	–0.000426	–0.000170	–0.000065
εMP4 - doubles	–0.000116	0.000065	0.000329	0.000338
εMP4 - triples	0	–0.000514	–0.000263	–0.000142
εMP4 - quadruples	0	–0.000325	–0.000281
εMP4 - total	–0.002112	–0.001200	–0.000385
ΔEINT (kcal)
SCFa		14.33 (−4%)	16.20 (+8%)	16.56 (+11%)
MP2a		17.29 (+15%)	18.28 (+22%)	18.20 (+22%)
MP3a		15.59 (+4%)	16.05 (+7%)	15.87 (+6%)
MP4a,b		15.02 (0%)	14.97 (0%)	(14.46) (−3%)

% errors are given in parentheses after the ΔEINT values, relative to the MP4(N = 6) value of 14.97 kcal.

MP4 values were not calculated for N = 8, and thus that total is placed within parentheses.

Since our method is based upon a closed form for the wave functions, it is a straightforward process to add in correlation through the use of the Møller–Plesset perturbation theory. For the FH···F interactions in our model studies, we found that, although the MP2 level corrections led to a correction to the hydrogen-bond energies of around 1 kcal/mol, truncating the process at the MP2 level leads to considerably larger errors than those seen in the Hartree–Fock results. The optimal level appears to be at the MP3 level for these systems, since at all levels of calculation, these results were all well within 1 kcal of the MP4 results.

Conclusions

There are several advantages to using our cyclic periodic wave function (CPWF) approach for the treatment of infinitely periodic systems. First, this is a wave function-based method. Thus, it can be applied at any level of approximation—from semiempirical to full ab initio, with or without the inclusion of correlation and with any desired basis set.

Second, it is a chemically intuitive approach. The repeat units are chosen to be either the individual molecules within molecular solids or the individual monomer units within polymers,[60,61] and are placed at their exact translationally periodic positions, as found in the infinitely periodic crystal, without any built-in periodic boundary conditions (PBCs). For chemists, it is a more intuitive way to study large organic crystals. In this way, the CPWF approach differs from the cyclic cluster model (CCM), which places the atoms or molecules at cyclic periodic positions in space, and our approach is also different from all of the methods that apply PBCs. The zeroth-order wave function corresponds to a single isolated molecule or monomer unit. Increasing the periods for the repeat lengths (N, N, N), then corresponds to the inclusion of exact interactions with all nearest neighbors, when N = 4, and with all next-nearest neighbors as well, when N = 6, etc., and the ground-state wave function would ultimately converge to the exact Hartree–Fock ground state in the limit of N = ∞. Longer-range interactions are at least partially recovered by the use of our Madelung correction scheme,[40] which corrects both the individual Fock matrix elements as well as the total energy by approximately recovering all interactions between farther neighbors.[40] We have found that using N = 4 suffices for most solid-state systems,[60,61] as long as the intermonomer interactions are weaker than covalent interactions. Even when covalent interactions occur between monomers, we have found that N = 4 often suffices, though using N = 6 in the direction of the covalent interaction will provide greater accuracy for some systems. Without the inclusion of the Madelung corrections, higher values for the repeat length would be necessary to obtain equivalent accuracy.

Next, the electron density is treated as being infinitely periodic. Thus, bulk properties of solid-state systems or infinite polymers are directly calculable from the contribution of any single repeat unit in the system. Since every repeat unit will contribute an identical electron density as every other repeat unit, the chemical equivalence of the repeat units is maintained. Note that this is not the case for calculations that use a cluster approach. In cluster models, there will be many different chemical environments within the cluster, making it difficult, if not impossible, to determine contributions from any particular unit within the cluster.

The other advantage is that the number of integrals which must be calculated and the dimensions of matrices which must be diagonalized are tremendously reduced from those necessary using either cluster approaches or N-body approaches.

Because our wave functions use individual molecules or monomers as the repeat unit, instead of the unit cell of the crystal,[60,61] and because our wave functions are cyclically periodic, instead of translationally periodic, our method is origin independent and does not depend in any way upon how the repeat unit is chosen. The CPWF approach builds the symmetry relations between identical molecules or monomers directly into the wave functions, ensuring that the electron density at each individual molecule within the infinite system is identical.

An extension of this approach allows accurate inclusion of delocalization, delocalized systems can also be treated with the same method. Although the final equations were not derived here for DFT calculations, our wave function in eq can be used for DFT studies, as well as for ab initio and semiempirical calculations.

Finally, the CPWF approach is applied at the AM1 and PM3 semiempirical levels of approximation to study different three-dimensional infinitely periodic solid-state systems stabilized by different kind of weak interactions between repeat units, all intra- and intermonomer geometry parameters are optimizable in this CPWF method.[60,61] Thus, our approach leads to complete flexibility in the determination of the three-dimensional geometry of the solid-state system.