James A Platts
email: platts@cardiff.ac.uk

J Grant Hill

Interactions using Local Correlation Methods: Energy Partitioning, Geometry Optimisation and Harmonic Frequency Calculations

Keywords: Energy Partitioning, Geometry Optimisation and Harmonic Frequency Calculations Local correlation, Non-covalent interactions, Hydrogen bonding, Pistacking, Zero point enrgy

Non-covalent interactions play an important part in a wide range of chemical and biological phenomena, for instance folding of proteins and nucleic acids into their correct 3-dimensional shapes. 1 Intermolecular forces acting within molecular clusters are widely used as models of such non-covalent interactions, since their effects can be more easily separated from intra-molecular forces. Ab initio theoretical methods are increasingly used to complement experimental study of such clusters. The majority of theoretical treatments of non-covalent interactions employ the supermolecular approach, in which interaction energy is calculated as the difference in energy between the cluster and its constituent monomers. However, two significant drawbacks of this approach are the need for large orbital basis sets, including polarisation and diffuse functions, for proper description of all interactions between molecules, and the related problem of basis set superposition error (BSSE). 2 These problems can be ameliorated, to a large extent, by use of local correlation methods, which effectively eliminate BSSE by restricting excitations from occupied orbitals to those virtual orbitals that are in close spatial proximity. 3 Moreover, in combination with the density fitting (DF) approximation, methods such as DF-LMP2 4 and DF-LCCSD 5,6,7,8 make significant efficiency savings over conventional calculations, approaching linear scaling in favourable cases, hence allowing application of suitable orbital basis sets for description of non-covalent interactions. A further advantage of the local correlation approach is that intermolecular correlation energy can be partitioned into physically meaningful contributions, such as dispersion and ionic terms, on the basis of the occupied and virtual orbitals involved. 3 This approach is complementary to decomposition schemes, and is based solely on the local correlation treatment of the cluster and so requires no further calculations or significant computational resources. BSSE can also be accounted for by use of the counterpoise method, 9 but for a cluster of n molecules, this procedure requires (n+1) energy calculations in the full basis set of the cluster, substantially increasing the computational resources required. The situation becomes still less satisfactory if geometry optimisation of the cluster is a goal, as separate gradient calculations for the cluster and all fragments are necessary. Without such corrections, geometry optimisation is likely to result in underestimation of intermolecular separations and overestimation of stabilities.

Alternative approaches, such as symmetry adapted perturbation theory (SAPT) eliminate BSSE by construction, but are currently limited to pairs of molecules. 10 advantage in calculating the necessary data. Recently, Hrenar et al. described the use of DF-LMP2 methods for numerical evaluation of the Hessian using finite differences of analytical gradients, made possible by using a highly parallelised algorithm. 11 In doing so, they demonstrated the utility of this approach, and identified scaling factors that correct for at least some of the deficiencies of the harmonic approach.

Grimme's spin-component scaling (SCS) method, 12 and variants thereof, have been shown to correct for deficiencies in MP2 for stacked complexes, but may actually degrade performance for hydrogen bonded complexes. 13 Recently, Marchetti et al. proposed a dispersion-weighted SCS scheme that combines the accuracy of conventional MP2 for electrostatic complexes and of SCS-MP2 for dispersion bound ones. 14 In the current investigation, density fitted, local correlation methods are applied to a range of intermolecular interactions, encompassing complexes bound mainly by electrostatic effects, by dispersion forces, or by a combination of these. These complexes are taken from the recent review of Hobza and co-workers, 15 who suggest that they are representative of most important classes of non-covalent interaction. DF-LMP2, and DF-LCCSD where possible, methods are used to partition the interaction energy, to optimise geometries, and to calculate harmonic frequencies. In doing so, the aim is to show that DF-LMP2 and its SCS variant are capable of accurately and efficiently predicting the geometry and harmonic vibrational motion of complexes bound by non-covalent interactions, paving the way for studies of larger systems for which benchmark data is not available.

Computational Methods

All calculations were performed with the MOLPRO suite of programs, 16 employing the correlation consistent aug-cc-pVTZ orbital basis set throughout. 17,18 Density fitting of MP2 energies employed the corresponding MP2FIT basis set, 19 while the underlying Hartree-Fock calculation also employed density fitting 20 with the corresponding JKFIT basis. 21 DF in LCCSD for integrals over three and four external orbitals used the larger aug-cc-pVQZ/MP2FIT auxiliary basis sets to minimize fitting errors. The localised orbitals required for local correlation methods were generated via the Pipek-Mezey method, 22 while the orbital domain selection followed the procedure of Boughton and Pulay, 23 with merging of rotationally invariant π-domains where appropriate. Initial test calculations indicated that using domains calculated at large intermolecular separation and then frozen for calculations on the interacting system produced almost negligible differences in geometry and harmonic frequencies when compared with domains computed solely for the interacting system (so-called floating domains). For example, with the small methane dimer system a maximum difference of 6 cm -1 between normal modes was observed, while for the larger parallel-displaced benzene dimer system the maximum difference was 0.2 cm -1 . As deviations of this size are likely to be small compared with both the basis set and method errors, floating domains were used for all geometry optimisations and calculation of numerical harmonic frequencies.

Analytical gradients for DF-LMP2 methods 24 are available in MOLPRO, enabling efficient geometry optimisation. All geometry optimisations reported commenced from the reported literature structure, and employed default criteria for convergence. Calculation of root mean square (RMS) deviation between initial and final structures used Chemcraft. 25 Following geometry optimisation, harmonic frequencies were calculated by numerical differentiation of analytic gradients using central differences, employing the default step size of 0.01 a.u. This is a sizeable task for the larger molecules considered, but was made feasible by use of a parallel algorithm, 26 as implemented in MOLPRO, in which individual energy calculations on different geometry displacements used in the formation of gradients and hessians are performed independently on separate processors. Individual DF-LMP2 calculations did not take account of molecular symmetry (as this usually leads to symmetry equivalent localised orbitals), but this numerical algorithm can use symmetry to reduce the number of displacements required. The local character of the orbitals allows for the partitioning of the intermolecular interaction energy into intramolecular, dispersive, and ionic components of the correlation energy. 3 This is achieved by dividing the system into monomers and separating the excitations from localised orbitals into different classes based on the domains the pair of electrons are excited into.

Results and Discussion

Table 1 presents the results of partitioning the DF-LMP2 correlation energy into its intermolecular components, along with the interaction energy at the SCF level, for the S22 set of complexes of Hobza and co-workers. 15 Following reference 3, dispersion refers to simultaneous excitations on different monomers, ionic to excitation from an occupied orbital on one monomer to a virtual orbital on the other along with a single excitation on the acceptor monomer, and dispersionexchange to cross-excitations. Such definitions stem naturally from the local correlation approach, but may not be directly comparable to assignments from methods such as SAPT. 10 Intramolecular correlation can also be calculated in this fashion, but does not show significant differences from the sum of monomer correlation energies. The three intermolecular correlation terms are taken directly from the DF-LMP2 cluster calculation, while the SCF interaction energy was calculated using the supermolecular approach. In all cases, dispersion-exchange energies are essentially negligible, and are included only for completeness. The next eight complexes were denoted as predominantly bound through dispersion, 15 an assignment supported by the present analysis. In all eight clusters, SCF energies are repulsive, with all stabilisation stemming from correlation effects. Unlike the hydrogen bonded complexes, dispersion effects are much greater (typically by a factor of two to five) than ionic ones in these complexes, more than counteracting the repulsive SCF term. This is exemplified by the archetypal stacked, parallel-displaced benzene dimer, for which dispersion effects are 2.5 times greater than ionic ones. The final seven complexes are of so-called "mixed" type, with no single contribution dominating. This designation is again supported by the current analysis; SCF energies are attractive in three cases and repulsive in the remaining four. Dispersion contributions are generally larger than ionic ones, but not to the same extent as for the dispersion-bound clusters. Complexes that could be designated as containing X-H...π (X = C, N, O) non-covalent bonds, such as benzene...water, benzene...ammonia and benzene...HCN, show significant differences from the conventional hydrogen bonds. Most notably, dispersion contributions are significantly greater than SCF ones. The clear contrast in the origin of stabilisation between the three classes of complex is illustrated in Figure 1. 1. Vertical dashed lines show separation of classes of complex used in ref. 15. The numbering on the abscissa corresponds to the ordering of the complexes in Table 1.

The analysis of Table 1 is based on MP2 methods, known to perform adequately for H-bonding but to substantially over-estimate interaction energies of stacked complexes. Analogous partitioning of correlation energy with the DF-LCCSD method (energy partitioning of perturbative triples contributions is not currently possible) has also been performed. For inclusion of all relevant electron pairs in the LCCSD treatment, modification of the default distance based domain selection criteria is required: after some trial and error, all orbital pairs within a distance of 8 a.u. are treated as strong pairs, i.e., at the CCSD level. Unfortunately, this renders DF-LCCSD calculations of the larger complexes unfeasible with a moderately large basis set such as aug-cc-pVTZ.

