- Email: [email protected]
Proton transfer in ice
Volume 210, number 1,2,3
CHEMICAL PHYSICS LETTERS
23 July 1993
Proton transfer in ice Changyol Lee ’ Departmenr
qf Phvsics.
Harvard
Universit.v. Carnhridge,
MA 02138,
IJSA
and David Vanderbilt Department
of Physics and Astronomy,
Rutgrrs
Llnirersity,
Piscataway.
NJ 08855-0849,
USA
Received I8 May 1993
Born-Oppenheimer energy surfaces for proton transfer in Hz0 cubic ice are calculated using gradient-corrected density-functional theory. Recombination of a pair of H&I* and OH- ionic defects is found to occur without barrier for both first-neighbor and second-neighbor pairs. The energy barriers for proton transfer in the configurations (H20-H-H20)+ and (OH-H-OH)are calculated to be only 0.069 and 0.058 eV, respectively. Relaxation of the atomic coordinates is found to be critical in determining the magnitude of the energy barriers. The relaxed geometries are not consistent with a solitonic description of the proton motion
The energetics of proton transfer along hydrogen bonds in water and other hydrogen-bonded substances is obviously of fundamental importance for understanding a wide variety of chemical reactions. However, direct experimental information about the energy surfaces is difficult to obtain, because the energy involved during the proton transfer is small compared with the hydrogen bond formation energy. Several ab initio theoretical studies [l-4] on the Born-Oppenheimer energy surfaces of the transferring proton have been carried out, using either molecular-orbital methods [l-3] or valence-bond methods [4]. In most of these calculations, the energy surfaces of the proton transfer were calculated for a dimer or chain of water molecules wirh all oxygen and all but one proton positions frozen, and the effects of relaxation were studied by comparing the energy surfaces for different oxygen separations. A controversial point has been whether the migration of the proton is cooperative (solitonic) or not (diffusive) [ $6 1. ’ Present address: Laboratory of Atomic and Solid State Physics, Cornell University, Ithaca, NY 14853-2501, USA. Elsevier Science Publishers B.V.
In this Letter, we present realistic calculations of energy surfaces for proton transfer in a realistic environment, as calculated from ab initio calculations using gradient-corrected density-functional theory [ 7-91. The total energies are calculated in a supercell geometry as a function of a reaction coordinate which describes the progress of the transferring proton along an O-O “bond”. The atomic coordinates are fully relaxed, subject only to the constraint that this generalized reaction coordinate have a specified value. We find that the relaxation is crucial to the determination of accurate potential energy surfaces, contrary to Scheiner [ 21, Relaxation is found to reduce the magnitude of the energy barriers. and Ieads to smaller barriers than those reported in earlier work [ l-41. Moreover, our calculations of the relaxed geometries do not support a solitonic picture of the proton motion. The details of our method of calculation are as follows. Density-functional total energies and forces are calculated via a steepest-descent version of the CarParrinello method [ IO]. The electronic wavefunctions are expanded in a plane-wave basis with a kinetic-energy cutoff of only 20 Ry; this is achieved by 279
Volume 210, number 1,2,3
CHEMICAL PHYSICS LETTERS
making use of Vanderbilt ultrasoft pseudopotentials [ 111. Local-density approximation [ 71 expressions for the exchange-correlation energy and potential in the pseudopotential generation and in the total-energy calculation are augmented with the gradient corrections by Perdew [ 81 and Becke [ 91. These corrections were previously found to be necessary for a proper description of the hydrogen bonding [ 12,131. When relaxation is allowed, an efficient algorithm is used in which the wavefunctions and atomic coordinates are updated simultaneously. Further details of the calculations and parameters of the pseudopotentials can be found in ref. [ 141. The details of the Car-Parrinello method with Vanderbilt ultrasoft pseudopotentials are discussed in ref.
[151In this work, we consider a simple-cubic supercell containing 8 HZ0 molecules in the structure of cubic ice. We have chosen cubic ice instead of hexagonal ice for reasons of technical simplicity. However, cubic and hexagonal ice have such a similar nearneighbor environment that we expect the results to apply equally well to both cases. While the geometry is not realistic for the case of liquid water, it is much less unrealistic than the geometries studied in previous theories [ l-41. Thus, the results may also have implications for aqueous chemistry. The lattice constant of 12.0 au for the eight-molecule unit cell is taken from experiment at - 130°C [ 161. Only the T point is used for the Brillouin-zone integration. The oxygens constitute a diamond structure and the hydrogens are randomly located within the constraints of Pauling’s ice rules [ 171. We use a proton configuration in which four of the Hz0 molecules have their dipole moment along f and the other four have their dipole moments along 3, -1, i, and -z^. This is the “most disordered” proton configuration consistent with the periodicity of an eightmolecule supercell. We first obtain the relaxed configuration of this defect-free reference structure. A nearest-neighbor defect pair is generated by displacing one of the 16 protons along the O-O bond on which it sits, so that it becomes part of an H,OC ionic defect and leaves an OH- behind. A second-neighbor defect pair is then made by following with a displacement of one of the other protons of the H30+ unit along its O-O bond, thereby shifting the H,O’ defect to the next molecule. In either case, all other 280
23 July 1993
atoms are initially held fixed; later, relaxations are allowed subject to a constraint which fixes the relative proximity of the proton to its two closest oxygen neighbors. We consider first the case of the nearest-neighbor defect pair, formed as above. The “unrelaxed” BornOppenheimer energy surface is obtained simply by calculating the total energy of the unit cell as the transferring proton moves along the O-O bond, keeping all other coordinates fixed. The “relaxed” energy surface, on the other hand, is obtained from a series of calculations in which a generalized coordinate is constrained to have a given value, while all other coordinates are relaxed via steepest-descent minimization of energy using the calculated forces as a guide. The generalized coordinate Q is designed to reflect the progress of the transferring proton along the O-O bond, and is defined in our case as Q= -do/d,
ifd,
Q=do/d2 if d2 cd,, , Q= (d, -dz)/
(d, +d, - 2d,) otherwise.
(1)
Here do is the nominal nearest-neighbor OH distance, d, = IRH, - Roa 1, and d2 = 1RH, - Roa 1, where H,, Od, and 0, are the transferring proton, the proton donor oxygen, and the proton acceptor oxygen, respectively. In this way, Q is continuously defined and has a value I or - 1 when H, is at the nominal distance from acceptor or donor oxygen, respectively. In order to keep the system constrained to a chosen value of Q during the steepest-descent minimization, we introduce a constraint function g({RJ)=O, where (R,) are the positions of the atoms. From eq. ( I), we have g=Qd, +do for d,
(2) where G, is the constraint force, Gj= ag/tlR,. When the atom coordinates are updated using this modificd force, the new coordinates automatically continue to satisfy the constraint. The energy surfaces from unrelaxed and relaxed geometries as a function of the generalized coordinate Q are shown in fig. la. Both curves have one
(a)
Table I Structure of the relaxed geometry of a pair of H,O+-OH- ionic defects at first (upper panel) and second (lower panel) nearestneighbor separation as a function of the generalized coordinate Q. Odl 0,. H, are donor, acceptor and the transferring proton, respecGvely, and Hs are other protons associated with donor or acceptor. All lengths are in au
I
‘i.d@D 0.2
I
:
\
:
X”
\ l-LzL.l I
I
I
-1.0
I.0
I@)
I
0.1
I
\ \ b ‘,
:
:
-O.-J
1
I
-1.0
1.0
-1.0
1.0
Q
0,-O,
0*-H,
0,-H,
0*-H
0,-H
-1.156 -0.922 -0.46 I 0.000 0.461 0.922 I.156
5.28 5.15 4.89 4.78 4.74 4.81 5.12
1.66 1.97 2.21 2.39 2.58 2.85 3.45
3.62 3.18 2.69 2.39 2.17 I .96 1.66
1.93 1.92 I.91 I .90 1.89 1.89 1.89
1.92, 1.92 1.93, 1.93 1.94, I.95 1.96, 1.97 1.99, 2.00 2.02,2.04 2.10, 2.13
-1.156 -0.922 -0.46 I 0.000 0.461 0.922 I.156
5.20 5.09 4.86 4.77 4.8 I 5.01 5.19
I .66 1.97 2.20 2.38 2.63 3.04 3.53
3.53 3.12 2.67 2.38 2.18 1.97 1.66
3.02, I.95 2.96, 1.94 2.84, 1.92 2.57, 1.92 2.20, I .93 2.03, 1.92 1.97, I.91
1.92, 1.92 1.93, 1.93 1.94, 1.94 1.97, 1.97 2.02, 2.02 2.10, 2.1 I 2.12, 2.18
ir:_,.
(e)
0.1
\ \
-. J
b, I
-1.0
23 July 1993
CHEMICAL PHYSICS LETTERS
Volume 2 IO, number 1,2,3
I
I
1.0
Fig. 1. The Born-Oppenheimer energy surfaces for proton transfer as a function of the generalized coordinate Q when the atoms are fixed (dashed) or relaxed (solid). Horizontal axes are Q and vertical axes are energy in eV. Zero of the energy is arbitrary. (a) Proton transfer from H20-Hz0 to H,O+-OH-. (b) Proton transfer from H,O-H,O+-OHto H,O+-H20-OH-. (c) Proton transfer in ( H20-H-H*O)+. (d) Proton transfer in (OHH-OH)-. (e) Proton transfer for ( H20-H-HZO)+ in vacuum. ( f ) Proton transfer for (OH-H-OH ) - in vacuum.
minimum, showing that the pair of positive and negative defects is unstable toward recombination, regardless of whether relaxation of atomic coordinates is allowed. Relaxed geometries are listed at table 1. The 0,-O, distance ranges from 4.74 to 5.28 au, and takes its minimum when the transferring proton is nearly midway between the oxygen!% As expected, as the transferring proton approaches (recedes), the other 0,-H or 0,-H bondlengths increase (decrease). However, the 0,-H bond distance only changes from 1.92 to 2.13 au, and the Od-H bond distances only change from 1.89 to 1.93 au, as the proton transfer takes place. This indicates that the proton movement is not cooperative enough to sup-
port a solitonic picture for the mechanism, since each H remains uniquely associated with its parent O,, or 0,.
Next, we further isolate the pair of H,O+ and OHionic defects to be second neighbors, and consider a proton transfer from the H,O+ (Q= 1) to the central Hz0 (Q= - 1) in such a way that the final configuration is that of a nearest-neighbor defect pair. Potential energy curves are shown in fig. lb and relaxed geometries are listed at table 1. As we move the proton away from the H30+ ion, the transferring proton is flanked by two neutral Hz0 molecules, i.e. we have OH--H20,-H: -H,O,. However, the electrostatic potential on H: due to the OH- ion introduces a strong asymmetry in the potential, so that the transferring proton prefers being around the H20d molecule, and the latter wants to release one proton back to OH- when relaxation is allowed. This can be seen by noting that as Q varies from 1.156 to - 1.156, the distance from Od to the released H changes from 1.97 to 3.02 au, which is beyond the midpoint of the O-O bond. Finally, the other 0,-H bond distances only change from 1.92 to 2.18 au during the proton transfer. Again, this is too small to support the solitonic picture. Having studied the recombination of oppositely 281
Volume 2 10, number I ,2,3
CHEMICAL PHYSICS LETTERS
charged defect pairs, we now turn to an investigation of the migration of isolated charge ionic defects by proton transfer. That is, we would like to study an isolated ( H20-H-HZO)+ or (OH-H-OH)-. Unfortunately, it is topologically impossible to arrange for a periodic supercell containing just one such ionic defect. Moreover, we have already seen that a supercell containing a pair of oppositely charged defects is unstable toward recombination without barrier, even for the optimal separation consistent with our eight-molecule unit cell, because of the strong electrostatic asymmetry. Therefore, we introduce the following technical artifice: we replace OH-- ( H20H-H,O)+ by HF-(H20-H-H20)+, and H30f(OH-H-OH)by NH,-(OH-H-OH)-. The resulting structures are shown in fig. 2. The neutral HF and NH3 species are introduced in order to terminate the network of hydrogen bonds in a way consistent with the Pauling ice rules, without introducing a second charge defect. In this way, the asymmetry seen by the transferring proton in ( H,O-H-H20 )’ or (OH-H-OH)is greatly reduced, because the large electrostatic contribution is eliminated. Nevertheless, some asymmetry remains. Moreover, we must keep in mind that the effect of the F or N impurity is not fully known, and that the unit cell now acquires a net charge.
(a)
Potential energy surfaces calculated in this way for the cases of (H,O-H-H,O)+ and (OH-H-OH)are shown in figs. I c and Id. Assuming that the perturbation due to the presence of F or N provides a linear gradient of potential energy in the region of interest, we calculate the common tangents of the two minima and subtract the resulting linear contribution from the tilted potential energy curve. From the potential energy curves which have been “symmctrized” in this way, we calculate the energy barrier for the proton transfer. For the unrelaxed geometry, the energy barriers in the cases of positive and negative defects are 0.104 and 0.142 eV, respectively. When relaxatidns are included, they become 0.069 and 0.058 eV, respectively. The relaxed geometries are listed in table 2. (One of the Od-H distances is longer, i.e. from 1.92 to 2.24 au, because the H in question is next to F.) The values for the energy barriers from previous work and the present work are summarized in table 3. Since we allow the ions to relax as the proton transfers, it is understandable that we have smaller energy barriers than those obtained in previous work, which assumed rigid structures. We find the energy barrier in ( H20-H-H20)+ is a little higher than that in (OH-H-OH ) - when relaxation is allowed. Some of the earlier work [ 1,4] suggested that the barrier Table 2 Structure of the relaxed geometry in the HJO+ (upper panel) and OH- (lower panel) defect as a function of the generalized coordinate Q. O,, O,, I-l,are donor, acceptor and the transferring proton, respectively, and Hs are other protons associated with donor or acceptor. All lengths are in au 0,-H,
0,-H,
0,-H
0,-H
5.13 -1.262 - 1.063 4.93 4.78 -0.623 0.000 4.76 0.623 4.81 1.063 4.96 1.262 5.13
I .66 I .98 2.21 2.38 2.60 2.99 3.47
3.41 2.96 2.58 2.38 2.21 1.98 1.66
I .94, 2.24 1.93,2.12 1.91,2.03 1.91, 1.99 1.90, 1.96 1.89, 1.94 1.89, 1.92
1.92, I .92 1.93, 1.93 1.95, I.95 1.97, 1.97 2.00, 2.00 2.05, 2.04 2.10, 2.09
5.21 -1.262 - 1.063 5.01 4.80 -0.623 0.000 4.75 0.623 4.77 1.063 4.94 1.262 5.16
1.66 1.98 2.21 2.38 2.51 2.97 3.50
3.55 3.04 2.60 2.38 2.21 1.98 I.66
1.93 I.91 1.90
1.89 1.89 I .90 1.90 I.91 1.93 1.95
Q
Fig. 2. (a) Structure of HF-(H20-H-H20)+ unit; black, grey, and white atoms are hydrogen, oxygen, and fluorine, respectively. (b) Same for NH>-(OH-H-OH)-, except the white atom 1snitrogen.
282
23 July I993
Od-0.
1.89 1.89 1.89 I .88
CHEMICAL PHYSICS LETTERS
Volume 210, number 1,2,3
Table 3 Energy barriers in positive and negative defects from previous work. Barrier in positive defect (eV)
Barrier in negative defect
4.71 5.22 5.71
0.19 0.98 2.20
0.62 1.90 3.50
4.72 4.91 5.22
0.08
0.38
0.31
0.71
0.65
[41
5.46
1.09
141 (21”
5.67 4.82
1.53 0.06
1.06 1.44 1.81
5.20 5.57
0.33 0.73
5.18 5.18
0.32 0.33
0.31 0.33
5.20 5.20
0.27 0.104
0.23 0.142
relaxed
0.069
0.058
Ref.
(11 (11 [II [41 141 141
i21"' [21" [31b' 131"
present (vacuum ) present (ice) present (ice)
O,,-0, bond length (au)
'IRigid staggered ( H&-H-H20)
(ev)
+ structure. b, SCF. ‘I MP3.
for the negative defect should be much higher than that in positive defect, because in the former case the proton moves between incomplete OH- radicals, while in the latter it moves between complete Hz0 molecules. More recent calculations by Scheiner [ 31 showed, however, that the energy barriers for the two cases are similar. From our calculations, we find that it is not important to the size of the barrier whether or not the transferring proton moves between complete molecules. In fact, the energy barrier of the positive defect is found to be larger than that of the negative defect, because the atomic relaxations are more important for the negative defect. However, this result is not necessarily inconsistent with a higher mobility of the positive defects, since the required relaxations (including motions of heavy oxygen atoms) would be expected to inhibit the mobility of the negative defects. We also study the energy surface for proton transfer in (H*O-H-H20)+ and (OH-H-OH)in vacuum, for comparison to the results of earlier work. Such configurations were considered in earlier work largely because the methods used were not practical for larger clusters. Potential energy curves are shown in figs. le and If. The energy barriers in (H20-H-
23 July 1993
HzO) + and (OH-H-OH) - from the unrelaxed geometry with Od-0, separation 5.20 au are 0.27 and 0.23 eV, respectively. These values are only 20°&30% smaller than those obtained in the most recent of the previous work [2,3]. It is notable that the hydrogens, not having other neighboring oxygens, do not relax appreciably; the Od-H or 0,-H bondlengths in the positive defect remain 1.88 au, and the Od-H or 0,-H bondlengths in the negative defect change bnly from 1.88 to 1.9 1 au, as Q changes from - 1.134 to 1,134. Most importantly, however, we find that when relaxation is allowed, the energy of the system is minimized when the transferring proton is at the center of the Od-0, bond. Thus, relaxation completely eliminates the energy barrier. Therefore, we believe it is not very meaningful to study the proton transfer for such a pair of molecules in vacuum. In conclusion, it is shown that the pair of H30+ and OH- ionic defects do not have a metastable state against recombination for first- and second-nearest neighbor separation. The transferring proton in (H*O-H-H20)+ and (OH-H-OH)is subject to an energy barrier of only 0.069 and 0.058 eV, respectively, when the ions are allowed to relax along with the motion of the transferring proton. The motion of protons is not cooperative enough to support a solitonic mechanism for the migration of defects. The calculations on ( H20-H-H20) ’ and (OH-HOH) - in vacuum are not physically meaningful, when relaxation is allowed. The authors acknowledge Professor Victor F. Petrenko for helpful discussions. This work is supported by NSF Grant No. DMR-91-15342. Supercomputer time was provided by the National Center for Supercomputing Applications under grant DMR920003N.
References [ I] M. Weissmann and N.V. Cohan, J. Chem. Phys. 43 ( 1965) 124. [2]S.Scheiner,J.Am.Chem.Soc.103(1981)315. [ 31 S. Scheiner, Accounts Chem. Res. 18 ( 1985) 174. [4] M. Hasegawa, K. Daiyasu and S. Yomosa, J. Phys. Sot. Japan 27 ( 1969) 999. [S] A. Godzik, Chem. Phys. Letters 17 I (1990) 217. [ 6 ] ES. Nylund and G.P. Tsironis, Phys. Rev. Letters 66 ( 1991) 1886.
283
Volume 2 IO, number 1,2,3
CHEMICAL PHYSICS LETTERS
[7] W. Kohn and L.J. Sham, Phys. Rev. 140 (1965) Al 133. [S] J.P. Perdew, Phys. Rev. B 33 (1986) 8822. [9] A.D. Becke, Phys. Rev. A 38 (1988) 3098. [IO] R. Car and M. Parrinello, Phys. Rev. Letters 55 (1985) 2471. [ 1I ] D. Vanderbilt, Phys. Rev. B 4 I ( 1990) 7892. [ 121C. Lee. D. Vanderbilt, K. Laasonen, R. Car and M. Parrinello, Phys. Rev. Letters 69 ( I992 ) 462. [ 131 K. Laasonen, M. Parrinello, R. Car, C. Lee and D. Vanderbilt, Chem. Phys. Letters 207 ( 1993) 208.
284
23 July 1993
[ 141C. Lee, D. Vanderbilt, K. Laasonen, R. Car, and M. Parrinello, Phys. Rev. B 47 (1993) 4863. [ 15 ] K. Laasonen, A. Pasquarello, R. Car, C. Lee and D. Vanderbilt, Phys. Rev. B 47 (1993) 10142. [ 161 M. Blackman and N.D. Lisgarten, Advan. Phys. 7 (1958) 189. [ 171 L. Pauling, Proc. Nat. Acad. Sci. US 14 ( 1928) 359.
CHEMICAL PHYSICS LETTERS
23 July 1993
Proton transfer in ice Changyol Lee ’ Departmenr
qf Phvsics.
Harvard
Universit.v. Carnhridge,
MA 02138,
IJSA
and David Vanderbilt Department
of Physics and Astronomy,
Rutgrrs
Llnirersity,
Piscataway.
NJ 08855-0849,
USA
Received I8 May 1993
Born-Oppenheimer energy surfaces for proton transfer in Hz0 cubic ice are calculated using gradient-corrected density-functional theory. Recombination of a pair of H&I* and OH- ionic defects is found to occur without barrier for both first-neighbor and second-neighbor pairs. The energy barriers for proton transfer in the configurations (H20-H-H20)+ and (OH-H-OH)are calculated to be only 0.069 and 0.058 eV, respectively. Relaxation of the atomic coordinates is found to be critical in determining the magnitude of the energy barriers. The relaxed geometries are not consistent with a solitonic description of the proton motion
The energetics of proton transfer along hydrogen bonds in water and other hydrogen-bonded substances is obviously of fundamental importance for understanding a wide variety of chemical reactions. However, direct experimental information about the energy surfaces is difficult to obtain, because the energy involved during the proton transfer is small compared with the hydrogen bond formation energy. Several ab initio theoretical studies [l-4] on the Born-Oppenheimer energy surfaces of the transferring proton have been carried out, using either molecular-orbital methods [l-3] or valence-bond methods [4]. In most of these calculations, the energy surfaces of the proton transfer were calculated for a dimer or chain of water molecules wirh all oxygen and all but one proton positions frozen, and the effects of relaxation were studied by comparing the energy surfaces for different oxygen separations. A controversial point has been whether the migration of the proton is cooperative (solitonic) or not (diffusive) [ $6 1. ’ Present address: Laboratory of Atomic and Solid State Physics, Cornell University, Ithaca, NY 14853-2501, USA. Elsevier Science Publishers B.V.
In this Letter, we present realistic calculations of energy surfaces for proton transfer in a realistic environment, as calculated from ab initio calculations using gradient-corrected density-functional theory [ 7-91. The total energies are calculated in a supercell geometry as a function of a reaction coordinate which describes the progress of the transferring proton along an O-O “bond”. The atomic coordinates are fully relaxed, subject only to the constraint that this generalized reaction coordinate have a specified value. We find that the relaxation is crucial to the determination of accurate potential energy surfaces, contrary to Scheiner [ 21, Relaxation is found to reduce the magnitude of the energy barriers. and Ieads to smaller barriers than those reported in earlier work [ l-41. Moreover, our calculations of the relaxed geometries do not support a solitonic picture of the proton motion. The details of our method of calculation are as follows. Density-functional total energies and forces are calculated via a steepest-descent version of the CarParrinello method [ IO]. The electronic wavefunctions are expanded in a plane-wave basis with a kinetic-energy cutoff of only 20 Ry; this is achieved by 279
Volume 210, number 1,2,3
CHEMICAL PHYSICS LETTERS
making use of Vanderbilt ultrasoft pseudopotentials [ 111. Local-density approximation [ 71 expressions for the exchange-correlation energy and potential in the pseudopotential generation and in the total-energy calculation are augmented with the gradient corrections by Perdew [ 81 and Becke [ 91. These corrections were previously found to be necessary for a proper description of the hydrogen bonding [ 12,131. When relaxation is allowed, an efficient algorithm is used in which the wavefunctions and atomic coordinates are updated simultaneously. Further details of the calculations and parameters of the pseudopotentials can be found in ref. [ 141. The details of the Car-Parrinello method with Vanderbilt ultrasoft pseudopotentials are discussed in ref.
[151In this work, we consider a simple-cubic supercell containing 8 HZ0 molecules in the structure of cubic ice. We have chosen cubic ice instead of hexagonal ice for reasons of technical simplicity. However, cubic and hexagonal ice have such a similar nearneighbor environment that we expect the results to apply equally well to both cases. While the geometry is not realistic for the case of liquid water, it is much less unrealistic than the geometries studied in previous theories [ l-41. Thus, the results may also have implications for aqueous chemistry. The lattice constant of 12.0 au for the eight-molecule unit cell is taken from experiment at - 130°C [ 161. Only the T point is used for the Brillouin-zone integration. The oxygens constitute a diamond structure and the hydrogens are randomly located within the constraints of Pauling’s ice rules [ 171. We use a proton configuration in which four of the Hz0 molecules have their dipole moment along f and the other four have their dipole moments along 3, -1, i, and -z^. This is the “most disordered” proton configuration consistent with the periodicity of an eightmolecule supercell. We first obtain the relaxed configuration of this defect-free reference structure. A nearest-neighbor defect pair is generated by displacing one of the 16 protons along the O-O bond on which it sits, so that it becomes part of an H,OC ionic defect and leaves an OH- behind. A second-neighbor defect pair is then made by following with a displacement of one of the other protons of the H30+ unit along its O-O bond, thereby shifting the H,O’ defect to the next molecule. In either case, all other 280
23 July 1993
atoms are initially held fixed; later, relaxations are allowed subject to a constraint which fixes the relative proximity of the proton to its two closest oxygen neighbors. We consider first the case of the nearest-neighbor defect pair, formed as above. The “unrelaxed” BornOppenheimer energy surface is obtained simply by calculating the total energy of the unit cell as the transferring proton moves along the O-O bond, keeping all other coordinates fixed. The “relaxed” energy surface, on the other hand, is obtained from a series of calculations in which a generalized coordinate is constrained to have a given value, while all other coordinates are relaxed via steepest-descent minimization of energy using the calculated forces as a guide. The generalized coordinate Q is designed to reflect the progress of the transferring proton along the O-O bond, and is defined in our case as Q= -do/d,
ifd,
Q=do/d2 if d2 cd,, , Q= (d, -dz)/
(d, +d, - 2d,) otherwise.
(1)
Here do is the nominal nearest-neighbor OH distance, d, = IRH, - Roa 1, and d2 = 1RH, - Roa 1, where H,, Od, and 0, are the transferring proton, the proton donor oxygen, and the proton acceptor oxygen, respectively. In this way, Q is continuously defined and has a value I or - 1 when H, is at the nominal distance from acceptor or donor oxygen, respectively. In order to keep the system constrained to a chosen value of Q during the steepest-descent minimization, we introduce a constraint function g({RJ)=O, where (R,) are the positions of the atoms. From eq. ( I), we have g=Qd, +do for d,
(2) where G, is the constraint force, Gj= ag/tlR,. When the atom coordinates are updated using this modificd force, the new coordinates automatically continue to satisfy the constraint. The energy surfaces from unrelaxed and relaxed geometries as a function of the generalized coordinate Q are shown in fig. la. Both curves have one
(a)
Table I Structure of the relaxed geometry of a pair of H,O+-OH- ionic defects at first (upper panel) and second (lower panel) nearestneighbor separation as a function of the generalized coordinate Q. Odl 0,. H, are donor, acceptor and the transferring proton, respecGvely, and Hs are other protons associated with donor or acceptor. All lengths are in au
I
‘i.d@D 0.2
I
:
\
:
X”
\ l-LzL.l I
I
I
-1.0
I.0
I@)
I
0.1
I
\ \ b ‘,
:
:
-O.-J
1
I
-1.0
1.0
-1.0
1.0
Q
0,-O,
0*-H,
0,-H,
0*-H
0,-H
-1.156 -0.922 -0.46 I 0.000 0.461 0.922 I.156
5.28 5.15 4.89 4.78 4.74 4.81 5.12
1.66 1.97 2.21 2.39 2.58 2.85 3.45
3.62 3.18 2.69 2.39 2.17 I .96 1.66
1.93 1.92 I.91 I .90 1.89 1.89 1.89
1.92, 1.92 1.93, 1.93 1.94, I.95 1.96, 1.97 1.99, 2.00 2.02,2.04 2.10, 2.13
-1.156 -0.922 -0.46 I 0.000 0.461 0.922 I.156
5.20 5.09 4.86 4.77 4.8 I 5.01 5.19
I .66 1.97 2.20 2.38 2.63 3.04 3.53
3.53 3.12 2.67 2.38 2.18 1.97 1.66
3.02, I.95 2.96, 1.94 2.84, 1.92 2.57, 1.92 2.20, I .93 2.03, 1.92 1.97, I.91
1.92, 1.92 1.93, 1.93 1.94, 1.94 1.97, 1.97 2.02, 2.02 2.10, 2.1 I 2.12, 2.18
ir:_,.
(e)
0.1
\ \
-. J
b, I
-1.0
23 July 1993
CHEMICAL PHYSICS LETTERS
Volume 2 IO, number 1,2,3
I
I
1.0
Fig. 1. The Born-Oppenheimer energy surfaces for proton transfer as a function of the generalized coordinate Q when the atoms are fixed (dashed) or relaxed (solid). Horizontal axes are Q and vertical axes are energy in eV. Zero of the energy is arbitrary. (a) Proton transfer from H20-Hz0 to H,O+-OH-. (b) Proton transfer from H,O-H,O+-OHto H,O+-H20-OH-. (c) Proton transfer in ( H20-H-H*O)+. (d) Proton transfer in (OHH-OH)-. (e) Proton transfer for ( H20-H-HZO)+ in vacuum. ( f ) Proton transfer for (OH-H-OH ) - in vacuum.
minimum, showing that the pair of positive and negative defects is unstable toward recombination, regardless of whether relaxation of atomic coordinates is allowed. Relaxed geometries are listed at table 1. The 0,-O, distance ranges from 4.74 to 5.28 au, and takes its minimum when the transferring proton is nearly midway between the oxygen!% As expected, as the transferring proton approaches (recedes), the other 0,-H or 0,-H bondlengths increase (decrease). However, the 0,-H bond distance only changes from 1.92 to 2.13 au, and the Od-H bond distances only change from 1.89 to 1.93 au, as the proton transfer takes place. This indicates that the proton movement is not cooperative enough to sup-
port a solitonic picture for the mechanism, since each H remains uniquely associated with its parent O,, or 0,.
Next, we further isolate the pair of H,O+ and OHionic defects to be second neighbors, and consider a proton transfer from the H,O+ (Q= 1) to the central Hz0 (Q= - 1) in such a way that the final configuration is that of a nearest-neighbor defect pair. Potential energy curves are shown in fig. lb and relaxed geometries are listed at table 1. As we move the proton away from the H30+ ion, the transferring proton is flanked by two neutral Hz0 molecules, i.e. we have OH--H20,-H: -H,O,. However, the electrostatic potential on H: due to the OH- ion introduces a strong asymmetry in the potential, so that the transferring proton prefers being around the H20d molecule, and the latter wants to release one proton back to OH- when relaxation is allowed. This can be seen by noting that as Q varies from 1.156 to - 1.156, the distance from Od to the released H changes from 1.97 to 3.02 au, which is beyond the midpoint of the O-O bond. Finally, the other 0,-H bond distances only change from 1.92 to 2.18 au during the proton transfer. Again, this is too small to support the solitonic picture. Having studied the recombination of oppositely 281
Volume 2 10, number I ,2,3
CHEMICAL PHYSICS LETTERS
charged defect pairs, we now turn to an investigation of the migration of isolated charge ionic defects by proton transfer. That is, we would like to study an isolated ( H20-H-HZO)+ or (OH-H-OH)-. Unfortunately, it is topologically impossible to arrange for a periodic supercell containing just one such ionic defect. Moreover, we have already seen that a supercell containing a pair of oppositely charged defects is unstable toward recombination without barrier, even for the optimal separation consistent with our eight-molecule unit cell, because of the strong electrostatic asymmetry. Therefore, we introduce the following technical artifice: we replace OH-- ( H20H-H,O)+ by HF-(H20-H-H20)+, and H30f(OH-H-OH)by NH,-(OH-H-OH)-. The resulting structures are shown in fig. 2. The neutral HF and NH3 species are introduced in order to terminate the network of hydrogen bonds in a way consistent with the Pauling ice rules, without introducing a second charge defect. In this way, the asymmetry seen by the transferring proton in ( H,O-H-H20 )’ or (OH-H-OH)is greatly reduced, because the large electrostatic contribution is eliminated. Nevertheless, some asymmetry remains. Moreover, we must keep in mind that the effect of the F or N impurity is not fully known, and that the unit cell now acquires a net charge.
(a)
Potential energy surfaces calculated in this way for the cases of (H,O-H-H,O)+ and (OH-H-OH)are shown in figs. I c and Id. Assuming that the perturbation due to the presence of F or N provides a linear gradient of potential energy in the region of interest, we calculate the common tangents of the two minima and subtract the resulting linear contribution from the tilted potential energy curve. From the potential energy curves which have been “symmctrized” in this way, we calculate the energy barrier for the proton transfer. For the unrelaxed geometry, the energy barriers in the cases of positive and negative defects are 0.104 and 0.142 eV, respectively. When relaxatidns are included, they become 0.069 and 0.058 eV, respectively. The relaxed geometries are listed in table 2. (One of the Od-H distances is longer, i.e. from 1.92 to 2.24 au, because the H in question is next to F.) The values for the energy barriers from previous work and the present work are summarized in table 3. Since we allow the ions to relax as the proton transfers, it is understandable that we have smaller energy barriers than those obtained in previous work, which assumed rigid structures. We find the energy barrier in ( H20-H-H20)+ is a little higher than that in (OH-H-OH ) - when relaxation is allowed. Some of the earlier work [ 1,4] suggested that the barrier Table 2 Structure of the relaxed geometry in the HJO+ (upper panel) and OH- (lower panel) defect as a function of the generalized coordinate Q. O,, O,, I-l,are donor, acceptor and the transferring proton, respectively, and Hs are other protons associated with donor or acceptor. All lengths are in au 0,-H,
0,-H,
0,-H
0,-H
5.13 -1.262 - 1.063 4.93 4.78 -0.623 0.000 4.76 0.623 4.81 1.063 4.96 1.262 5.13
I .66 I .98 2.21 2.38 2.60 2.99 3.47
3.41 2.96 2.58 2.38 2.21 1.98 1.66
I .94, 2.24 1.93,2.12 1.91,2.03 1.91, 1.99 1.90, 1.96 1.89, 1.94 1.89, 1.92
1.92, I .92 1.93, 1.93 1.95, I.95 1.97, 1.97 2.00, 2.00 2.05, 2.04 2.10, 2.09
5.21 -1.262 - 1.063 5.01 4.80 -0.623 0.000 4.75 0.623 4.77 1.063 4.94 1.262 5.16
1.66 1.98 2.21 2.38 2.51 2.97 3.50
3.55 3.04 2.60 2.38 2.21 1.98 I.66
1.93 I.91 1.90
1.89 1.89 I .90 1.90 I.91 1.93 1.95
Q
Fig. 2. (a) Structure of HF-(H20-H-H20)+ unit; black, grey, and white atoms are hydrogen, oxygen, and fluorine, respectively. (b) Same for NH>-(OH-H-OH)-, except the white atom 1snitrogen.
282
23 July I993
Od-0.
1.89 1.89 1.89 I .88
CHEMICAL PHYSICS LETTERS
Volume 210, number 1,2,3
Table 3 Energy barriers in positive and negative defects from previous work. Barrier in positive defect (eV)
Barrier in negative defect
4.71 5.22 5.71
0.19 0.98 2.20
0.62 1.90 3.50
4.72 4.91 5.22
0.08
0.38
0.31
0.71
0.65
[41
5.46
1.09
141 (21”
5.67 4.82
1.53 0.06
1.06 1.44 1.81
5.20 5.57
0.33 0.73
5.18 5.18
0.32 0.33
0.31 0.33
5.20 5.20
0.27 0.104
0.23 0.142
relaxed
0.069
0.058
Ref.
(11 (11 [II [41 141 141
i21"' [21" [31b' 131"
present (vacuum ) present (ice) present (ice)
O,,-0, bond length (au)
'IRigid staggered ( H&-H-H20)
(ev)
+ structure. b, SCF. ‘I MP3.
for the negative defect should be much higher than that in positive defect, because in the former case the proton moves between incomplete OH- radicals, while in the latter it moves between complete Hz0 molecules. More recent calculations by Scheiner [ 31 showed, however, that the energy barriers for the two cases are similar. From our calculations, we find that it is not important to the size of the barrier whether or not the transferring proton moves between complete molecules. In fact, the energy barrier of the positive defect is found to be larger than that of the negative defect, because the atomic relaxations are more important for the negative defect. However, this result is not necessarily inconsistent with a higher mobility of the positive defects, since the required relaxations (including motions of heavy oxygen atoms) would be expected to inhibit the mobility of the negative defects. We also study the energy surface for proton transfer in (H*O-H-H20)+ and (OH-H-OH)in vacuum, for comparison to the results of earlier work. Such configurations were considered in earlier work largely because the methods used were not practical for larger clusters. Potential energy curves are shown in figs. le and If. The energy barriers in (H20-H-
23 July 1993
HzO) + and (OH-H-OH) - from the unrelaxed geometry with Od-0, separation 5.20 au are 0.27 and 0.23 eV, respectively. These values are only 20°&30% smaller than those obtained in the most recent of the previous work [2,3]. It is notable that the hydrogens, not having other neighboring oxygens, do not relax appreciably; the Od-H or 0,-H bondlengths in the positive defect remain 1.88 au, and the Od-H or 0,-H bondlengths in the negative defect change bnly from 1.88 to 1.9 1 au, as Q changes from - 1.134 to 1,134. Most importantly, however, we find that when relaxation is allowed, the energy of the system is minimized when the transferring proton is at the center of the Od-0, bond. Thus, relaxation completely eliminates the energy barrier. Therefore, we believe it is not very meaningful to study the proton transfer for such a pair of molecules in vacuum. In conclusion, it is shown that the pair of H30+ and OH- ionic defects do not have a metastable state against recombination for first- and second-nearest neighbor separation. The transferring proton in (H*O-H-H20)+ and (OH-H-OH)is subject to an energy barrier of only 0.069 and 0.058 eV, respectively, when the ions are allowed to relax along with the motion of the transferring proton. The motion of protons is not cooperative enough to support a solitonic mechanism for the migration of defects. The calculations on ( H20-H-H20) ’ and (OH-HOH) - in vacuum are not physically meaningful, when relaxation is allowed. The authors acknowledge Professor Victor F. Petrenko for helpful discussions. This work is supported by NSF Grant No. DMR-91-15342. Supercomputer time was provided by the National Center for Supercomputing Applications under grant DMR920003N.
References [ I] M. Weissmann and N.V. Cohan, J. Chem. Phys. 43 ( 1965) 124. [2]S.Scheiner,J.Am.Chem.Soc.103(1981)315. [ 31 S. Scheiner, Accounts Chem. Res. 18 ( 1985) 174. [4] M. Hasegawa, K. Daiyasu and S. Yomosa, J. Phys. Sot. Japan 27 ( 1969) 999. [S] A. Godzik, Chem. Phys. Letters 17 I (1990) 217. [ 6 ] ES. Nylund and G.P. Tsironis, Phys. Rev. Letters 66 ( 1991) 1886.
283
Volume 2 IO, number 1,2,3
CHEMICAL PHYSICS LETTERS
[7] W. Kohn and L.J. Sham, Phys. Rev. 140 (1965) Al 133. [S] J.P. Perdew, Phys. Rev. B 33 (1986) 8822. [9] A.D. Becke, Phys. Rev. A 38 (1988) 3098. [IO] R. Car and M. Parrinello, Phys. Rev. Letters 55 (1985) 2471. [ 1I ] D. Vanderbilt, Phys. Rev. B 4 I ( 1990) 7892. [ 121C. Lee. D. Vanderbilt, K. Laasonen, R. Car and M. Parrinello, Phys. Rev. Letters 69 ( I992 ) 462. [ 131 K. Laasonen, M. Parrinello, R. Car, C. Lee and D. Vanderbilt, Chem. Phys. Letters 207 ( 1993) 208.
284
23 July 1993
[ 141C. Lee, D. Vanderbilt, K. Laasonen, R. Car, and M. Parrinello, Phys. Rev. B 47 (1993) 4863. [ 15 ] K. Laasonen, A. Pasquarello, R. Car, C. Lee and D. Vanderbilt, Phys. Rev. B 47 (1993) 10142. [ 161 M. Blackman and N.D. Lisgarten, Advan. Phys. 7 (1958) 189. [ 171 L. Pauling, Proc. Nat. Acad. Sci. US 14 ( 1928) 359.