Introduction
The theoretical investigation of intermolecular interactions, particularly those involving hydrogen bonds, is of fundamental importance to chemistry, biology, and solid state and surface physics. Because of its importance in many chemical and biological processes, the study of the properties of water in all of its phases has been an area of intensely active research for the past decade. The earlier theoretical studies of ice (Weissmann and Cohan, 1965) were confined to a model of the lattice in which there were only two or three interacting molecules. Often, only the electrons involved in the hydrogen bond were included in the calculations. However, since 1972, several authors (Crowe and Santry, 1973; Santry, 1973; Cotterill and others, 1973; Scott, 1975; Plummer and Stein, 1976) have shown the power of semi-empirical molecular-orbital calculations for studying larger systems of ice and ice-like lattices. These calculations concentrated on determinations of structural parameters, molecular electronic charge distributions, and stabilization energy for idealized infinite ice lattices or small ice-like aggregates. It was found that these semi-empirical quantum-mechanical calculations were in surprisingly close agreement with those properties which could be compared with experimental values. Approximate molecular-orbital calculations of the CNDO/INDO type (Pople and Beveridge, 1970) give good agreement (Plummer, 1972 and unpublished results) with the results of ab initio calculations of a small basis set of small clusters of water molecules, and even predict the non-additive effects of hydrogen bonding found near the Hartree-Fock limit of the extended basis calculation (Hankins and others, 1970), thus giving confidence that the use of these techniques is justified and does yield physically meaningful results.
The nature of the surface of ice itself has long been a subject of controversy (Fletcher, 1970, 1973; Jellinek, 1972; Nason and Fletcher, 1975). The only semi-quantitative theoretical treatment of the ice surface has been developed by Fletcher (1970, 1973) but this addressed the question of the presence or absence of a liquid-like layer on the surface and not that of the stability of the different planes of ice nor the energetics of the surface-monomer interaction. Since the interaction of the surface with an incoming or adsorbed water monomer is of major importance in understanding the processes of nucleation and crystal growth, this report will concentrate on these aspects of the ice problem. The monomer-surface interaction is modelled by using a small section of an ice surface and then applying semi-empirical quantum-mechanical calculations of the CNDO/INDO type to investigate the energetics of a gas-phase molecule approaching the surface. Preliminary investigations of the water-ice surface using a ten-molecule cluster to represent the ice surface have been reported (Plummer and Stein, 1976) and compared with the results from an empirical pair-potential model (Plummer and others, 1976). The latter model is discussed in more detail in Kiefer and Hale (1977). Additional comparisons between the quantum-mechanical calculations and the pair potential will be made in a future report (Plummer, in press).
In an attempt to simplify the system as much as possible while retaining the essential features of the "real" physical situation, an idealized section of the ice surface was represented by thirteen water molecules arranged in the lattice positions of ice Ih. Questions into the effects of surface relaxation will be dealt with in subsequent work in which molecular dynamics will be used. The problem of the nature of "real" ice surface is not specifically addressed, but these studies together with the molecular-dynamics studies can be very instructive and provide valuable data on a molecular level.
Potential-energy curves were calculated for a monomer approaching different sites on both the basal and prism planes. Three-dimensional electron-density plots of the ice surfaces were generated. The specific geometries and surface structures considered are discussed. Then, details of the calculations are discussed together with a brief description of the CNDO technique, and we discuss the results of these calculations. Conclusions and implications for future work are included in the last section.
When a small group of molecules are used to represent the ice surface, we adopt the assumptions that the oxygen sublattice of ice Ih extends to the surface without re-arrangement, this appears to be experimentally justified for sufficiently low temperatures, (Firment and Somorjai, 1975),* and that a molecule approaching the surface is influenced primarily by
first and second nearest neighbors of the approach site. This latter point is amplified and justified later.
The Model
For the calculation of the interaction energies of a water monomer with the basal plane, a group of thirteen water molecules was arranged on a two-level ice Ih lattice as illustrated in Figure I. The positions of the oxygens are indicated by large circles while the hydrogen positions for a specific basal surface are shown as small filled circles. The c-axis is normal to the plane and the oxygens represented by filled circles lie 0.9 Å (1 Å = 10-10 m) below the plane. The oxygen-oxygen distances were fixed at 2.76 Å. The intramolecular O-H separation was 1.0 Å, slightly greater than in the vapor (0.957 Å), and no distinction was made in O-H separation between "free" hydrogens and those participating in hydrogen bonds. The nearest-neighbor distance between molecules in the same plane (dashed lines) was 4.5 Å, the lattice constant a 0. All H-O-H angles were opened up to 109.50°, the tetrahedral angle. The area of this section of the surface is approximately 80 ÅZ. For the basal face, four distinct proton configurations have been examined: (I) a configuration in which all the molecules in the first layer (open circles) had protons pointing out of the plane parallel to the c-axis, (2) a configuration in which no protons pointed out of the plane, (3) a configuration in which the protons on the central site pointed away from the c-axis while the remaining protons alternated up and down (referred to as BA) and (4) a second alternating configuration of surface protons in which the remaining protons were arranged so as to give the maximum surface stability (Plummer, 1972 and unpublished results). The calculations reported here were for the BA surface as it appeared to be the best representation of an "average" basal surface of ice. Additionally, entropy considerations would make a geometry like BA, with more proton disorder than the optimal surface, physically more probable. It has also been suggested (Pruppacher and Pflaum, 1975) that such a surface would be more conducive to the nucleation of new planes of ice than would a more polar surface.
For the basal plane, the interaction energies for the monomer approaching four specific types of surface sites were calculated. The site labelled 1 in Figure I was considered as the primary bonding position and the potential surface was then calculated for a variety of rotational orientations of the approaching monomers. The remaining three "sites" were not epitaxial positions but were of interest in "mapping out" the energy contours of the surface. These non-epitaxial positions should play a role in surface diffusion and in any liquid-like layer on the surface. One such position was in the center of one of the six-membered rings. A second position was directly above a second-layer oxygen. In these two cases as well as for the primary site, the incoming monomer was oriented so as to donate one of its hydrogens for bond formation with the surface. Thus the dipole of the approaching monomer was oriented approximately 52° from the surface. The fourth surface site was located half-way between site 1 and a first-layer nearest-neighbor molecule. The monomer approaching this site was oriented so both hydrogens pointed toward the surface.
For evaluating the interaction energy of a monomer with a prism face of ice Ih, the molecules represented the surface were arranged as shown in Figure 2. The c-axis is in the plane of the molecules and the plane of the oxygens represented by the closed circles lies 1.38 Å below those represented by open circles. Only one configuration of surface protons was considered for the prism face. Again the primary prism site is designated as 1 (Fig. 2).
The only other prism site examined had the monomer approaching with both hydrogens available for bonding to the surface. It should be noted that, since the hydrogen bonds which extend from the surface are not perpendicular to the surface for the prism face, nearest-neighbor surface molecules are less suitable for bifurcated bonding with an approaching monomer.
In addition to the several neutral basal and prism surface geometries, ion-pair defects were introduced into the surface and the interaction energies of these determined with an approaching monomer. The ion-pair defects were produced by moving a proton from one end of the hydrogen bond to the other (see Fig. 3). The central site was always one molecule of the ion pair. In one such defect pair, the surface molecule (first layer) had an H3O+ configuration while its companion in the second layer was OH-. In a second configuration the surface molecule was negative OH- while the H3O+ was in the second layer. A third ion-pair configuration has the ion pair separated by a neutral water and both ions were in the first layer of molecules (see Fig. 4). Since such defects are present in concentrations of 1 in 106 (Hobbs, 1974, chapters 2 and 6) their effects on the interaction should be studied. If the effects are sufficiently small, the ion defects could be neglected in any statistical study of the ice surface- vapor- surface cluster system
The Calculations
For these studies of the monomer-ice surface interaction all valence electrons, both of the monomer and of the molecules used to represent the surface, were included. The atomic orbital basis for these calculations was of the conventional Slater valence-set type, that is, is orbitals on all the hydrogens and 2S, 2p x , 2p y orbitals centered on each oxygen. The problem then consisted of calculating the expectation value of the energy operator for a variety of monomer positions using product wave functions describing the 42 nuclei and 112 electrons. In order to make the problem tractable (and not too expensive) a number of the multicenter integrals were neglected (neglect of differential overlap) and the two electron repulsion integrals were parameterized from experimental atomic-spectra data. The specific technique used has been widely applied to molecular systems and systems of molecular aggregates. These approximate molecular-orbital theories are known by the acronyms CNDO, INDO, or NDDO depending on whether the neglect is "complete", "intermediate", or limited to "neglect of diatomic differential overlap" (Pople and Beveridge, 1970).
We have used both the CNDO and INDO methods and have found the qualitative predictions to be the same from both. The results reported here are CNDO/2 calculations (the 2 identifies the specific parametrization used for the integrals) which could be more readily compared with and calibrated against previous results for water systems (Hoyland and Kier, 1969; Dobosh, 1974; Scott, 1975).
The calculations were performed on an IBM 370/168 computer and required approximately 1.5 minutes of CPU time per monomer-surface geometry. The calculation of the electron 'density for a specific surface configuration required approximately 1.2 min of CPU time, with an additional 0.15 min to generate a three-dimensional plot.
The stability of a given surface geometry was determined by subtracting 13 times the energy of an individual surface monomer from the energy for the 13-molecule aggregate. The binding energy of the monomer E b with the surface was determined by subtracting the energy of the surface E s plus the energy of the monomer E n , from the total energy for the 14-molecule aggregate E t
The net atomic charge on any atom I a was equal to the core charge on the atom ?a minus the gross electron population on that atom P a
For all of the calculations, the geometry of the incoming monomer was fixed at an O-H separation of 1.0 Å and a H-O-H angle of 105°, and these were not allowed to change upon bond formation. Additionally, no relaxation of the surface was allowed, either for a change in surface geometry or upon interaction with the approaching monomer.
For each basal or prism site the interaction energies were calculated for 12 to 20 monomer surface distances surrounding the minimum separation. A least-squares fit to the resulting potential curves allowed the calculation of vibrational frequencies for a monomer-surface stretch. For the principal bonding sites on each face, the energy as a function of rotation about the monomer-surface bond was calculated and barriers to rotation estimated.
Results and Discussion
These calculations have used a small cluster of 13 water molecules to model the surface of ice Ih. The use of clusters of atoms to model the solid substrate in chemisorption investigations (Johnson and Messmer, 1974 ; Messmer, 1977) produced results in reasonable agreement with experiment. Since we were attempting to model the interaction of a water monomer with an infinite ice surface, we needed to determine whether the properties sought could be attributed to localized properties of the extended solid and whether the number of molecules used to represent the surface was sufficient to reproduce the short-range effects of the monomer-surface interaction.
It has been found in simulations of ice (Santry, 1973) that intermolecular electron exchange and delocalization contributions are quite short range, dominated by nearest-neighbor interactions and convergent by second-nearest-neighbor inclusion. We found the same trends for the monomer-surface interactions. For the basal plane, 94% of the interaction energy could be attributed to nearest-neighbor interactions with an additional 5.4% from second-nearest-neighbor interactions. For the prism face, convergence was not quite so rapid, and third-nearest-neighbor interactions were needed to provide greater than 99% of the interaction energy. With the cluster employed here the monomer interaction included fourthnearest neighbors for the principal bonding sites on each face and third-nearest neighbors elsewhere. Thus, it appears reasonable to assume that the energetics of the molecule close to the surface could be modelled using this small section of the surface. The principal limitation imposed by the use of the small cluster is that it cannot be expected faithfully to mimic relaxation of the extended crystal when defects are incorporated. Thus no relaxation was allowed in these structures, and our attention was directed to the perfect ice surface.
Monomer-basal-surface interactions
To visualize more easily how the surface would appear to an incoming water molecule, the electron densities at a number of points on the surface were calculated. Then the density in any plane could be displayed as a three-dimensional plot in which the relative magnitude of the density is represented by the height of the peaks in the plot. One such plot for the basal (BA) face is shown in Figure 5. This plot illustrates the electron density perpendicular to the c-axis in the plane of the upper layer of molecules (refer to Fig. I). The area of the plot is 180 Å2 and the distance between grid points is 0.25 Å. The viewer is oriented 30° above the plane and 45° clockwise from the orientation depicted in Figure 1. The slight asymmetries in the plot are due to the disorder in the proton configurations. As expected, the lone pair of the oxygen atoms produced the peak densities, and these together with the effect of molecular dipoles served to orient the incoming monomer. However, the magnitude of these effects was quite small beyond 4 or 5 Å. To illustrate the result that nearest-neigh b or interactions dominate and how short-range many of the surface features are, Figure 6 shows the electrondensity plot I Å above the surface.
The change in electron density with distance above the plane has been examined more closely in the vicinity of the approach sites. In Figure 7 the density is plotted against distance above the plane for points along the line connecting the primary bonding site and the site for bifurcated bonding. The center peak is located at the primary site; bifurcated bonding occurs midway between the center and right-most peaks, with the two hydrogens in the plane and pointing toward the peaks. Figure 8 illustrates the electron density for an approach trajectoryover a second-layer molecule. The approach site is at the center peak; the other points lie on a line parallel to that of Figure 7. The scale of peak heights has been increased in Figure 8 to show more detail.
The most stable configuration of basal (BA) and monomer had the monomer 2.55 Å above the primary site with the non-bonded monomer hydrogen oriented 180° from the surface molecule H-O-H bisector. The potential curve for this interaction is shown in Figure 9 and has an energy minimum of -33.3 kJ/mol (34.7 x 10-2 eV). The barrier to rotation about the monomer- surface bond was 2. 1 kJ/mol (2.17 x 10-2 eV) which is approximately kT. The stretching frequency for the bond, calculated via a least-squares fit to the potential curve, was 280 cm-1.
The second most stable monomer-surface geometry was the one in which the monomer was in a bifurcated position and had both hydrogens participating in bonds to the surface. The distance of the monomer oxygen from the surface at the energy minimum was 1.7 Å and the energy was -28.4 kJ/mol (0.30 eV). The potential curve for the monomer approaching this site together with that for the monomer over a second-layer molecule are shown in Figure 10. Figure 11 shows the potential barrier for rotation between the primary site with one hydrogen bond and the bifurcated site with two monomer-surface bonds. The barrier height was 11.3 kJ/mol (0.117 eV) above the primary site. Other orientations of the monomer dipole above the bifurcated site were investigated but all had weaker bonding to the surface.
The remaining two basal sites for which the monomer-surface interaction was calculated were not true bonding configurations since no monomer-surface hydrogen bond was formed. However these are both local attractive minima in the potential surface. The attraction interaction is approximately one-half that of the primary bonding site. The magnitude and exact location of these local minima are highly dependent on the orientation of both the surface protons and those of the incoming water monomer.
Monomer-prism-surface interactions
For the prism face, electron-density plots of the surface were again generated in order to assess in a better way the effects of the surface on an incoming monomer. The plot obtained for the prism plane is shown in Figure 12 where the density is plotted parallel to the c-axis and through the upper layer of oxygen atoms (see Fig. 2). The magnitude of the interaction of the monomer with the primary site is found to be slightly less than for the basal face (-33 kJ/mol) although the difference is at the limit of the accuracy of the calculations. However, bonding to the bifurcated site was substantially less for the prism face amounting to an energy of 16.3 kJ/mol. The barrier for rotation about the monomer surface bond at the primary site was about the same as for the basal face: 2.09 kJ/mol. The vibrational frequency calculated for elongation of this bond (266 em-1) is slightly less than for the basal plane.
Monomer-ion-pair-defect interactions
Ordinarily the Bernal-Fowler rules apply in ice, that is, each water molecule is surrounded by two near hydrogens and two hydrogens at approximately twice the distance of the near hydrogens. These conditions guarantee the neutrality of each molecule as well as assuring the presence of one and only one hydrogen per bond. However, if this were always the case it would be hard to explain the observed electrical properties of ice. Bjerrum (1951) suggested that certain types of point defects, which violated the Bernal-Fowler rules, could be present and that their migration could provide an explanation for many of the electrical properties of ice (Camplin and Glen, 1973; Glen and Paren, 1975). The types of point defects suggested were ion-pair defects and orientational defects. In the former, a proton is transferred from one water molecule to its neighbor producing an H30+ and an OH- pair. If the negative ion then captures a proton from an adjacent molecule before recombination occurs, the ionic defect has been somewhat stabilized and we have the first step in charge migration. In the case of the orientational defects, rotation of a single molecule can put two hydrogens on the same bond (a D-defect) and simultaneously create a bond having no hydrogens (an L-defect). In this report we will discuss only the ionic defects. The interactions of the monomer with a surface containing orientational defects as well as isolated ions will be reported later (Plummer, in press).
As described earlier, a study of ion-pair defects on the basal surface was performed. An upper bound for the energy required to produce the defects is given by
where E d is the defect energy, E s n is the energy of a neutral surface, and E s i is the energy associated with an ionic surface. This number will be an over-estimate, since it includes only the energy required for ionization and charge separation but no compensating relaxation of the surface. The values obtained ranged from 600 to 850 kJ/mol (6.3 to 8.8 eV), about the same magnitude as the proton affinity of water, 670 kJ/mol (7 eV).
The primary basal site for monomer interaction was always one of the ion pairs and was the only site of monomer approach included in these studies. In the first case, a proton was moved from a second-layer molecule producing an H30+ ion at the surface having all protons pointing into the surface. This orientation was not conducive to hydrogen-bond formation with an incoming monomer but nevertheless exhibited an attractive interaction in excess of 20.9 kJ/mol. When the proton was shifted from the primary site to a second-layer molecule producing an OH- primary site, a very strong bond with the incoming monomer resulted. The potential energy as a function of distance for these interactions is illustrated in Figure 13. The negative primary site also exhibited a higher barrier to rotation about the monomer-ion bonds than found with the neutral site.
A third ion-pair configuration was produced by shifting the proton from the second-layer H30+ to an adjacent upper-layer molecule (see Fig. 4). This had the effect of separating the ion pair and having both ion defects in the upper layer of molecules. The energy required to produce the separated defect was 210 kJ/mol (2.2 eV). When the ions were separated in this fashion, the strength of the interaction with an incoming monomer was essentially the same as for a monomer with a "free" OH- ion and greater than in the second type of configuration.
Conclusions
In this study we have modelled the interaction of a water monomer with an extended, perfect ice surface represented by a cluster of 13 water molecules whose oxygens were arranged at the lattice sites for hexagonal ice. Summaries of the results appear in Tables I and II. We have shown that the interaction energies are dominated by nearest-neighbor effects and hence this model appears quite reasonable.
The attraction energies were calculated for a variety of possible attachment sites on both prism and basal faces. It was found that epitaxial sites, those which allowed a hydrogen bond to be formed between the monomer and the surface, produced the largest interaction energies. These energies were approximately the same (-33.4 kJ/mol) for the two ice crystal surfaces.
The second-most-stable configuration has two "bonds" formed between the monomer and the basal surface. This bifurcated structure was separated from the most-stable configuration by an energy barrier of 6.3 kJ/mol. This orientation of the monomer is probably most frequently observed as a metastable or intermediate structure for a molecule in transition between two of the epitaxial bonding sites. This position can be regarded as a local minimum in the potential for surface diffusion. Its existence and the moderate barrier between it and a more stable epitaxial structure, makes reasonable the consideration of a "bipedal walk" process for the surface diffusion of water on the basal surface of ice. The situation is quite different on the prism face where a similar bifurcated geometry was found to be less stable than "nonbonding" configurations on the basal ice surface.
The other two approach sites considered for the basal surface were not true bonding configurations, but they proved to be attractive sites with interaction energies in the range of - (17-21) kJ/mol. The magnitude of this attraction will be highly dependent on the orientation of both the monomer and the surface hydrogens and will be investigated more extensively in future work. Since these are local minima in the surface potential, they can attract and physically adsorb a monomer if its incoming energy is not too great. The monomer could be quasi-bound for some time before being re-evaporated or diffusing across the surface to be incorporated into an epitaxial position. These sites can also provide added stabilization for the addition of a second molecule to a monomer already bound to a primary site thereby enhancing the initial step in cluster formation.
From these calculations we can estimate a "static" diffusion rate and residence time for a monomer on the basal surface. For these evaluations, we use the expressions of Burton and others (1951).
and
where χs is the diffusion time, α is the distance between nearest neighbors in the same plane, v is the stretching frequency, v' is the diffusion frequency, E A and E D are the average binding and diffusion energies, and τs is the residence time. For the average binding energy E A we used that associated with a primary site, 32 kJ/mol. The diffusion barrier was taken as to kJ/mol and v' and v were 450 cm-1 and 28o cm-1 respectively. For a temperature of 260 K, we find the mean path length to be approximately 270α, where α is 4.5 Å, the distance between nearest neighbors in the same plane. The residence time is estimated as 7 x 10-7 s. It should be noted that these are estimates of dynamic quantities from static calculations. Additionally, they should be regarded as an upper limit to these quantities as we have taken the bottom of the potential well as our value for E A so that zero-point vibrational energies have not been excluded. Probably vibrational excitation should be included as well. However, these refinements did not seem justified in the light of the other approximations inherent in Equations (4) and (5). The estimate of χs = 270α may be compared with the estimate of Burton and others of χs = 400α for the (111) face of an f c.c. crystal.
We also examined the effect on the monomer-surface interaction of incorporating ion-pair defects in the surface. As one might expect, we found the interaction energies to be greater, ranging from three to five times as strong as with the neutral site. The potential wells were also steeper, yielding higher stretching frequencies. The energies required to produce the defects compare well with the experimental proton affinity of water and also with calculated barriers to proton transfer (Weissmann and Cohan, 1965). Future studies will address the question of lattice relaxation and defect mobility. It would be expected that the presence of ion-defect pairs or isolated ions in the surface would significantly alter diffusion rates near those sites. Indeed, because of the strength of the interaction, the results of these calculations predict no diffusion from the ionic sites.
In summary, I emphasize the following points: First, while the strengths of the bondings between the basal and prism planes are roughly the same, the studies reported here suggest that the diffusion rates and residence times may be substantially different for the two crystal faces. Therefore, factors such as supersaturation and flux of monomers to the surface, as well as the temperature of the surface and vapor, will play a large role in determining the ice-crystal growth habit. Such factors should be studied in a dynamical fashion and such calculations are under way. Secondly, the ice surface should be viewed not as having isolated strong-bonding regions separated by repulsive regions but rather as consisting of attractive regions separated by low energy barriers. The non-bonding local minima provide the energy mechanism whereby a molecule approaching a surface does not have to attach initially to a primary bonding site but can be attracted and remain in contact with the surface long enough to either move to a more favorable bonding position on the surface or become attached to a molecule or cluster already bonded to the surface. This, in turn, again suggests that surface diffusion would be a primary mechanism for surface cluster formation and growth. Third, the presence of these non-bonding sites together with the essentially free rotation about monomer-surface bonds provide motions not present in the crystal, and hence may be the molecular origin of a liquid-like layer on the ice surface.
While these calculations are only a model for the real interaction, they represent a first attempt to model this interaction from a molecular point of view using quantum-mechanical calculations. The results reported do provide considerable information about the monomer - surface interaction and within the limits of the model allow us to draw preliminary conclusions about processes on the surface of ice.
Acknowledgements
The author would like to acknowledge the help of E. Stein and K. Whisnant for their computational assistance and work with the ten-molecule surfaces. She wishes to thank the Computer Facilities at the University of Missouri for much of the computer time needed for these studies. She also thanks Dr James L. Kassner, Jr for his continued encouragement of this work and Dr 0. R. Plummer, for many useful discussions and suggestions.
Discussion
W. B. KAMB: What you call a "bifurcated site" does not involve bifurcated hydrogen bonds as that term is conventionally used; you could perhaps better call it a "doubly-bonded site".
Also, the particular interaction between the H30+ ion and monomer that you described is one where no H-bond is formed. If this is definitely more stable than the H-bonded configuration, it might suggest that, in the ice crystal, H30+ ions and L-defects should be permanently associated.
P. L. M. PLUMMER: Perhaps it would be better to refer to the single-proton-down configuration as a "single donor" and the two-proton-down configuration as a "double donor" configuration, then an oxygen-down orientation could be a "single" or "double" acceptor depending on its attachment to the surface.
I cannot address your second question as to whether the H30+ and L-defects should be permanently associated in the bulk, but for the surface-monomer case, where the monomer is free to rotate so as to optimize the interaction when all the protons on H3O+ are directed into the solid, the L-defect configuration is more stable. It should be recalled, however, that the all-proton-down orientation is only one of four for H30+, the others are orientated for hydrogen-bond formation with an incoming monomer which would, of course, be a much stronger interaction.
J. W. GLEN: Whilst agreeing with Dr Kamb's second remark, I disagree with the first as I believe the "bifurcated" approach was to two molecules both in the top layer and thus tending to a four-oxygen ring. As the apparent failure of the monomer to give hydrogen bonding in one orientation of the H30+ ion on the surface is potentially important, as Kamb has pointed out, I am interested to know whether you tried OH- ions in both orientations (proton out of the surface as well as proton in the surface).
PLUMMER: No, I have not as yet, but that calculation will be done. However, for the OH-, either orientation is favourable for hydrogen-bond formation with the monomer.
J. E. BERTIE: Your calculations will certainly be used by those interested in the formation and growth of ice crystals from the vapour. It is particularly important, therefore, that an estimate of the effect of the approximations in your treatment be given. I note that you used the CNDO and INDO approximations to the ab initio molecular-orbital treatment. These are very approximate in principle and are nowadays only used in cases which cannot be treated by the full molecular-orbital method. Since the energies that you calculated are small (≲40 kJ/mol), it is particularly necessary to give an estimate of the effect of your approximations and their values. Can you give this estimate please?
PLUMMER: It is certainly true that the CNDO/INDO techniques are semi-empirical in nature and must be used with great care. The results are much more reliable for comparing different structures within the same types of systems where you are looking at energy differences, and are not nearly as good for either absolute total energies or for comparisons between systems (say water systems against organics). It is my feeling that the values for bonding energies are probably good to 10 or 15% and that differences greater than 4 kJ/mol are real. Those less than 1.2 kJ/mol may well not be. We have compared the CNDO with other semi-empirical methods, empirical potentials, and with ab initio calculations, when possible, and these are the ranges of accuracy which seem reasonable from these comparisons. The separations are surely about to% too small.