The partitioning of the DF-LCCSD correlation energy is presented in Table 2 for selected, smaller complexes, and shows that the trends observed using DF-LMP2 methods are largely preserved at this improved level of theory. For dispersion bound and mixed complexes, the tendency of MP2 to overestimate interaction energies is evident in both dispersion and ionic terms, with rather larger reductions in the former than the latter. This is exemplified in data for the ethene dimer, for which dispersion contributions are reduced by 0. This can be a laborious process, and might be significantly eased by the advantages of the DF-LMP2 method discussed above, particularly the removal of BSSE. Table 3 reports the results of full, unconstrained geometry optimisation of the same S22 complexes using DF-LMP2/aug-cc-pVTZ, starting from the literature geometry. Proper application of local correlation methods requires symmetry to be disabled, but the optimisation algorithm retains any elements of symmetry present in the initial geometry.

Optimisation generally proceeded smoothly, in most cases reaching the default convergence criteria in MOLPRO (maximum gradient below 3 x 10 -4 a.u.) in less than 10 steps. Table 3 demonstrates that the resulting optimised geometries are very similar to their literature equivalents in all but one case, that being the phenol dimer. Performance for hydrogen bonded and dispersion bound complexes is similar, whereas the slightly larger RMS deviation for the mixed complexes is almost entirely due to the phenol dimer (omitting this point gives a mean RMS of 0.018 Å for the remaining six complexes). No significant difference in the quality of agreement is observed between MP2 and CCSD(T) optimised literature geometries. aromatic rings into slightly closer contact (distance between centroids is 5.12 Å vs. 5.39 Å in the literature structure). This subtle change notwithstanding, the overall conclusion is that this efficient and counterpoise-free approach to geometry optimisation satisfactorily and consistently reproduces conventional, corrected geometries. Both conventional MP2 and DF-LMP2 are well known to overestimate the interaction energy of stacked complexes, and might therefore be expected to underestimate intermolecular separation in such complexes. This is apparent in potential energy curves for the stacked benzene dimer, which reach a minimum at vertical separation of ca. 3.5 Å at CBS(T) level and ca. 3.35 Å at MP2. 27 Recently, it was demonstrated that a new parameterisation of SCS designed specifically for noncovalent interactions and denoted SCSN, 13a overcomes this deficiency, resulting in high-accuracy interaction energies for all classes of complex with errors of less than 0.5 kcal mol -1 . In this approach, correlation energy from same-and opposite-spin electron pairs are weighted by scaling factors in order to correct for systematic errors in standard MP2. The same scaling factors can be applied to calculation of analytic gradients, allowing optimisation to minima on this scaled potential energy surface. Table 3 also reports the RMS deviations from literature geometry using this approach.

In general, SCSN optimised structures are slightly further from literature values than unscaled DF-LMP2, both for MP2 and CCSD(T) reference data. Overall performance is similar to unscaled optimisation for hydrogen bonded and mixed complexes, but differences are much larger for dispersion bound complexes. The relatively small RMS deviations reported in A second SCSN geometry optimisation was performed for the T-shaped benzene dimer, following the recent study by Wang et al. that showed the lowest energy geometry of this complex to deviate from the ideal C 2v geometry. 28 Starting from a slightly offset geometry, optimisation proceeded smoothly to a point just 0.06 kcal mol -1 lower in energy than that reported in Table 3. Shown in Figure 3, the SCSN geometry is similar to that reported in ref. 28, with the donor C-H directed towards a C-C bond rather than the centre of the acceptor ring. optimised structures are true minima at this theoretical level, as confirmed by the lack of imaginary frequencies. As noted above, any molecular symmetry present in the literature geometry was retained in these calculations, such that our data lends support to the symmetries employed in the original study of Hobza and co-workers. This is in marked contrast to a recent DFT-based study of the same set of molecules, in which one or more imaginary frequencies were observed for several complexes. 29 The only exception is the T-shaped benzene dimer, for which the high-symmetry (C 2v ) form is found to be a transition state, with a single imaginary frequency of -10.8 cm -1 . Reoptimisation with lower (C s ) symmetry results in a geometry similar to that shown in Figure 3, which is found to be a true minimum on the potential energy surface.

A similar procedure was followed for the SCSN-optimised geometries. In this case, both the C 2v and C s optimised forms of the T-shaped benzene dimer are found to be true minima, echoing the findings of Wang et al, 28 who used dispersion-corrected DFT to show that both forms are stable and are connected by a transition state with a barrier of just 0.01 kcal mol -1 . In this regard, SCSN again seems to be slightly more accurate than unscaled DF-LMP2 in finding both forms to be stable, but the inherent errors associated with this approach are much larger than these energy differences, such that any conclusion based on this data must be tentative at best. All remaining structures are found to be minima using this approach. Two geometries of the stacked (paralleldisplaced) benzene dimer were tested, one with a C-H bond situated over the partner benzene and a second with a C-C bond over the other benzene. Both are found to be minima using SCSN, with effectively identical energies within the expected accuracy of the method (less than 0.1 kcal mol -1 difference).

As well as confirming structures as minima, harmonic frequency calculations can be used to estimate the change in zero-point vibrational energy (ZPVE) that accompanies formation of a complex from constituent monomers. Such changes stem largely from conversion of translational and rotational degrees of freedom of individual molecules into intermolecular vibrations. Such vibrations are likely to be strongly anharmonic in nature, such that harmonic data can only be a first-order approximation to the true ZPVE correction. Nonetheless, obtaining even harmonic vibrational data for the larger complexes considered here is a considerable challenge, and the data presented in general rather smaller than for hydrogen bonded ones. For example, ∆ZPVE for the stacked uracil dimer is less than half the value for the hydrogen bonded isomer, with a similar pattern for the adenine-thymine complexes. A slightly surprising exception to this pattern is the methane dimer, whose ∆ZPVE value of 1.16 kcal mol -1 is more than twice the binding energy. This may represent the breakdown of the harmonic approximation, or be a result of noise introduced by the use of numerical differences in such a weakly bound case.

In addition to raw DF-LMP2 data, Table 4 also contains ∆ZPVE values scaled by the value of 0.956 recommended by Hrenar et al, 11 and values from numerical frequency calculation based on SCSN energies. Changes using the single scaling factor are generally very small, and do not significantly alter the conclusions drawn from unscaled data. A similar pattern emerges from most of the SCSN data, with one exception: in the stacked indole-benzene complex, the SCSN value is almost double that from unscaled DF-LMP2. In the absence of benchmark data, it is not possible to decide which of these values is the more reliable, but in the context of other stacked complexes the raw DF-LMP2 value seems rather small.

Conclusions

The utility of local correlation methods, when combined with density fitting, has been demonstrated in the study of non-covalent interactions. By construction, these methods allow for decomposition of correlation energy into contributions from dispersion and ionic excitations, thus giving insight into the origin of stabilisation at no additional computational cost to the correlation calculation. This analysis, at both DF-LMP2 and DF-LCCSD levels, clearly distinguishes complexes bound by hydrogen bonding, dispersion, or a mixture of these effects. The lack of basis set superposition error in local methods, coupled with their computational efficiency, makes them ideally suited for simple geometry optimisation and harmonic frequency calculation. Geometry optimisation leads to structures in close agreement with published, counterpoise-corrected geometries, demonstrating the utility of this approach. Harmonic frequency analysis, using both unscaled and SCSN-scaled MP2 data, indicates that most complexes within the S22 set are at least local minima on their respective potential energy surfaces. Non-covalent interactions play an important part in a wide range of chemical and biological phenomena, for instance folding of proteins and nucleic acids into their correct 3-dimensional shapes. 1 Intermolecular forces acting within molecular clusters are widely used as models of such non-covalent interactions, since their effects can be more easily separated from intra-molecular forces. Ab initio theoretical methods are increasingly used to complement experimental study of such clusters. The majority of theoretical treatments of non-covalent interactions employ the supermolecular approach, in which interaction energy is calculated as the difference in energy between the cluster and its constituent monomers. However, two significant drawbacks of this approach are the need for large orbital basis sets, including polarisation and diffuse functions, for proper description of all interactions between molecules, and the related problem of basis set superposition error (BSSE). 2 These problems can be ameliorated, to a large extent, by use of local correlation methods, which effectively eliminate BSSE by restricting excitations from occupied orbitals to those virtual orbitals that are in close spatial proximity. 3 Moreover, in combination with the density fitting (DF) approximation, methods such as DF-LMP2 4 and DF-LCCSD 5,6,7,8 make significant efficiency savings over conventional calculations, approaching linear scaling in favourable cases, hence allowing application of suitable orbital basis sets for description of non-covalent interactions. A further advantage of the local correlation approach is that intermolecular correlation energy can be partitioned into physically meaningful contributions, such as dispersion and ionic terms, on the basis of the occupied and virtual orbitals involved. 3 This approach is complementary to decomposition schemes, and is based solely on the local correlation treatment of the cluster and so requires no further calculations or significant computational resources.

