J. Phys. Chem. B 2009, 113, 8475–8480
8475
Aspects of the Proton Transfer in Liquid Phosphonic Acid Jan-Ole Joswig* and Gotthard Seifert Physikalische Chemie, Technische UniVersita¨t Dresden, 01062 Dresden, Germany ReceiVed: NoVember 7, 2008; ReVised Manuscript ReceiVed: April 30, 2009
Born-Oppenheimer molecular dynamics simulations have been performed to investigate the proton-transfer mechanism in liquid phosphonic acid at different temperatures. The transfers are observed if additional charge carriers are present. The liquid phosphonic acid has been characterized by pair-distribution functions, selfdiffusion coefficients, and hopping rates. Moreover, we find that the proton transfer is prepared by the following necessary local geometrical criteria: (1) an O-H-O angle close to 180° and (2) an O-O distance lower than 2.5 Å. However, we observe in many cases a Zundel-like species, in which an excess proton is shared, complexed, and stabilized by two adjacent phosphonic acid molecules. These unsuccessful attempts to transfer a proton fulfill the same local criteria as successful transfers, and thus, we conclude that nonlocal requirements may play an important role for the success of the transfer too. Introduction Phosphonic-acid-based polymers are currently intensively investigated as possible candidates for electrolyte membranes in proton-exchange membrane fuel cells (PEM-FCs). In a fuel cell, the membrane has to separate the gaseous compounds H2 and O2, to transport the protons in one direction, and to prevent transport of electrons in the opposite direction. Thus, high demands have to be made on materials used as proton-exchange membranes. Presently, the most commonly used proton conductor is Nafion, which is a sulfonated tetrafluoroethylene copolymer. It does not only fulfill the above-mentioned criteria but is additionally thermally and mechanically very stable. This makes it an ideal material for proton-exchange membranes. However, proton-exchange membranes may have certain disadvantages that lower the fuel-cell performance; the operating temperature of the fuel cell is limited to approximately 100 °C if the proton transport is water-based (as in Nafion). Moreover, methanol, if used as a hydrogen source, easily supports poisoning of platinum-based electrodes by carbon monoxide. To elaborate on new and different proton-conducting polymers for fuel cells, other compounds have been investigated, for example, imidazole-based,1,2 sulfonic,3 or phosphonicacid-based4-6 polymers. Especially, the latter are very promising. Therefore, we are currently focusing on the proton-transfer mechanisms in phosphonic-acid-based polymers and phosphonic acid (PA) itself. Although phosphonic acid and PA-based polymers are experimentally studied extensively for fuel cell applications, fundamental information about the proton-transfer mechanism is still missing. Therefore, we have chosen to study the protontransfer mechanism in liquid phosphonic acid first, in order to pass over to larger systems including PA-based polymers and the proton-transfer mechanism within these. Phosphonic acid has not drawn very much attention from the theoretical point of view so far. Some structural investigations,7,8 exist and the gas-phase acidity, protonation energies, transfer barriers, and geometries9-11 have been studied. A small number * To whom correspondence should
[email protected].
be
addressed.
E-mail:
of studies have put effort into aspects related to ours, for example, the investigations of the hydration of propyl phosphonic acid,12 the structure, and potential energy surfaces of methyl and heptyl phosphonic acid.13 Moreover, several studies investigated proton-transfer mechanisms or proton mobility in systems not related to phosphonic acid by means of molecular dynamics. Some of the studied systems are related to the electrolyte in fuel cells as well, for example, imidazole chains14 or liquid methanol15 and heptyl phosphonic acid.16 Recent studies employing molecular dynamics simulations have dealt with superprotonic solid-state conductors such as CsHSO4 and CsH2PO4.17,18 Both are especially interesting as electrolytes for higher operating temperatures. Another solid-state material for separator applications in fuel cells are perovskite-type oxides. Some experimental and theoretical studies on these were reviewed by Kreuer19 earlier. Finally, important work has been published on the proton transport and mobility in water. Here, the pioneering work of Tuckerman and co-workers20-22 has to be mentioned. These studies investigate the structure and dynamics of charge defects in terms of excess or missing protons in water using ab initio molecular dynamics simulations. For an excess proton, fluctuations between the so-called Zundel (H5O2+) and Eigen (H9O4+) ions are observed. This work is also the basis of the present study since we use a similar molecular dynamics setup as will be described below. In an earlier study,23 we have compared the structural data of phosphonic acid obtained with different levels of theory. The study included the bonding information, infrared frequencies, proton affinity, and finally, the energetics of the proton transfer in a protonated dimer. Thereupon, we have performed molecular dynamics simulations of liquid phosphonic acid at different temperatures, which we present in this study. For this purpose, different simulation setups have been used; the simulation box either contains phosphonic acid molecules only, or in order to introduce a charge carrier, one PA molecule is either protonated or deprotonated. In the following section, we will briefly revise our computational methods, afterward discuss our results, and finally, summarize and conclude.
10.1021/jp8098335 CCC: $40.75 2009 American Chemical Society Published on Web 05/29/2009
8476
J. Phys. Chem. B, Vol. 113, No. 25, 2009
Joswig and Seifert
Computational Methods Section The Density Functional Tight-Binding Method. A density functional tight-binding method (DFTB)24,25 has been used in the present study. The method is based on the density functional theory of Hohenberg and Kohn26 in the formulation of Kohn and Sham.27 The single-particle Kohn-Sham eigenfunctions ψi(r) are expanded in a set of localized atom-centered basis functions φm(r). These functions are determined by selfconsistent density functional calculations on the isolated atoms employing a large set of Slater-type basis functions. The total energy of the system of interest relative to the isolated atoms is written as
Etot )
∑ εi - ∑ εkj i
kjk
k
+
1 2
∑ Ukl(|Rk - Rl|)
(1)
k*l
where εi is the Kohn-Sham eigenvalue of the ith orbital of the system of interest, εkjk is the eigenvalue of the isolated atom (k being an atom index, jk an orbital index), and Ukl is a pair potential between atoms k and l, which is fitted to parameterfree density functional calculations of the several molecules28 employing a hybrid B3LYP exchange-correlation functional.29,30 All electrons except for the 3s and 3p electrons of phosphorus, the 2s and 2p electrons of oxygen, and the 1s electrons of hydrogen were treated within a frozen core approximation. Furthermore, a self-consistent charge (SCC) scheme31 corrects the results by adding Coulomb interactions of the atomic charges to the Hamiltonian, so that the total energy reads
Etot )
∑ cimcinHmn0 + 21 ∑ ∆qk∆qlγkl + Erep imn
(2)
kl
0 where Hmn are the Hamilton matrix elements, ∆q are atomic charge fluctuation calculated via the Mulliken charges, and the γ function includes the second-order contribution of the exchange-correlation energy for small distances, where it is important, and goes for 1/r for large interatomic distances representing the Coulomb interaction between two point charges. The repulsive term is defined as in eq 1. We have validated the approximations made in the DFTB approach in comparison to ab initio calculations in an earlier study on phosphonic acid.23 Moreover, the method has been successfully applied in several studies on proton mobility and proton-transfer reactions in different systems.14,19,28,32 Computational Details. Born-Oppenheimer molecular dynamics (MD) simulations have been performed using the DFTB method implemented in the deMon computer code.33 The cubic simulation box contained 32 PA molecules (224 atoms). Its size (a ) 13.97 Å) was chosen according to the density of phosphonic acid. Periodic boundary conditions were applied, and the simulations were restricted to the Γ-point. A CarParrinello treatment of the employed DFTB method34 is currently being tested for future application. A fcc (face-centered cubic) arrangement of the molecules has been chosen as a starting configuration, from which the molecules have been randomly displaced and oriented. The systems were equilibrated for 35 ps at different temperatures using the DFTB Hamiltonian. After the initial equilibration, the data of the 30 ps production runs within a microcanonical (NVE) ensemble were analyzed. The MD time step was set to 0.25 fs, and the trajectories were propagated in time using the velocity Verlet algorithm.
Figure 1. Pair-distribution functions of different pairs of atoms (P-P, P-O, P-H). The sharp peaks up to 2.5 Å result from the ordered molecular structure, whereas the PDF at larger distances results from different molecules in the liquid. Inset: Pair-distribution functions of all present pairs of atoms (containing at least one P atom) for two different temperatures.
Three setups have been used, one containing 32 PA molecules (224 atoms), one in which one PA molecule has been protonated additionally (225 atoms), and one setup in which a proton was abstracted from one PA molecule (223 atoms). In the following, these three setups will be referred to as the neutral, positive, and negative setup, respectively. For each setup, the MD simulations have been performed at six different temperatures spanning the range between melting and boiling point of phosphonic acid (350-470 K). Results Liquid Phosphonic Acid. To characterize the liquid phosphonic acid, we have calculated the pair-distribution functions (PDFs) from the MD trajectories. The PDF is defined as
g(r) )
V N2
〈∑ ∑ N
N
i
j*i
〉
δ(r - rij)
(3)
The normalization factor accounts for the volume V of the simulation box and the number of atoms N. The interatomic distances between atoms i and j are denoted as rij. In Figure 1, we display the PDF for the atomic pairs P-P, P-O, and P-H. Up to a distance of approximately 2.5 Å, the PDFs show a quite ordered structure resulting from the fixed molecular geometry of the PA molecule. The P-O curve (green) has two very sharp peaks, which result from the single- and double-bonded oxygen atoms. The P-H curve (blue) shows a sharp peak at 1.4 Å representing the P-H bond. The much broader peak centered around 2.3 Å gives the distance between the P atom and the hydroxyl H atoms. This second peak has a larger width at half-maximum, which shows that the hydroxyl H atoms move more strongly with respect to the P atom compared to the P-bonded hydrogen atoms. One reason is certainly that they are not directly bonded, so that a change in the P-O-H angle changes the P-H distance significantly. The P-H PDF (blue curve) does not decay to zero at 2.5 Å. This indicates a small amount of elongated P-H distances, either through elongated O-H bonds or through large P-O-H angles. Both may occur during a proton transfer, but from the PDF, we can only assume that the nonzero values reflect the transfers.
Proton Transfer in Liquid Phosphonic Acid
J. Phys. Chem. B, Vol. 113, No. 25, 2009 8477
The red curve in Figure 1 representing the intermolecular P-P distances has a very broad peak between 4.0 and 5.5 Å. This is the first PA shell around each PA molecule. In almost the same manner, the P-O curve shows a small maximum at 4 Å, which results from the intermolecular P-O distances. Here, we can also see that the system is in the liquid state since the intermolecular PDF curves are oscillating around 1. In contrast, the PDFs have zero values within the molecule (intramolecular PDF), which is a sign for a very ordered (here, molecular) structure. Comparing the PDFs obtained from simulations at different temperatures or with different charge carriers present in the simulation box (neutral, positive, negative setup) does not give any significant differences. Not surprisingly, the overall liquid structure is similar in these cases. However, we observe a thermal broadening of the intramolecular distance peaks resulting from thermal motion (inset of Figure 1). From the MD trajectories, we also have access to the selfdiffusion coefficients. Using Einstein’s relation
1 Dt ) 〈|r(t) - r(0)| 2〉 6
(4)
we compute the diffusion coefficients using the position of the phosphorus atoms as an approximate center of mass. For the neutral sets of simulations we have obtained self-diffusion coefficients of 1.4-2.9 × 10-7 cm2/s in the range of 350-470 K. Present charge carriers increase these values only slightly. Some experimental values for different species containing phosphonic acid groups are available. The self-diffusion coefficients of pure phosphonic acid (1.0-6.0 × 10-6 cm2/s at 350-430 K)5,35 are about 1 order of magnitude larger than the values given above. The same holds for the self-diffusion coefficients of 1-heptylphosphonic acid measured by pulsed field-gradient NMR.3 The fact that our simulations underestimate the diffusion coefficients may have two reasons. First, finite-size effects resulting from a too small simulation box may yield to selfinteractions as discussed by Yeh and Hummer.36 A second reason could be that the simulation time is not sufficient. For selected setups, we have increased the simulation time to over 100 ps with no significant change in the resulting diffusion coefficients. Most likely, the underestimation therefore will be a result of finite-size effects. Lee and Tuckerman37 have found these to be the reason for the underestimation of the diffusion coefficients of water by a factor of 4 as well. It is a general and well-known problem in MD simulations of liquids. The MD trajectories give also access to the hopping rates. To calculate these, we have defined that a proton transfer takes place if a proton is transferred from the vicinity of one phosphorus atom to the vicinity of another phosphorus atom. We take a radius of rPH ) 2.4 Å between the hydroxyl H and the P atoms. Figure 2 shows the number of proton transfers as a function of time for the neutral, positive, and negative setups. If we take the number of transfers in the neutral simulation box as a reference, an additional proton increases the transfer rate, whereas a negative charged setup results in lower values. This holds for all investigated temperatures. We observe hopping rates of 0.4-1.5 × 1013 s-1 in the investigated temperature range of 350-470 K. The proton-transfer criterion may certainly be changed; if we decrease the radius to rPH ) 2.3 Å, all hopping rates decrease by 1 to 2 orders of magnitude (∼1011 s-1). The values obtained in such a way are, therefore, only approximate. This is especially
Figure 2. Number of proton-transfer events during 5 ps at 450 K. The red, green, and blue curves count the events in the neutral, positive, and negative setups, respectively. Inset: Each proton transfer in the positive setup is marked by a vertical red line. The green curve and the scale of the ordinate are the same as those in the main figure.
true for the proton transfers in the neutral setups. Here, actual proton transfers, which would indicate a self-dissociation, are not observed. The nonzero hopping rates result from shared hydrogen atoms with small hydrogen bridge distances that are in the vicinity of two phosphorus atoms at the same time. In Figure 2, we observe additionally that at certain times, the event curve exhibits large jumps. The inset of the figure shows here many events occurring within a short time interval. These result from a proton jumping back and forth between two adjacent PA molecules (proton rattling). In other words, it is shared by two PA molecules, the positive charge is distributed over both molecules, and the whole complex is comparable to a Zundel ion (H5O2+). This situation may be called an unsuccessful proton transfer. In the neutral setup, no additional charge carrier is complexed, but two PA molecules exhibit a strong hydrogen bridge. The lower hopping rates in the negative setup are due to a charge delocalization of the negative charge over the entire H2PO3- molecule. In order to investigate the proton-transfer mechanism, we need to have an excess charge carrier. Therefore, we will use the positive setups for our discussion in the following, if not stated otherwise. Single Proton Transfers. After characterizing the liquid phosphonic acid, we will now turn toward the proton-transfer event itself. For a single proton-transfer reaction, we have tracked several properties as functions of time. Figure 3 displays different interatomic distances between atoms of the two molecules (labeled 1 and 2). The transfer takes place at t ) 5.2 ps. Here, we witness a successful transfer. From the graph, we can observe the following: (1) The distance between the two phosphorus atoms is of minor importance for the transfer. During the transfer, it does not have a particularly low value (approximately 5 Å in a range of 3.9-5.4 Å during the displayed time). (2) The O-O distance should be low (