BSSE can also be accounted for by use of the counterpoise method, 9 but for a cluster of n molecules, this procedure requires (n+1) energy calculations in the full basis set of the cluster, substantially increasing the computational resources required. The situation becomes still less satisfactory if geometry optimisation of the cluster is a goal, as separate gradient calculations for the cluster and all fragments are necessary. Without such corrections, geometry optimisation is likely to result in underestimation of intermolecular separations and overestimation of stabilities.

Alternative approaches, such as symmetry adapted perturbation theory (SAPT) eliminate BSSE by construction, but are currently limited to pairs of molecules. methods for numerical evaluation of the Hessian using finite differences of analytical gradients, made possible by using a highly parallelised algorithm. 11 In doing so, they demonstrated the utility of this approach, and identified scaling factors that correct for at least some of the deficiencies of the harmonic approach.

Grimme's spin-component scaling (SCS) method, 12 and variants thereof, have been shown to correct for deficiencies in MP2 for stacked complexes, but may actually degrade performance for hydrogen bonded complexes. 13 Recently, Marchetti et al. proposed a dispersion-weighted SCS scheme that combines the accuracy of conventional MP2 for electrostatic complexes and of SCS-MP2 for dispersion bound ones. 14 In the current investigation, density fitted, local correlation methods are applied to a range of intermolecular interactions, encompassing complexes bound mainly by electrostatic effects, by dispersion forces, or by a combination of these. These complexes are taken from the recent review of Hobza and co-workers, 15 who suggest that they are representative of most important classes of non-covalent interaction. DF-LMP2, and DF-LCCSD where possible, methods are used to partition the interaction energy, to optimise geometries, and to calculate harmonic frequencies. In doing so, the aim is to show that DF-LMP2 and its SCS variant are capable of accurately and efficiently predicting the geometry and harmonic vibrational motion of complexes bound by non-covalent interactions, paving the way for studies of larger systems for which benchmark data is not available.

Computational Methods

All calculations were performed with the MOLPRO suite of programs, 16 employing the correlation consistent aug-cc-pVTZ orbital basis set throughout. 17,18 Density fitting of MP2 energies employed the corresponding MP2FIT basis set, 19 while the underlying Hartree-Fock calculation also employed density fitting 20 with the corresponding JKFIT basis. 21 DF in LCCSD for integrals over three and four external orbitals used the larger aug-cc-pVQZ/MP2FIT auxiliary basis sets to minimize fitting errors. The localised orbitals required for local correlation methods were generated via the Pipek-Mezey method, 22 while the orbital domain selection followed the procedure of Boughton and Pulay, 23 with merging of rotationally invariant π-domains where appropriate. Initial test calculations indicated that using domains calculated at large intermolecular separation and then frozen for calculations on the interacting system produced almost negligible differences in geometry and harmonic frequencies when compared with domains computed solely for the interacting system (so-called floating domains). For example, with the small methane dimer system a maximum difference of 6 cm -1 between normal modes was observed, while for the larger parallel-displaced benzene dimer system the maximum difference was 0.2 cm -1 . As deviations of this size are likely to be small compared with both the basis set and method errors, floating domains were used for all geometry optimisations and calculation of numerical harmonic frequencies.

Analytical gradients for DF-LMP2 methods 24 are available in MOLPRO, enabling efficient geometry optimisation. All geometry optimisations reported commenced from the reported literature structure, and employed default criteria for convergence. Calculation of root mean square (RMS) deviation between initial and final structures used Chemcraft. 25 Following geometry optimisation, harmonic frequencies were calculated by numerical differentiation of analytic gradients using central differences, employing the default step size of 0.01 a.u. This is a sizeable task for the larger molecules considered, but was made feasible by use of a parallel algorithm, 26 as implemented in MOLPRO, in which individual energy calculations on different geometry displacements used in the formation of gradients and hessians are performed independently on separate processors. Individual DF-LMP2 calculations did not take account of molecular symmetry (as this usually leads to symmetry equivalent localised orbitals), but this numerical algorithm can use symmetry to reduce the number of displacements required. The local character of the orbitals allows for the partitioning of the intermolecular interaction energy into intramolecular, dispersive, and ionic components of the correlation energy. 3 This is achieved by dividing the system into monomers and separating the excitations from localised orbitals into different classes based on the domains the pair of electrons are excited into.

Results and Discussion

Table 1 presents the results of partitioning the DF-LMP2 correlation energy into its intermolecular components, along with the interaction energy at the SCF level, for the S22 set of complexes of Hobza and co-workers. 15 Following reference 3, dispersion refers to simultaneous excitations on different monomers, ionic to excitation from an occupied orbital on one monomer to a virtual orbital on the other along with a single excitation on the acceptor monomer, and dispersionexchange to cross-excitations. Such definitions stem naturally from the local correlation approach, but may not be directly comparable to assignments from methods such as SAPT. 10 Intramolecular correlation can also be calculated in this fashion, but does not show significant differences from the sum of monomer correlation energies. The three intermolecular correlation terms are taken directly from the DF-LMP2 cluster calculation, while the SCF interaction energy was calculated using the supermolecular approach. In all cases, dispersion-exchange energies are essentially negligible, and are included only for completeness. The first seven clusters, denoted as hydrogen bonded in the original work, all display significant stabilisation at the SCF level, reflecting the importance of electrostatics in such interactions.

Nonetheless, dispersion and ionic correlation effects contribute significantly to the overall stabilisation, with approximately equal importance. In all seven cases, the SCF stabilisation is the largest single term, but in three (ammonia dimer, 2-pyridoxine...2-aminopyridine and adenine...thymine) the combined correlation terms outweigh the SCF contribution. For example, in Watson-Crick hydrogen bonded adenine...thymine the SCF, dispersion and ionic contributions are 39%, 35% and 29% of the overall binding energy, respectively.

The next eight complexes were denoted as predominantly bound through dispersion, 15 an assignment supported by the present analysis. In all eight clusters, SCF energies are repulsive, with all stabilisation stemming from correlation effects. Unlike the hydrogen bonded complexes, dispersion effects are much greater (typically by a factor of two to five) than ionic ones in these complexes, more than counteracting the repulsive SCF term. This is exemplified by the archetypal stacked, parallel-displaced benzene dimer, for which dispersion effects are 2.5 times greater than ionic ones. The final seven complexes are of so-called "mixed" type, with no single contribution dominating. This designation is again supported by the current analysis; SCF energies are attractive in three cases and repulsive in the remaining four. Dispersion contributions are generally larger than ionic ones, but not to the same extent as for the dispersion-bound clusters. Complexes that could be designated as containing X-H...π (X = C, N, O) non-covalent bonds, such as benzene...water, benzene...ammonia and benzene...HCN, show significant differences from the conventional hydrogen bonds. Most notably, dispersion contributions are significantly greater than SCF ones. The clear contrast in the origin of stabilisation between the three classes of complex is illustrated in Figure 1.

Figure 1

Illustration of energy partitioning data from Table 1. Vertical dashed lines show separation of classes of complex used in ref. 15. The numbering on the abscissa corresponds to the ordering of the complexes in Table 1. The analysis of Table 1 is based on MP2 methods, known to perform adequately for H-bonding but to substantially over-estimate interaction energies of stacked complexes. Analogous partitioning of correlation energy with the DF-LCCSD method (energy partitioning of perturbative triples contributions is not currently possible) has also been performed. For inclusion of all relevant electron pairs in the LCCSD treatment, modification of the default distance based domain selection criteria is required: after some trial and error, all orbital pairs within a distance of 8 a.u. are treated as strong pairs, i.e., at the CCSD level. Unfortunately, this renders DF-LCCSD calculations of the larger complexes unfeasible with a moderately large basis set such as aug-cc-pVTZ.

The partitioning of the DF-LCCSD correlation energy is presented in Table 2 for selected, smaller complexes, and shows that the trends observed using DF-LMP2 methods are largely preserved at this improved level of theory. For dispersion bound and mixed complexes, the tendency of MP2 to overestimate interaction energies is evident in both dispersion and ionic terms, with rather larger reductions in the former than the latter. This is exemplified in data for the ethene dimer, for which dispersion contributions are reduced by 0.1 kcal mol -1 and ionic by 0.02 kcal mol -1 . This compares well with the ∆CCSD(T) correction to the overall interaction of +0.09 reported in ref 15. This can be a laborious process, and might be significantly eased by the advantages of the DF-LMP2 method discussed above, particularly the removal of BSSE. Table 3 reports the results of full, unconstrained geometry optimisation of the same S22 complexes using DF-LMP2/aug-cc-pVTZ, starting from the literature geometry. Proper application of local correlation methods requires symmetry to be disabled, but the optimisation algorithm retains any elements of symmetry present in the initial geometry.

Optimisation generally proceeded smoothly, in most cases reaching the default convergence criteria in MOLPRO (maximum gradient below 3 x 10 -4 a.u.) in less than 10 steps. Table 3 demonstrates that the resulting optimised geometries are very similar to their literature equivalents in all but one case, that being the phenol dimer. Performance for hydrogen bonded and dispersion bound complexes is similar, whereas the slightly larger RMS deviation for the mixed complexes is almost entirely due to the phenol dimer (omitting this point gives a mean RMS of 0.018 Å for the remaining six complexes). No significant difference in the quality of agreement is observed between MP2 and CCSD(T) optimised literature geometries. Both conventional MP2 and DF-LMP2 are well known to overestimate the interaction energy of stacked complexes, and might therefore be expected to underestimate intermolecular separation in such complexes. This is apparent in potential energy curves for the stacked benzene dimer, which reach a minimum at vertical separation of ca. interaction energies for all classes of complex with errors of less than 0.5 kcal mol -1 . In this approach, correlation energy from same-and opposite-spin electron pairs are weighted by scaling factors in order to correct for systematic errors in standard MP2. The same scaling factors can be applied to calculation of analytic gradients, allowing optimisation to minima on this scaled potential energy surface. Table 3 also reports the RMS deviations from literature geometry using this approach.

In general, SCSN optimised structures are slightly further from literature values than unscaled DF-LMP2, both for MP2 and CCSD(T) reference data. Overall performance is similar to unscaled optimisation for hydrogen bonded and mixed complexes, but differences are much larger for dispersion bound complexes. The relatively small RMS deviations reported in Table 3 mask some rather more significant changes; for instance, in the stacked benzene dimer, the separation of mean planes is 3.361 Å in the literature (MP2/TZ-CP optimised) geometry, whereas our SCSN optimisation results in a separation of 3.487 Å, much closer to the minimum of the CBS(T) curve. 27 A similar trend is apparent in other stacked complexes, albeit to a lesser extent: in the pyrazine dimer, the molecular planes are separated by 3.261 Å in literature geometry and 3.270 Å in the SCSN geometry. In stacked adenine-thymine, molecular planes are not quite parallel, so for the current investigation the perpendicular distance is defined as that from the mean plane of adenine to the centroid of the 6-membered ring of thymine, which is 3.174 Å in literature and 3.205 Å in SCSN optimised geometry. Thus, it is evident that the excess stabilisation of stacked complexes at MP2 level carries through into potential energy surfaces, and can be simply and efficiently corrected by the empirical SCSN scheme.

A second SCSN geometry optimisation was performed for the T-shaped benzene dimer, following the recent study by Wang et al. that showed the lowest energy geometry of this complex to deviate from the ideal C 2v geometry. 28 Starting from a slightly offset geometry, optimisation proceeded smoothly to a point just 0.06 kcal mol -1 lower in energy than that reported in Table 3. Shown in Figure 3, the SCSN geometry is similar to that reported in ref. 28, with the donor C-H directed towards a C-C bond rather than the centre of the acceptor ring. The fully optimised DF-LMP2 geometries were then used to calculate harmonic frequencies using numerical differentiation. Firstly, and most importantly, this procedure confirms that all but one optimised structures are true minima at this theoretical level, as confirmed by the lack of imaginary frequencies. As noted above, any molecular symmetry present in the literature geometry was retained in these calculations, such that our data lends support to the symmetries employed in the original study of Hobza and co-workers. This is in marked contrast to a recent DFT-based study of the same set of molecules, in which one or more imaginary frequencies were observed for several complexes. 29 The only exception is the T-shaped benzene dimer, for which the high-symmetry (C 2v ) form is found to be a transition state, with a single imaginary frequency of -10.8 cm -1 . Reoptimisation with lower (C s ) symmetry results in a geometry similar to that shown in Figure 3, which is found to be a true minimum on the potential energy surface.

A similar procedure was followed for the SCSN-optimised geometries. In this case, both the C 2v and C s optimised forms of the T-shaped benzene dimer are found to be true minima, echoing the findings of Wang et al, 28 who used dispersion-corrected DFT to show that both forms are stable and are connected by a transition state with a barrier of just 0.01 kcal mol -1 . In this regard, SCSN again seems to be slightly more accurate than unscaled DF-LMP2 in finding both forms to be stable, but the inherent errors associated with this approach are much larger than these energy differences, such that any conclusion based on this data must be tentative at best. All remaining structures are found to be minima using this approach. Two geometries of the stacked (paralleldisplaced) benzene dimer were tested, one with a C-H bond situated over the partner benzene and a second with a C-C bond over the other benzene. Both are found to be minima using SCSN, with effectively identical energies within the expected accuracy of the method (less than 0.1 kcal mol -1 difference).

As well as confirming structures as minima, harmonic frequency calculations can be used to estimate the change in zero-point vibrational energy (ZPVE) that accompanies formation of a complex from constituent monomers. Such changes stem largely from conversion of translational and rotational degrees of freedom of individual molecules into intermolecular vibrations. Such vibrations are likely to be strongly anharmonic in nature, such that harmonic data can only be a first-order approximation to the true ZPVE correction. Nonetheless, obtaining even harmonic vibrational data for the larger complexes considered here is a considerable challenge, and the data presented in a Scaled from raw DF-LMP2 data using the recommended scaling factor of 0.956 ∆ZPVE values on formation of complexes from monomers are largest for the hydrogen bonded complexes, and are in general rather smaller for dispersion bound and mixed complexes. The single largest change is found for the formamide dimer, closely followed by water and formic acid dimers. The value for the water dimer is similar to previous estimates using MP2 methods with the harmonic approximation (2.1 -2.3 kcal mol -1 ), and as expected rather larger than values that take anharmonicity into account (1.8 -2.1 kcal mol -1 ). 30 The larger aromatic systems containing two general rather smaller than for hydrogen bonded ones. For example, ∆ZPVE for the stacked uracil dimer is less than half the value for the hydrogen bonded isomer, with a similar pattern for the adenine-thymine complexes. A slightly surprising exception to this pattern is the methane dimer, whose ∆ZPVE value of 1.16 kcal mol -1 is more than twice the binding energy. This may represent the breakdown of the harmonic approximation, or be a result of noise introduced by the use of numerical differences in such a weakly bound case.

In addition to raw DF-LMP2 data, Table 4 also contains ∆ZPVE values scaled by the value of 0.956 recommended by Hrenar et al, 11 and values from numerical frequency calculation based on SCSN energies. Changes using the single scaling factor are generally very small, and do not significantly alter the conclusions drawn from unscaled data. A similar pattern emerges from most of the SCSN data, with one exception: in the stacked indole-benzene complex, the SCSN value is almost double that from unscaled DF-LMP2. In the absence of benchmark data, it is not possible to decide which of these values is the more reliable, but in the context of other stacked complexes the raw DF-LMP2 value seems rather small.

Conclusions

The utility of local correlation methods, when combined with density fitting, has been demonstrated in the study of non-covalent interactions. By construction, these methods allow for decomposition of correlation energy into contributions from dispersion and ionic excitations, thus giving insight into the origin of stabilisation at no additional computational cost to the correlation calculation. This analysis, at both DF-LMP2 and DF-LCCSD levels, clearly distinguishes complexes bound by hydrogen bonding, dispersion, or a mixture of these effects. The lack of basis set superposition error in local methods, coupled with their computational efficiency, makes them ideally suited for simple geometry optimisation and harmonic frequency calculation. Geometry optimisation leads to structures in close agreement with published, counterpoise-corrected geometries, demonstrating the utility of this approach. Harmonic frequency analysis, using both unscaled and SCSN-scaled MP2 data, indicates that most complexes within the S22 set are at least local minima on their respective potential energy surfaces. The only exception to this conclusion is the C 2v T-shaped benzene dimer, which is a minimum with SCSN scaling but a saddle point with unscaled MP2. Zero-point vibrational energy corrections to binding energies have also been

Figure 1

 1 Figure 1 Illustration of energy partitioning data from Table 1. Vertical dashed lines show

Figure 2

 2 shows the literature and optimised geometries of the phenol dimer, which by inspection are highly similar. The hydrogen bond O...H distances in the two cases are almost identical (1.937 vs. 1.935 Å). However, the dihedral angle between the mean planes of each molecule are different, with values of 60.5º in the literature geometry and 56.3º in the DF-LMP2 optimised structure. Thus, DF-LMP2 brings the

aFigure 2

 2 Figure 2 Phenol dimer a) literature geometry; b) DF-LMP2 optmised.

Figure 3

 3 Figure 3 SCSN optimised geometry of T-shaped benzene dimer

10

 10

  Molecular clusters are dynamic objects, and understanding of their vibrational properties is key to understanding their thermodynamic and spectroscopic properties. The near-universal starting point for theoretical description of vibrational behaviour is the harmonic approximation, which, despite its known shortcomings, is widely used to interpret thermodynamic and spectroscopic properties of interest. The efficiency and lack of BSSE in local correlation methods is again a significant advantage in calculating the necessary data. Recently, Hrenar et al. described the use of DF-LMP2

Figure 2

 2 shows the literature and optimised geometries of the phenol dimer, which by inspection are highly similar. The hydrogen bond O...H distances in the two cases are almost identical (1.937 vs. 1.935 Å). However, the dihedral angle between the mean planes of each molecule are different, with values of 60.5º in the literature geometry and 56.3º in the DF-LMP2 optimised structure. Thus, DF-LMP2 brings the aromatic rings into slightly closer contact (distance between centroids is 5.12 Å vs. 5.39 Å in the literature structure). This subtle change notwithstanding, the overall conclusion is that this efficient and counterpoise-free approach to geometry optimisation satisfactorily and consistently reproduces conventional, corrected geometries.

Figure 2

 2 Figure 2 Phenol dimer a) literature geometry; b) DF-LMP2 optmised.

Figure 3

 3 Figure 3 SCSN optimised geometry of T-shaped benzene dimer

  show much smaller ∆ZPVE, falling to ca. 1 kcal mol -1 for 2-pyridoxine...2aminopyridine and adenine...thymine. Values for dispersion bound and mixed complexes are in

Table 1

 1 SCF interaction energies and partitioned DF-LMP2 intermolecular energies (kcal mol-1 )

		SCF	Dispersion	Ionic	Dispersion-
					Exchange
	Ammonia dimer	-1.27	-1.23	-0.77	-0.01
	Water dimer	-3.72	-1.06	-1.23	+0.01
	Formic acid dimer	-15.49	-4.20	-5.63	+0.01
	Formamide dimer	-12.32	-3.73	-4.09	+0.08
	Uracil dimer HB	-9.31	-4.82	-4.79	+0.07
	2-Pyridoxine...2-aminopyridine	-5.81	-5.52	-4.76	-0.03
	Adenine…Thymine HB	-6.40	-5.68	-4.79	+0.01
	Methane dimer	+0.35	-0.73	-0.14	0.00
	Ethene dimer	+0.79	-1.88	-0.60	+0.01
	Benzene…methane	+1.60	-2.19	-0.86	+0.02
	Benzene dimer stack	+6.33	-7.22	-2.78	+0.02
	Pyrazine dimer	+7.60	-7.77	-3.11	+0.03
	Uracil dimer stack	+5.20	-9.15	-3.13	+0.01
	Indole…benzene stack	+8.62	-10.60	-4.46	+0.11
	Adenine…Thymine stack	+5.72	-14.05	-5.42	+0.05
	Ethene…ethyne	-0.52	-1.01	-0.58	0.01
	Benzene…water	-0.35	-2.01	-1.05	+0.01
	Benzene…ammonia	+0.74	-2.10	-0.95	+0.01
	Benzene…cyanide	-1.55	-2.75	-1.63	+0.04
	Benzene dimer T	+2.54	-3.84	-1.57	+0.05
	Indole…benzene T	+2.00	-5.13	-2.55	+0.05

Table 2

 2 Partitioned DF-LCCSD interaction energies for selected complexes (kcal mol-1 )

	1 kcal mol -1 and ionic by 0.02 kcal mol -1 . This compares

CCSD(T) optimisation with triple-or quadruple-ζ basis sets and including counterpoise correction.

Table 3

 3 RMS deviations from literature geometry (Å)

	Name	DF-LMP2	SCSN
	Ammonia dimer a	0.025	0.027
	Water dimer a	0.030	0.037
	Formic acid dimer a	0.017	0.016
	Formamide dimer a	0.025	0.028
	Uracil dimer HB	0.008	0.026
	2-Pyridoxine...2-aminopyridine	0.009	0.034
	Adenine…Thymine HB	0.011	0.028
	Methane dimer a	0.005	0.050
	Ethene dimer a	0.017	0.049
	Benzene…methane	0.005	0.039
	Benzene dimer stack	0.029	0.068
	Pyrazine dimer	0.021	0.048
	Uracil dimer stack	0.023	0.053
	Indole…benzene stack	0.017	0.039
	Adenine…Thymine stack	0.018	0.035
	Ethene…ethyne a	0.014	0.017
	Benzene…water	0.026	0.012
	Benzene…ammonia	0.019	0.021
	Benzene…cyanide	0.019	0.016
	Benzene dimer T	0.015	0.027
	Indole…benzene T	0.012	0.023
	Phenol dimer	0.183	0.173
	Mean	0.025	0.039
	Mean H-bond	0.018	0.028
	Mean dispersion	0.017	0.048

Table 3

 3 rather more significant changes; for instance, in the stacked benzene dimer, the separation of mean planes is 3.361 Å in the literature (MP2/TZ-CP optimised) geometry, whereas our SCSN optimisation results in a separation of 3.487 Å, much closer to the minimum of the CBS(T) curve.27 A similar trend is apparent in other stacked complexes, albeit to a lesser extent: in the pyrazine dimer, the molecular planes are separated by 3.261 Å in literature geometry and 3.270 Å in the SCSN geometry. In stacked adenine-thymine, molecular planes are not quite parallel, so for the current investigation the perpendicular distance is defined as that from the mean plane of adenine to the centroid of the 6-membered ring of thymine, which is 3.174 Å in literature and 3.205 Å in SCSN optimised geometry. Thus, it is evident that the excess stabilisation of stacked complexes at MP2 level carries through into potential energy surfaces, and can be simply and efficiently corrected by the empirical SCSN scheme.

	mask some

Table 4

 4

	is a step towards establishing accurate ZPVE corrections for systems of this
	size.

Table 4

 4 Harmonic zero point energy corrections for binding (kcal mol -1 ) Scaled from raw DF-LMP2 data using the recommended scaling factor of 0.956 ∆ZPVE values on formation of complexes from monomers are largest for the hydrogen bonded complexes, and are in general rather smaller for dispersion bound and mixed complexes. The single largest change is found for the formamide dimer, closely followed by water and formic acid dimers. The value for the water dimer is similar to previous estimates using MP2 methods with the harmonic approximation (2.1 -2.3 kcal mol -1 ), and as expected rather larger than values that take aminopyridine and adenine...thymine. Values for dispersion bound and mixed complexes are in

	Name	DF-LMP2	Scaled a	SCSN
	Ammonia dimer	1.70	1.62	1.66
	Water dimer	2.37	2.27	2.33
	Formic acid dimer	2.08	1.99	2.06
	Formamide dimer	2.64	2.53	2.52
	Uracil dimer HB	1.21	1.16	1.15
	2-Pyridoxine...2-aminopyridine	1.01	0.97	0.98
	Adenine…thymine HB	0.95	0.90	0.89
	Methane dimer	1.16	1.11	1.14
	Ethene dimer	0.83	0.80	0.82
	Benzene…methane	0.64	0.62	0.70
	Benzene dimer stack	0.53	0.51	0.56
	Pyrazine dimer	0.61	0.58	0.51
	Uracil dimer stack	0.58	0.56	0.56
	Indole…benzene stack	0.25	0.24	0.43
	Adenine…thymine stack	0.72	0.68	0.61
	Ethene…ethyne	0.71	0.68	0.73
	Benzene…water	1.06	1.01	1.04
	Benzene…ammonia	0.70	0.67	0.79
	Benzene…cyanide	0.60	0.57	0.71
	Benzene dimer T	0.52	0.49	0.67
	Indole…benzene T	0.62	0.60	0.59
	Phenol dimer	0.95	0.91	0.99

a anharmonicity into account (1.8 -2.1 kcal mol -1 ).

30

The larger aromatic systems containing two strong hydrogen bonds show much smaller ∆ZPVE, falling to ca. 1 kcal mol -1 for 2-pyridoxine...2-

  The only exception to this conclusion is the C 2v T-shaped benzene dimer, which is a minimum with SCSN scaling but a saddle point with unscaled MP2. Zero-point vibrational energy corrections to binding energies have also been estimated, using raw data and two distinct scaling methods. In general, all three methods are in reasonable agreement, with corrections of between 20% and 50% of uncorrected binding energies.

	F o r F o r F o r F o r
	P P P P
	e e e e
	e r e r e r e r
	R R R R
	e e e e
	v i e v i e v i e v i e
	w w w w
	O n l O n l O n l O n l
	y y y y

Table 1

 1 SCF interaction energies and partitioned DF-LMP2 intermolecular energies (kcal mol-1 )

		SCF	Dispersion	Ionic	Dispersion-
					Exchange
	Ammonia dimer	-1.27	-1.23	-0.77	-0.01
	Water dimer	-3.72	-1.06	-1.23	+0.01
	Formic acid dimer	-15.49	-4.20	-5.63	+0.01
	Formamide dimer	-12.32	-3.73	-4.09	+0.08
	Uracil dimer HB	-9.31	-4.82	-4.79	+0.07
	2-Pyridoxine...2-aminopyridine	-5.81	-5.52	-4.76	-0.03
	Adenine…Thymine HB	-6.40	-5.68	-4.79	+0.01
	Methane dimer	+0.35	-0.73	-0.14	0.00
	Ethene dimer	+0.79	-1.88	-0.60	+0.01
	Benzene…methane	+1.60	-2.19	-0.86	+0.02
	Benzene dimer stack	+6.33	-7.22	-2.78	+0.02
	Pyrazine dimer	+7.60	-7.77	-3.11	+0.03
	Uracil dimer stack	+5.20	-9.15	-3.13	+0.01
	Indole…benzene stack	+8.62	-10.60	-4.46	+0.11
	Adenine…Thymine stack	+5.72	-14.05	-5.42	+0.05
	Ethene…ethyne	-0.52	-1.01	-0.58	0.01
	Benzene…water	-0.35	-2.01	-1.05	+0.01
	Benzene…ammonia	+0.74	-2.10	-0.95	+0.01
	Benzene…cyanide	-1.55	-2.75	-1.63	+0.04
	Benzene dimer T	+2.54	-3.84	-1.57	+0.05
	Indole…benzene T	+2.00	-5.13	-2.55	+0.05

Table 2

 2 Partitioned DF-LCCSD interaction energies for selected complexes (kcal mol-1 )

	Name	Dispersion	Ionic	Dispersion-
				Exchange
	Ammonia dimer	-1.29	-0.79	-0.01
	Water dimer	-1.20	-1.24	-0.03
	Methane dimer	-0.72	-0.16	0.00
	Ethene dimer	-1.78	-0.58	0.00
	Benzene…methane	-2.03	-0.82	+0.01
	Ethene…ethyne	-0.95	-0.54	0.00
	Benzene…water	-1.90	-0.99	-0.01
	Benzene…ammonia	-1.96	-0.85	-0.01
	This analysis was performed at the literature geometries, which were obtained with either MP2 or
	CCSD(T) optimisation with triple-or quadruple-ζ basis sets and including counterpoise correction.

Table 3

 3 RMS deviations from literature geometry (Å)

	Name	DF-LMP2	SCSN
	Ammonia dimer a	0.025	0.027
	Water dimer a	0.030	0.037
	Formic acid dimer a	0.017	0.016
	Formamide dimer a	0.025	0.028
	Uracil dimer HB	0.008	0.026
	2-Pyridoxine...2-aminopyridine	0.009	0.034
	Adenine…Thymine HB	0.011	0.028

b

For literature geometries optmised with CCSD(T) only.

Table 4

 4

	F o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y

is a step towards establishing accurate ZPVE corrections for systems of this size.

Table 4

 4 Harmonic zero point energy corrections for binding (kcal mol -1 )

	Name	DF-LMP2	Scaled a	SCSN
	Ammonia dimer	1.70	1.62	1.66
	Water dimer	2.37	2.27	2.33
	Formic acid dimer	2.08	1.99	2.06
	Formamide dimer	2.64	2.53	2.52
	Uracil dimer HB	1.21	1.16	1.15
	2-Pyridoxine...2-aminopyridine	1.01	0.97	0.98
	Adenine…thymine HB	0.95	0.90	0.89
	Methane dimer	1.16	1.11	1.14
	Ethene dimer	0.83	0.80	0.82
	Benzene…methane	0.64	0.62	0.70
	Benzene dimer stack	0.53	0.51	0.56
	Pyrazine dimer	0.61	0.58	0.51
	Uracil dimer stack	0.58	0.56	0.56
	Indole…benzene stack	0.25	0.24	0.43
	Adenine…thymine stack	0.72	0.68	0.61
	Ethene…ethyne	0.71	0.68	0.73
	Benzene…water	1.06	1.01	1.04
	Benzene…ammonia	0.70	0.67	0.79
	Benzene…cyanide	0.60	0.57	0.71
	Benzene dimer T	0.52	0.49	0.67
	Indole…benzene T	0.62	0.60	0.59
	Phenol dimer	0.95	0.91	0.99

  estimated, using raw data and two distinct scaling methods. In general, all three methods are in reasonable agreement, with corrections of between 20% and 50% of uncorrected binding energies.

	Page 29 of 31	Molecular Physics
	F o r F o r
	P P
	e e
	e r e r
		R R
		e e
		v i e v i e
		w w
		O n l O n l
		y y

URL: http://mc.manuscriptcentral.com/tandf/tmph

URL: http://mc.manuscriptcentral.com/tandf/tmph Molecular Physics

Acknowledgements

This work was funded by the UK EPSRC, grant ref. EP/C013328/1. Acknowledgements This work was funded by the UK EPSRC, grant ref. EP/C013328/1.


D Coffey

J L Diez-Ferrer

E C Corredor

J I Arnaudas

M Ciria
email: ciria@unizar.es

Manifold domain structure of double films with perpendicular magnetic anisotropy

Keywords: numbers: 75.60.Ch, 75.70.-i, 75.70.Ak, 75.70.Kw

We present epitaxial structures made of twin nickel blocks with perpendicular magnetic anisotropy separated by a copper layer which, for some values of this interleaving layer, show domain structures with four levels of contrast in magnetic force microscopy images. This manifold domain structure implies that the magnetization in the Ni blocks, besides the parallel orientation, undergoes a non-collinear configuration with respect to each other. To explain this result we consider a magnetoelastic domain structure with M in the plane that can elude the clamping done by the substrate with an average strain of -42 • 10 -6 (≈ 70% of the bulk value). Thus, the out-of-plane anisotropy is balanced and a biquadratic exchange coupling can stabilizes non-collinear domain configurations between the Ni blocks.

Introduction

The development of nanostructures including layers with perpendicular magnetic anisotropy has allowed keeping the pace of increasing the areal density in recording media [1]. These kind of layers aim to improve the performance of spintronic devices like planar nanowires [2] and nanopillar spin valves [3] since reducing the critical current for spin-transfer switching in these kind of systems appears to be more feasible than in layers with in-plane anisotropy. Patterned media for magnetic recording improves areal density by reducing the bit boundary noise of granular films, and, possibly, by means of multiple storage states [4]. Using this strategy 2 n available states can be achieved, where n is the number of magnetic blocks, by using layers with different values of the effective perpendicular magnetic anisotropy constant K ef f .

Most of these studies have been done in structures where the magnetic layers with perpendicular anisotropy are made of Co/Pt or Co/Pd blocks, and interleaving layers suitable for large antiferromagnetic coupling, Ru [5,6,7,8] or NiO [9,10,11]. The balance between interlayer interactions, namely magnetostatic and exchange coupling, dictates the domain configuration formed of these metastable states and tunning the oscillating exchange coupling can satisfactory stabilize remanent states [12]. A further step is the use of blocks with competing magnetic anisotropies (easy axis perpendicular to the surface and in the film plane) to obtain artificial magnetic structures with noncollinear and not orthogonal magnetic configurations [13].

A common factor of these structures is that a tiny difference in the energy, tuned by choosing structural parameters, can be favored by a small magnetic field, resulting in a remarkable change of the magnetic configuration, therefore high order contributions to the total magnetic energy may be important to explain the magnetic state.

Here we present magnetic force microscopy images in structures in which a copper block separates two twin nickel layers with out-of-plane magnetization due to a magnetoelastic (ME) effect. We observe that the number of levels of the MFM signal changes with the thickness of the copper block, and images with two, tree and up to four stable states are reported. Whereas the images with two and three levels are explained as result of dipolar and lineal exchange interaction that keep the magnetization in each Ni block perpendicular to the film plane but parallel or anti-parallel to each other, the fourfold contrast suggests the presence of domain structures with in-plane components of M in the Ni blocks. The stability of these four states is explained as a result of the existence of a biquadratic exchange coupling and the formation of magnetoelastic domains that elude the clamping done by the buffer layer. Thus, the metastable deformation of the nickel lattice in each domain is able to balance the perpendicular anisotropy and the biquadratic exchange interaction stabilizes the non collinear structure.

Experiment

Epitaxial Cu (5nm)/Ni (3nm)/Cu (100 nm) and Cu (5nm)/Ni (3nm)/Cu (t Cu )/Ni (3nm)/Cu (100 nm) films with t Cu = 3.5 nm, 4 nm, 4.5 nm and 6 nm were grown on Si (001) wafers at room temperature by electron-beam evaporation in a chamber with a base pressure below 2 x 10 -10 Torr, using a procedure reported elsewhere [14]. Hysteresis loops and domain images of the samples were taken at room temperature using a vibrating sample magnetometer (VSM) and a magnetic force microscope (MFM), respectively. Plane scans at constant average height from the sample surface are used to obtain the MFM images presented here. The short distance forces that produce the topographic profile can be minimized by adjusting the sample-tip distance, while the long distance magnetic force can still be detected in this configuration, greatly reducing the distortion of the domain structure that may happen in the more usual retrace mode.

Results and Discussion

Figure 1 displays hysteresis loops for double films with the magnetic field H applied parallel (thin continuous line) and perpendicular to the film plane (thick continuous line). The remanent magnetization M r is close to one for all the M-H loops with H perpendicular to the plane as expected for a loop taken along an easy direction; meanwhile the in-plane loops display no significant values for M r . These findings indicate the presence of an effective magnetic anisotropy perpendicular to the film, being the value of K ef f about 90 kJ/m 3 for all the films presented in this work [14].

The sharp change in the value of M at the value of the coercive field H c indicates that the inversion of M proceeds with the propagation of domain walls previously nucleated.

For the structures with t Cu = 3.5 nm and 4.5 nm a plateau is observed during the inversion of the magnetization that, on the contrary, is absent for the structures with t Cu = 4 nm and 6 nm.

The absence of a plateau in the M-H loop for the structure with t Cu = 4 nm, between 3.5 and 4.5 nm, suggests that the features around H c can be due to the presence of a oscillating interaction such as the interlayer exchange coupling (both bilinear and biquadratic). Nevertheless, the minor M-H loop, see dashed line in figure 1(a), done after saturating the sample in a positive H and applying a negative field of -185 Oe, shows a marginal shift, ≈ 10 Oe, suggesting a small value of a bilinear exchange coupling, as could be expected for structures with interleaving copper layers larger than 3 nm [15]. The in-plane M-H loops does not show anomalies related to the plateaus observed in the hysteresis loops taken with H perpendicular to the plane [see figure 1(b)-(d)].

The MFM technique is useful in the study of the magnetism in double films with perpendicular magnetization because it has revealed parallel (P) and antiparallel (AP) orientation of M between blocks [16]. Figure 2 shows images taken in the virgin state on the structures with t Cu = 4, 4.5 and 6 nm, and a 3 nm thick film. The structures with t Cu = 4 [figure 2 We note that these domain structures are unstable, being the stray field from the tip enough to modify the domain structures observed in figure 2. For the samples with t Cu = 4 and 6 nm and for the Ni thin film, the magnetic tip erases completely the domain structure in such a way that no contrast is observed after a few scans. This situation is different for the structure with t Cu = 4.5 nm. Figure 3(a) and (b) illustrate how a two states domain structure is obtained after the extreme hues areas turn into areas with the two intermediate hues due to the magnetic field from the tip. Enlarging the scan area to image untouched regions shows again the extreme contrast lost in the repeatedly scanned area, see figure 3(c). Note that while the white or black regions tend to disappear during the scanning, areas and bubbles with intermediate colors remain stable: see how the white bubble in figure 3 3(d)]. The persistence of metastable states seems to be associated with the hysteresis loop [see figure 1], since similar images, showing four states, have been taken for the t Cu = 3.5 nm structure.

Domain configurations have been explained by means of a magnetic state where the magnetization vector of the top and bottom Ni blocks, namely m top and m bot , are consider homogeneous along the normal direction because the nickel thickness is smaller than the nickel exchange length( A/K ≈ 10 nm, with A = 10 -11 J/m and K = 90 kJ/m 3 ). For MFM images with three hues m top and m bot are always perpendicular to the plane but may be parallel or antiparallel to each other. For the latter configuration the up-down or down-up domains provides the same contrast in the MFM image, because the dipolar field coming from the two AP configurations is much smaller (1/100 for scan distances in the range of 50-100 nm) than that coming from double films with P orientation [16] or from structures with an odd number of AP blocks [7].

For structures with t Cu = 3.5 and 4.5 nm, assigning to m top and m bot an orientation perpendicular to the plane and parallel to each other will produce domains with extreme signal (± 1.5), while a non collinear distribution with a net perpendicular component could explain the intermediate signals (± 0.6), see sketchs in figure 4(a). Nevertheless, the actual orientation of m top and m bot can not be established, due to the dependence of the magnetic signal with the lateral dimension of the domain structure [17], smaller domains giving rise to a larger signal; and due to the screening of the perpendicular field of the bottom block by a non perpendicular upper magnetic block.

The AP configuration is the result of a competition between antiferromagnetic exchange coupling and the dipolar interaction [6,10,18] or the result of two uncoupled films with different coercive field. We discard the latter consideration because the structure with the larger value of t Cu inverts M almost completely at the coercive field value. Therefore, the observation of images with larger number of levels in double films suggests the presence of structures where m top and m bot no longer are perpendicular to the plane and therefore additional magnetic contributions in the energy balance.

Model

In order to explain the manifold contrast observed in the MFM images we consider the relevant contributions to the energy density e mag that determine the orientation of M in the Ni blocks. It is well established that including in e mag the magnetostatic energy e ms term, an out-of-plane magnetic anisotropy Ksin 2 θ (θ is the angle between M and the film normal) and the domain wall energy, the lowest energy domain structure consists in regions with AP magnetization pointing perpendicularly to the plane, as is observed in Cu/Ni/Cu thin films [19]. Including an interleaving block between the magnetic blocks has required the presence of a bilinear exchange term in the magnetic energy, -J 1 m top m bot , (J 1 is the bilinear exchange constant) to explain the contrast in MFM images taken in Co/(Pt,Pd)/Ru based multilayers. In these films with AP orientation of m top and m bot , the observed contrast has been attributed to the domain wall shape [7]. A biquadratic term, -J 2 (m top m bot ) 2 , where J 2 is the biquadratic exchange coupling constant [20] favors a 90 o angle φ between m top and m bot . Thus competing interactions require, to find the equilibrium configuration, minimizing the total energy per block:

e mag /t Ni = e ms /t Ni + Ksin 2 θ/t Ni -J 1 m top m bot -J 2 (m top m bot ) 2 (1)
An estimation of the value of J 2 required to obtain angles of φ=90 o can be obtained assuming first that the biquadratic exchange contribution is able to balance the magnetic anisotropy, secondly that e ms is the value for a single domain (≈ (1/2)µ 0 M 2 cos 2 θ), and finally that J 1 is small. This gives J 2 ≈ -0.14 mJ/m 2 , a very large value considering the range of thicknesses of copper for which the presence of noncollinear structure is proposed, t Cu ≈ 3.5 -4.5 nm, and that typical values for |J 2 | ranges between 0.01 -0.1 mJ/m 2 , being the larger values found at thicknesses of the interlayer block no larger than 1 nm [21].

Noting that the image with up to four different levels of the MFM signal is observed as a metastable state since it is modified by the application of magnetic field, we consider the modification of the energy e mag arisen from a configuration of magnetic domains. Since the elastic and magnetoelastic energies are the most relevant contributions to understand the magnetic configuration in the Ni layers, we analyze the variation of these contributions if M would lie in the basal plane inducing a tetragonal magnetostrictive distortion in the magnetic domains, see figure 4(b), with the domains forming 90 o domain walls, figure 4(c). This mechanism is similar to the twin related domain formation that releases the elastic energy in tetragonally strained ferroelectric and ferroelastic epitaxial films [22], by the formation of structures of domains with the c axis of the tetragonal domain related to neighbor domains by a rotation of 90 degrees. This mechanism has been suggested to explain the domain configuration in magnetostrictive Co 50 Fe 50 /Co 80 B 20 films [23] and the magnetic order in Dy films (a crystal with hcp structure) grown on top of non-magnetic Y-Lu films [24]. These magnetoleastic domain structures gives rise to a in-plane strain distribution that was observed in Dy/Lu superlattices by the measuring of a three-fold splitting in the diffraction peaks due to domains with M pointing along the six in-plane easy directions [25].

In all the cases the domain pattering consists in volumes with alternating orientation of the magnetization or the polarization vector that due to the stress state associated with the M or P alternate compression and tension states inside each domain [see figure 4(c)]. Because the boundary for the in-plane lattice parameter imposed by the seed layer has to be fulfilled it has been suggested that the relaxation on the strained film can take place first at the domain boundaries that act as displacement dampers [22] and through the generation of anisotropic networks of linear densities of misfit dislocation [26].

The epitaxial strain in the film due to the energy balance between elastic and plastic contributions to the total energy is taken as reference. The thickness of the Cu seed layer, t seed , much larger than that value for Ni blocks, indicates that homogeneous stresses in the nickel layer have a negligible effect on the system because their weight on the total energy goes with a factor t Ni /(t Ni +t seed ). Thus, we evaluate the variation of the elastic density energy if the magnetostrictive strain is allowed to exist in the Ni film by the formation of a set of i domains, each one with a volume v i and with M in the plane along <10> directions, see figure 4(c). Inside each domain, to the epitaxial strain components ε 0 jk (with j, k = x, y, z defined along the <100> directions) the magnetostrictive deformations ε me jk are added, getting for the strain components ε jk = ε 0 jk + ε me jk . Taking the expression of the elastic energy for a crystal with cubic symmetry, inserting the strain components ε jk , keeping only contributions linear in the magnetostrictive strain, and considering that the nickel blocks are under a biaxial strain with the stress σ zz =0, we have that ε

0 xx = ε 0 xx = ε 0 ; ε 0 zz = -(2c 12 /c 11 )ε 0 ; ε 0 xy = ε 0 yz = ε 0 zx =0
, where c 11 and c 12 are elastic constants, the next equation is obtained for the energy released in each domain:

∆e i = (ε me xx + ε me yy )ε 0 (c 11 + c 12 ) -2 c 2 12 c 11 (2)
Note that the more usual magneto-elastic corrections effects, that go with (ε me jj ) 2 , result negligible compared with the terms of Equation ( 2) due to the strong internal deformation, ε 0 that for our films is ε 0 ≈ 2.1 10 -2 [14], while ε me ii ∼-6 10 -5 [27]. The relationship between magnetostrictive strain ε me , irreducible deformations of the cubic lattice, orientation of M and the measurement direction is,

ε me (α, β) = (3/2)λ 100 [α 2 x β 2 x + α 2 y β 2 y + α 2 z β 2 z -(1/3)] + 3λ 111 (α x α y β x β y + α y α z β y β z + α z α x β z β x ) ( 3 )
where the α i correspond to the cosines of M, the β i to the cosines of the measurement direction and λ 100 and λ 111 are the magnetostriction cubic coefficients that stand, respectively, for the tetragonal and rhombohedral distortions of the cubic lattice.

If M is along the <10> in-plane direction, we obtain, respectively, for the strain components parallel and transversal to M ε me 10 = λ 100 and ε me 01 = -λ 100 /2, getting for ∆e i =(1/2)ε 0 Cλ 100,i , where C includes the elastic constants combination in the braket of Equation ( 2). The total variation of the density energy is evaluated as a sum over the i domains with fraction volume v i in which the film is divided: ∆e = i v i ∆e i . Defining the average value of λ 100 as λ 100 = i (1/2)λ 100,i v i , the expression ∆e=(1/2)ε 0 C λ 100 is obtained.

For the Ni/Cu system ∆e release energy (∆e <0), since ε 0 ≈0.021 [14], and λ 100 is negative [27]. If ∆e were able to keep one of the magnetic layers with the domains in the film plane, it should overcome K ef f ≈ 90 10 3 J/m 3 . The value of λ 100 in the film in order to release enough elastic energy to match K ef f is ≈ -42 10 -6 (with c 11 = 2.5 10 11 J/m 3 and c 12 = 1.6 10 11 J/m 3 [28]), a value clearly smaller than the bulk λ 100 (-60 10 -6 [27]) and three orders of magnitude smaller than 0 . This calculation indicates that this configuration would balance the out-of-plane magnetic anisotropy because it is compatible with the reported values for λ 100 . Thus, beside domain configuration with m top and m bot perpendicular to the plane, either parallel or antiparallel to each other that can explain the MFM images for the structures with t Cu = 6 and 4nm, respectively, more sophisticated structures of domains with canted orientation of the magnetization with in-plane and out of plane components of M would be also possible if the coupling between the Ni blocks favors a non co-lineal orientation of m top and m bot .

We stress that when M is uniform the copper lattice clamps the nickel lattice to a certain value due to the energy balance between elastic and plastic contributions to the film total energy. Only during a process of inversion of the M the proposed domain configuration could take place. Thus, in a scenario where the perpendicular anisotropy is compensated by the formation of magnetoelastic domains, second order interactions favoring non-collinear states are required to stabilize the multilevel domain patterns proposed in figure 4(a). The biquadratic exchange interaction that favors a orientation of 90 degrees between m top and m bot , even being small in absolute value, could play the capital role of stabilizing non-collinear structures and be responsible for the plateaux observed in the MH loops since additional Zeman energy, supplied by increasing the magnetic field, is needed to unlock the non-collinear state. In order to estimate the strength of J 2 we consider that the energy associated to the area A in the M-H loops due to the presence of plateaux [see figure 1(c)] can be assigned to the biquadratic term, and A = -J 2 (m top m bot ) 2 /t Ni , getting J 2 ≈-0.005 mJ/m 2 , a reasonable value for this kind of interaction.

Summary

Magnetic structures with four well defined values have been observed by magnetic force microscopy in twin nickel blocks with an interleaving copper layer. This structure is explained as the result of a domain configuration that includes magnetoelastic domains with in-plane components and biquadratic exchange coupling. We note that engineering layers with metastable domain configurations based on the strain state of the magnetic film makes this kind of materials suitable to be used in conjunction with piezoelectric materials that modify their dimensions by means of an electric field and, therefore, the strain state and the ME anisotropy energy in the magnetic film.

Figure 1 .

 1 Figure 1. Detail of the M-H loops with H perpendicular to the plane (thick lines) for structures Ni(3 nm)/Cu(t Cu )/Ni(3 nm) with t Cu (a) 3.5 nm, (b) 4 nm, (c) 4.5 nm and (d) 6 nm. In panel (a) a minor loop is also shown (dashed line). M-H loops with H along the film plane are also shown (thin line) for the remaining structures. The grey area in panel (c) is used to estimate the strength of the biquadratic exchange contribution (see text).

  Figure1displays hysteresis loops for double films with the magnetic field H applied parallel (thin continuous line) and perpendicular to the film plane (thick continuous line). The remanent magnetization M r is close to one for all the M-H loops with H perpendicular to the plane as expected for a loop taken along an easy direction; meanwhile the in-plane loops display no significant values for M r . These findings indicate the presence of an effective magnetic anisotropy perpendicular to the film, being the value of K ef f about 90 kJ/m 3 for all the films presented in this work[14]. The sharp change in the value of M at the value of the coercive field H c indicates that the inversion of M proceeds with the propagation of domain walls previously nucleated. For the structures with t Cu = 3.5 nm and 4.5 nm a plateau is observed during the inversion of the magnetization that, on the contrary, is absent for the structures with t Cu = 4 nm and 6 nm.The absence of a plateau in the M-H loop for the structure with t Cu = 4 nm, between 3.5 and 4.5 nm, suggests that the features around H c can be due to the presence of a oscillating interaction such as the interlayer exchange coupling (both bilinear and biquadratic). Nevertheless, the minor M-H loop, see dashed line in figure1(a), done after saturating the sample in a positive H and applying a negative field of -185 Oe, shows a marginal shift, ≈ 10 Oe, suggesting a small value of a bilinear exchange coupling, as could be expected for structures with interleaving copper layers larger than 3 nm[15]. The in-plane M-H loops does not show anomalies related to the plateaus observed in the hysteresis loops taken with H perpendicular to the plane [see figure1(b)-(d)].The MFM technique is useful in the study of the magnetism in double films with perpendicular magnetization because it has revealed parallel (P) and antiparallel (AP) orientation of M between blocks[16]. Figure2shows images taken in the virgin state on the structures with t Cu = 4, 4.5 and 6 nm, and a 3 nm thick film. The structures with t Cu = 4 [figure2(a)] and 4.5 nm [figure 2(b)] show three and four levels in the contrast as shown the profile [figure 2(e)] taken along the white line in figure 2(a) and the histograms shown in figure 2(f). Notice that by taking the histograms on the areas marked in figure 2(b) up to four well defined peaks can be identified in the region delimited by a straight line by only two peaks if the selected area is the dashed rectangle. For the remaining structures [figure 2(c) and (d)] two hues are observed in the images, although the observation of two well defined peaks in the histograms taken on the MFM images depend on the signal to noise ratio [figure 2(g) and (h)].We note that these domain structures are unstable, being the stray field from the tip enough to modify the domain structures observed in figure2. For the samples with t Cu = 4 and 6 nm and for the Ni thin film, the magnetic tip erases completely the domain structure in such a way that no contrast is observed after a few scans. This situation is different for the structure with t Cu = 4.5 nm. Figure3(a) and (b) illustrate how a two states domain structure is obtained after the extreme hues areas turn into areas with the two intermediate hues due to the magnetic field from the tip.

Figure 2 .

 2 Figure 2. Magnetic force microscopy images taken on Ni(3)/Cu(t Cu )/Ni(3) double films with t Cu (a) 4 nm; (b) 4.5 nm (c) 6 nm and (d) a 3 nm thick Ni film, (e)profiles taken on image (a). (f-h) Histograms taken on the MFM images except for panel (f) where the continuous and dashed histograms correspond, respectively to the areas marked with continuous and dashed lines on image (b).

  (a) has disappeared in figure 3(b) (white ellipse) and the inverted U domain in figure 3(b) and 3(c) initially included the two extreme colors [see figure 3(a)]. Profiles taken along scan lines [see figure 3(d), continuous line] and histogram [figure 2(f), continuous line] show up to four levels or peaks if they include virgin areas and only two peaks if they are done in the scanned area [see figure 3(d), dashed line]. The domain structure achieved after the perturbation carried out by the tip has a signal that is in between the extremal values observed in the virgin state. The signal has maxima and minima at about 1.5 and -1.5 and the intermediate signal values are around 0.6 and -0.6 [see arrows in figure

Figure 3 .

 3 Figure 3. (a)-(c)MFM images taken on Ni(3)/Cu(4.5)/Ni(3) double film on similar areas. (d) Profiles taken on image (c)

Figure 4 .

 4 Figure 4. (a) Sketch of domain configurations that can provide 4 levels in the MFM signal. The parallel domains provide the extremal values while the domains forming 90 degrees provide two intermediate levels. The configuration on the left panel minimizes the dipolar coupling between the nickel blocks. (b) Sketch of the tetragonal distortion of a cubic material and the resulting lattices parameters. The values of c and a 1 are particularized for a irreducible tetragonal magnetostrictive deformation with value λ 100 .(c) Stress distribution for a film with a 90 o domain wall, the thick arrows stand for the magnetization or electric polarization vector. The small arrows indicate the in-plane stress component at a section of the film within each domain for a material with negative magnetostriction.

Acknowledgments

This work has been supported by Spanish MICINN (Grants No. MAT2006-07094 and MAT2009-10040) and DGA (Grants E81 and PI049/08) and Fondo Social Europeo.