**Molecular Dynamics Simulations of Proton Transport in Proton Exchange Membranes Based on Acid-Base Complexes**

Liuming Yan\* and Liqing Xie *Department of Chemistry, Shanghai University, Shanghai, China* 

#### **1. Introduction**

Proton exchange membrane fuel cells (PEMFCs) are promising energy conversion devices for the future society for their high energy conversion efficiency, environmental friendliness, and structural compactness. Therefore, PEMFCs have attracted great attention from academic and industry institutions; and billions of dollars of research funding from both private organizations and governmental agencies are invested into this field. However, the large-scale applications of PEMFCs are challenged by the high cost, shortage of infrastructure for the production and distribution of hydrogen, as well as performance inefficiency of the core materials including the proton exchange membranes and electrocatalysts. Though the presently most accepted proton exchange membranes based on poly(perfluorosulfonic acid), such as NAFION®, possess high proton conductivity, high chemical and electrochemical stability, excellent mechanical properties, and long service life; their performances seriously degrade under water deficiency environment or at temperatures above 80oC because of evaporation dehydration. PEMFCs operate at temperatures well above 120oC are essential as the electrocatalytic activity improves greatly at high temperature; and thus the loading of precious platinum on the electrocatalyst could be greatly reduced. Furthermore, the tolerance of electrocatalysts to impurities such as CO in fed gas is also greatly improved at high temperature and the purification cost for fed gas could be significantly reduced. Other benefits of high-temperature PEMFCs include simplified water and heat management systems, and improved overall energy conversion efficiency.

The high temperature performances of proton exchange membranes, specially the high temperature proton conductivity, must be greatly improved in order to elevate the operational temperature of PEMFCs to well above the boiling point of water. One of the methods to improve the high temperature performances of proton exchange membranes is to blend water retention materials into membranes fabricated from poly(perfluorosulfonic acid); however, the operational temperature of PEMFCs could only be increased to about

<sup>\*</sup> Corresponding Author

120oC. Since the further increase in operational temperature is difficult using proton exchange membranes based on poly(perfluorosulfonic acid) with the addition of water retention materials, proton exchange membranes being independent on hydration degree need to be developed. Recently, the proton exchange membranes based on acid-base complexes, especially blends of poly(phosphonic acids) and poly(heterocycles), are attracting more and more research attentions owing to the excellent high-temperature proton conductivity under anhydrous states. The proton conduction performances of proton exchange membranes based on acid-base complexes are determined by three essential characteristics: the concentration of transportable protons depending on the acid-base equilibrium, the morphological rearrangement depending on the structural flexibility of the polymeric molecules, and the formation of hydrogen-bond network.

Density functional theory calculations and molecular dynamics simulations of the acid-base complexes are important for the understanding of the proton transport mechanism and could provide essential insight for the rational design of high-temperature proton exchange membranes. In our previous study, we have studied the acid-base equilibrium in complexes of benzimidazole and model acids by density functional theory calculations; and the calculation results are verified by 1H NMR spectra (Zhang & Yan, 2010). We have also studied rotational flexibility of phosphonic groups for some of the phosphonic acids using density functional theory calculations (Yan, Feng et al., 2011), and the hydrogen-bond network in the pristine and phosphoric acid doped polybenzimidazole using molecular dynamics simulations (Zhu, Yan et al., 2011).

In this chapter, density functional theory calculations and molecular dynamics simulations will be applied to the study of the hydrogen bonding characteristics of acid-base complexes, as well as acidic and basic bi- and multifunctional molecules, and the formation of hydrogen-bond network in pristine, hydrated, and heterocycles solvated poly(vinylphosphonic acid). In the first section, we introduce the application background, as well as the major challenges to the proton exchange membranes. In the second section, we describe the accepted proton transport mechanisms including the essential steps for proton transport. In section three, the calculation methods, as well as molecular force field models are developed. In sections four and five, the density functional theory calculation and molecular dynamics simulation results are reported, respectively. In the last section, some of the important conclusions are summarized.

#### **2. The proton transport mechanisms in proton exchange membranes**

The proton mobility in liquid water surpasses that of any other cations since protons transport via both vehicle mechanism (molecular diffusion) and hopping mechanism (structural diffusion). However, this simple perspective complicates in proton exchange membranes because of the structural complexity and configurational restriction. For example, the structural diffusion mechanism dominates the proton conducting process in well hydrated sulfonated ionomers where continuous hydrophilic subphase exists. On the other hand, the molecular diffusion mechanism dominates in sulfonated ionomers at low hydration degree where continuous hydrophilic subphase is broken into segregated hydrophilic domains. In sulfonated ionomers at intermediate hydration degree, both vehicle and hopping mechanisms contribute substantially to the proton conducting process (Kreuer, Rabenau et al., 1982).

120oC. Since the further increase in operational temperature is difficult using proton exchange membranes based on poly(perfluorosulfonic acid) with the addition of water retention materials, proton exchange membranes being independent on hydration degree need to be developed. Recently, the proton exchange membranes based on acid-base complexes, especially blends of poly(phosphonic acids) and poly(heterocycles), are attracting more and more research attentions owing to the excellent high-temperature proton conductivity under anhydrous states. The proton conduction performances of proton exchange membranes based on acid-base complexes are determined by three essential characteristics: the concentration of transportable protons depending on the acid-base equilibrium, the morphological rearrangement depending on the structural flexibility of the

Density functional theory calculations and molecular dynamics simulations of the acid-base complexes are important for the understanding of the proton transport mechanism and could provide essential insight for the rational design of high-temperature proton exchange membranes. In our previous study, we have studied the acid-base equilibrium in complexes of benzimidazole and model acids by density functional theory calculations; and the calculation results are verified by 1H NMR spectra (Zhang & Yan, 2010). We have also studied rotational flexibility of phosphonic groups for some of the phosphonic acids using density functional theory calculations (Yan, Feng et al., 2011), and the hydrogen-bond network in the pristine and phosphoric acid doped polybenzimidazole using molecular

In this chapter, density functional theory calculations and molecular dynamics simulations will be applied to the study of the hydrogen bonding characteristics of acid-base complexes, as well as acidic and basic bi- and multifunctional molecules, and the formation of hydrogen-bond network in pristine, hydrated, and heterocycles solvated poly(vinylphosphonic acid). In the first section, we introduce the application background, as well as the major challenges to the proton exchange membranes. In the second section, we describe the accepted proton transport mechanisms including the essential steps for proton transport. In section three, the calculation methods, as well as molecular force field models are developed. In sections four and five, the density functional theory calculation and molecular dynamics simulation results are reported, respectively. In the last section, some of

**2. The proton transport mechanisms in proton exchange membranes** 

The proton mobility in liquid water surpasses that of any other cations since protons transport via both vehicle mechanism (molecular diffusion) and hopping mechanism (structural diffusion). However, this simple perspective complicates in proton exchange membranes because of the structural complexity and configurational restriction. For example, the structural diffusion mechanism dominates the proton conducting process in well hydrated sulfonated ionomers where continuous hydrophilic subphase exists. On the other hand, the molecular diffusion mechanism dominates in sulfonated ionomers at low hydration degree where continuous hydrophilic subphase is broken into segregated hydrophilic domains. In sulfonated ionomers at intermediate hydration degree, both vehicle and hopping mechanisms contribute substantially to the proton conducting process (Kreuer,

polymeric molecules, and the formation of hydrogen-bond network.

dynamics simulations (Zhu, Yan et al., 2011).

the important conclusions are summarized.

Rabenau et al., 1982).

Although there are many advantages using water as solvent including the fast proton transport dynamics, low cost, and environmental friendliness, the greatest disadvantage is the limitation to the operational temperature of fuel cells as the proton exchange membranes rapidly dehydrate at temperatures above 80oC. Therefore, it is proposed that the operational temperature could be significantly elevated if the solvent water is substituted by nonaqueous protic solvent with high boiling point, such as pyrrole, pyrazole, imidazole, phosphoric acid, and phosphonic acid. The nonaqueous protic solvent molecules usually possess negligible molecular diffusion dynamics and only structural diffusion contributes substantially to proton conductivity. For example, proton transports mainly via the hopping mechanism in concentrated phosphoric acid, where approximately 90% of the total available protons from both phosphoric acid and water (in 85% and 100% H3PO4) contribute to proton conducting, as revealed by the 1H and 31P pulsed gradient spin-echo NMR study (Chung, Bajue et al., 2000). In materials functionalized with phosphonic groups, the structural diffusion via hydrogen bonds of neighboring phosphonic groups facilitates proton transport under dehydrated state at temperatures well above the boiling point of water (Schuster, Rager et al., 2005). Therefore, the elucidation of transport mechanisms is essentially important in the design and development of better proton conducting materials applicable to fuel cells operated at high temperature under anhydrous state (Yan, Feng et al., 2011).

Fig. 1. 1-D schematic diagram for the proton hopping transport in acid-base complex: Firstly, a proton hops from phosphonic group **c** to a neighboring base group **b**. And then, the hydrogen-bond network in **c**, **d**, **e**, and **f** rearranges to accommodate the hopping. If the hydrogen-bond network rearrangement is hindered, the proton will hop back and forth without macroscopic transport. In many systems where pervasive hydrogen-bond networks exist, the hydrogen-bond network rearrangement is hindered, macroscopic proton transport is hampered. Thirdly, proton hops from **b** to **a**, **d** to **c**, … As a result, the protons hop from right to left without the macroscopic transport of any molecules.

In order to facilitate the structural diffusion of protons in nonaqueous proton conducting systems, the following criteria must be met: the existence of dissociable protons and charge defects in terms of excess or missing protons (Kreuer, Fuchs et al., 1998; Joswig and Seifert, 2009), the formation of pervasive hydrogen-bond network, and the rapid rearrangement of the hydrogen-bond network (Fig. 1). The first criterion is accomplished by the complexation of the acidic and basic functional groups, where the acidic functional groups provide the dissociable protons or excess protons, while the basic functional groups accept the dissociated protons and create the charge defects. The second criterion is accomplished by the formation of pervasive hydrogen-bond network between the hydrogen bond donors and acceptors allowing the transport of protons from hydrogen bond donors to hydrogen bond acceptors without macroscopic molecular transport. However, the hopping of protons may only represent the local vibration of protons, and the protons will hop back if no hydrogenbond network rearrangement occurs to accommodate the hopping of protons (Paddison, Kreuer et al., 2006; Brunklaus, Schauff et al., 2009). For example, in the so-called hop-turn mechanism, the effective proton transport must include the bond rotation or rearrangement of the proton defect groups (the deprotonated phosphonic group, or the hydrogen phosphonate anion) (Münch, Kreuer et al., 2001; Sevil & Bozkurt, 2004; Yamada & Honma, 2005).

#### **3. The calculation methods and molecular force field models**

#### **3.1 The density functional theory calculation method**

The density functional theory calculations are employed to study the interaction between the acidic and basic functional groups in terms of hydrogen bonding of acid-base complexes. The density functional theory calculations were performed at the B3LYP/6- 31G(d)/ 6-311++G(d, p) level of theory (Lee, Yang et al., 1988; Becke 1993) as implemented in the GAUSSIAN 03 suite of programs (Frisch, Trucks et al., 2003). Though the hybrid B3LYP functionals do not always correctly reproduce the polarizabilities and hyperpolarizabilities (Champagne, Perpete et al., 1998; Zhang, Xu et al., 2009), the charge-transfer excitation energies (Chai & Head-Gordon, 2008), the noncovalent bonding interactions (Zhao & Truhlar, 2005), and the hydrogen bonding of aromatic molecules (Elstner, Hobza et al., 2001), the B3LYP functionals have achieved the greatest success in terms of the number of published applications (Yanai, Tew et al., 2004) and are still among the best functionals that provide accurate predictions for geometries and thermochemistry of small covalent systems (Curtiss, Raghavachari et al., 1997; Song, Tokura et al., 2007; Rohrdanz, Martins et al., 2009), especially for the energetics and geometrical properties of the proton transfer and other ion-molecule reactions (Curtiss, Raghavachari et al., 1991, 1993; Merrill & Kass, 1996; Curtiss, Redfern et al., 1998; Friesner, Murphy et al., 1999). Since our major interests are inter- and intramolecular hydrogen bonding interactions and morphologies of acid-base complexes, and acidic and basic bi- and multifunctional molecules, the hybrid B3LYP functionals are still among the best choices of functionals. In addition, the calculations based on hybrid B3LYP functionals and basis set of the split-valence type with diffusive functions have good balance between accuracy and computational cost.

In the density functional theory calculations, the isolated phosphonic acids, heterocyclic compounds, phosphonic acid-heterocycle complexes, and as well as acidic and basic bi- and multifunctional molecules are optimized at the B3LYP/6-31(d) level of theory. And then,

In order to facilitate the structural diffusion of protons in nonaqueous proton conducting systems, the following criteria must be met: the existence of dissociable protons and charge defects in terms of excess or missing protons (Kreuer, Fuchs et al., 1998; Joswig and Seifert, 2009), the formation of pervasive hydrogen-bond network, and the rapid rearrangement of the hydrogen-bond network (Fig. 1). The first criterion is accomplished by the complexation of the acidic and basic functional groups, where the acidic functional groups provide the dissociable protons or excess protons, while the basic functional groups accept the dissociated protons and create the charge defects. The second criterion is accomplished by the formation of pervasive hydrogen-bond network between the hydrogen bond donors and acceptors allowing the transport of protons from hydrogen bond donors to hydrogen bond acceptors without macroscopic molecular transport. However, the hopping of protons may only represent the local vibration of protons, and the protons will hop back if no hydrogenbond network rearrangement occurs to accommodate the hopping of protons (Paddison, Kreuer et al., 2006; Brunklaus, Schauff et al., 2009). For example, in the so-called hop-turn mechanism, the effective proton transport must include the bond rotation or rearrangement of the proton defect groups (the deprotonated phosphonic group, or the hydrogen phosphonate anion) (Münch, Kreuer et al., 2001; Sevil & Bozkurt, 2004; Yamada & Honma,

**3. The calculation methods and molecular force field models** 

The density functional theory calculations are employed to study the interaction between the acidic and basic functional groups in terms of hydrogen bonding of acid-base complexes. The density functional theory calculations were performed at the B3LYP/6- 31G(d)/ 6-311++G(d, p) level of theory (Lee, Yang et al., 1988; Becke 1993) as implemented in the GAUSSIAN 03 suite of programs (Frisch, Trucks et al., 2003). Though the hybrid B3LYP functionals do not always correctly reproduce the polarizabilities and hyperpolarizabilities (Champagne, Perpete et al., 1998; Zhang, Xu et al., 2009), the charge-transfer excitation energies (Chai & Head-Gordon, 2008), the noncovalent bonding interactions (Zhao & Truhlar, 2005), and the hydrogen bonding of aromatic molecules (Elstner, Hobza et al., 2001), the B3LYP functionals have achieved the greatest success in terms of the number of published applications (Yanai, Tew et al., 2004) and are still among the best functionals that provide accurate predictions for geometries and thermochemistry of small covalent systems (Curtiss, Raghavachari et al., 1997; Song, Tokura et al., 2007; Rohrdanz, Martins et al., 2009), especially for the energetics and geometrical properties of the proton transfer and other ion-molecule reactions (Curtiss, Raghavachari et al., 1991, 1993; Merrill & Kass, 1996; Curtiss, Redfern et al., 1998; Friesner, Murphy et al., 1999). Since our major interests are inter- and intramolecular hydrogen bonding interactions and morphologies of acid-base complexes, and acidic and basic bi- and multifunctional molecules, the hybrid B3LYP functionals are still among the best choices of functionals. In addition, the calculations based on hybrid B3LYP functionals and basis set of the split-valence type with diffusive functions

In the density functional theory calculations, the isolated phosphonic acids, heterocyclic compounds, phosphonic acid-heterocycle complexes, and as well as acidic and basic bi- and multifunctional molecules are optimized at the B3LYP/6-31(d) level of theory. And then,

**3.1 The density functional theory calculation method** 

have good balance between accuracy and computational cost.

2005).

vibrational frequency calculations are conducted to verify if local minima have reached. The optimizations are satisfactory since local minima are achieved for all the isolated molecules and acid-base complexes without any imaginary frequencies. Thirdly, hydrogen bonding enthalpies are calculated at the B3LYP/6-311++(d, p) level of theory using the standard procedures taking into account of zero point energies, finite temperature corrections, and the pressure-volume work term (pV) calculated at the B3LYP/6-31(d) level of theory. Finally, various structural and energetic properties are analyzed based on the density functional theory calculations.

#### **3.2 The molecular dynamics simulation method**

The molecular dynamics simulations, applied to the study of hydrogen-bond network of the pristine, hydrated, and heterocycles solvated poly(vinylphosphonic acid), are carried out using the DL\_POLY program (Smith, Leslie et al., 2003). 3-D periodic boundary conditions are applied to all the simulation cells, and the simulations are carried out in the NPT ensemble with temperature and pressure maintained, respectively, by the Nosé-Hoover thermostat and barostat with the same relaxation time parameters of 0.2 ps (Nosé, 1984; Hoover, 1985). The Verlet integration scheme is used to solve the Newtonian equation of motion for atoms or united-atoms in the simulation cell with an integration step time of 1 fs (Verlet, 1967). The initial simulation cells are constructed by insertion of the molecules randomly. During the initial simulations, the loose initial simulation cells are gradually squeezed and reach to equilibrium densities. After reaching of the equilibrium density, another 1,000,000 steps are carried out for the generation of structural characteristics of the systems.

#### **3.3 The molecular dynamics simulation systems**

In order to evaluate the microstructure and formation of hydrogen-bond network in the acid-base complexes, eight molecular systems are simulated. The first molecular dynamics simulation system consists of 24 poly(vinylphosphonic acid) (PVPA) oligomers; and each PVPA oligomer is composed of eight repeat units (Fig. 2). The second molecular dynamics

simulation system is hydrated PVPA consisting of 24 poly(vinyl hydrogen phosphonate anion) (PVHPA) oligomers; 192 hydronium cations to compensate the PVHPA (Fig. 2), and 384 water molecules. This system has a hydration degree of three; and all the phosphonic groups are singly ionized.

The other molecular dynamics simulation systems are acid-base complexes consisting of PVPA oligomers and heterocyclic compounds frequently used as nonaqueous solvent for high-temperature proton exchange membranes. In this study, the heterocyclic compounds, including pyrrole, pyrazole, and imidazole, are used as basic components of the acid-base complexes. The solvation degrees, ratios of heterocyclic compound to phosphonic group, are one eighth and three eighths (table 1).


Table 1. The molecular dynamics simulation systems

#### **3.4 The molecular force field models**

In the molecular dynamics simulation, the total potential energy of the molecular system is partitioned into bonded and nonbonded energies. The bonded energy Ubonded is further partitioned into bond stretching, bond angle bending, and dihedral torsion energies,

$$\mathbf{U}\_{\text{bonded}} = \sum\_{l} \frac{1}{2} k\_l \left( l - l\_0 \right)^2 + \sum\_{\theta} \frac{1}{2} k\_\theta \left( \theta - \theta\_0 \right)^2 + \sum\_{\wp} A\_\wp \left( 1 + \cos(m\wp - \delta) \right) \tag{1}$$

In equation 1, the right-hand terms correspond to bond stretching, band angle bending, and dihedral torsion energies, respectively; and the summations run over all bonds, bond angles, and dihedral angles. The l, θ, and φ represent the bond length, bond angle, and dihedral angle; l0, θ0, and δ correspond to their equilibrium values; kl, kθ, and Aφ corresponds to the potential parameters; and m represents the rotational multiplicity of dihedral angle. The nonbonded energy Unonbonded is partitioned into pairwise Lennard-Jones potential and Coulombic potential,

$$\mathbf{U}\_{\text{nonbonded}} = \sum\_{\text{ij}} 4\varepsilon\_{\text{ij}} \left( \left( \frac{\mathbf{o}\_{\text{ij}}}{\mathbf{r}\_{\text{ij}}} \right)^{12} - \left( \frac{\mathbf{o}\_{\text{ij}}}{\mathbf{r}\_{\text{ij}}} \right)^{6} \right) + \sum\_{\text{ij}} \frac{\mathbf{q}\_{\text{i}} \mathbf{q}\_{\text{j}}}{\mathbf{r}\_{\text{ij}}} \tag{2}$$

In equation 2, the summations run over all atom and united-atom pairs; rij are interatomic or inter united-atomic distances; and εij and σij are Lennard-Jones potential parameters; qi and qj are residual atomic charges.

simulation system is hydrated PVPA consisting of 24 poly(vinyl hydrogen phosphonate anion) (PVHPA) oligomers; 192 hydronium cations to compensate the PVHPA (Fig. 2), and 384 water molecules. This system has a hydration degree of three; and all the phosphonic

The other molecular dynamics simulation systems are acid-base complexes consisting of PVPA oligomers and heterocyclic compounds frequently used as nonaqueous solvent for high-temperature proton exchange membranes. In this study, the heterocyclic compounds, including pyrrole, pyrazole, and imidazole, are used as basic components of the acid-base complexes. The solvation degrees, ratios of heterocyclic compound to phosphonic group, are

Systems PVPA PVHPA H2O H3O+ pyrrole pyrazole imidazole

groups are singly ionized.

I 24

one eighth and three eighths (table 1).

II 24 384 192

Table 1. The molecular dynamics simulation systems

*l*

**3.4 The molecular force field models** 

Coulombic potential,

qj are residual atomic charges.

III 24 24 IV 24 72

V 24 24 VI 24 72

VII 24 24 VIII 24 72

partitioned into bond stretching, bond angle bending, and dihedral torsion energies,

2 2

bonded 0 0

nonbonded ij

In the molecular dynamics simulation, the total potential energy of the molecular system is partitioned into bonded and nonbonded energies. The bonded energy Ubonded is further

1 1 <sup>U</sup> (1 cos( )) 2 2 *<sup>l</sup>*

In equation 1, the right-hand terms correspond to bond stretching, band angle bending, and dihedral torsion energies, respectively; and the summations run over all bonds, bond angles, and dihedral angles. The l, θ, and φ represent the bond length, bond angle, and dihedral angle; l0, θ0, and δ correspond to their equilibrium values; kl, kθ, and Aφ corresponds to the potential parameters; and m represents the rotational multiplicity of dihedral angle. The nonbonded energy Unonbonded is partitioned into pairwise Lennard-Jones potential and

σ σ q q U 4<sup>ε</sup>

In equation 2, the summations run over all atom and united-atom pairs; rij are interatomic or inter united-atomic distances; and εij and σij are Lennard-Jones potential parameters; qi and

*kl l k A m* 

> 12 6 ij ij i j

rr r

(2)

ij ij ij ij ij

 

(1)

The united-atom models are applied to the PVPA and PVHPA oligomers, where all the hydrogen atoms bonded to carbon atoms are represented implicitly by united-atoms of C1 (with one hydrogen atom) and C2 (with two hydrogen atoms) and the hydrogen atoms bonded to an oxygen atom are represented explicitly. For the hydrated poly(vinylphosphonic acid), all the phosphonic groups are singly ionized resulting hydrogen phosphonate anions. The rigid body models are applied to all the solvent molecules including water, hydronium, pyrrole, pyrazole, and imidazole. The TIP3P model is applied to water molecule (Jorgensen, Chandrasekhar et al., 1983), and the all-atom rigid body force field model is applied to hydronium cation (Urata, Irisawa et al. 2005). The molecular structure of the heterocyclic compounds are optimized by density functional theory at the B3LYP/6-31G(d) level of theory. And all the atomic and united-atomic site names are reported in figure 3.

Fig. 3. Atomic and united-atomic site names for the molecular force field models

In table 2, it summarizes the force field parameters for all the simulated molecules. The equilibrium bond lengths l0, bond angles θ0, and dihedral angles δ for PVPA and PVHPA are adapted from the optimized structures at the B3LYP/6-31G(d) level of theory (Mayo, Olafson et al., 1990; Cornell, Cieplak et al., 1995; Yan, Zhu et al., 2007; Roy, Ataol et al., 2008). The Lennard-Jones parameters of likely atomic and united-atomic pairs are also summarized in table 2. And the Lennard-Jones parameters for cross-term interactions are evaluated using the Berthelot mixing rules (Allen & Tildesley, 1990).


Table 2. Force field parameters for molecular force field models

The residual charges for atoms are evaluated by density functional theory calculations at the B3LYP/6-31G(d) of theory level using the CHELPG method (Breneman & Wiberg, 1990). The residual charges for united-atoms are summation of the carbon and bonded hydrogen atoms. All the residual charges are summarized in figure 4.

Bond stretching

Bond angle bending

Dihedral torsion Dihedral angles A*δ* (kcal· mol-1) m δ (deg) C1-C2-C1-C2 1.40 3 0 C1-C2-C1-P 1.40 3 0 C2-C1-P-O2 1.25 3 -12 C2-C1-P-OH 1.25 3 -12 C1-P-OH-HO 2.30 1 38 O2-P-OH-HO 0.25 3 0 Lennard Jones parameters Atomic or united-atomic pairs εii (kcal· mol-1) σii (Ǻ)

H-H (HN-HN, HO-HO, HW-HW) 0.0000 2.5000

The residual charges for atoms are evaluated by density functional theory calculations at the B3LYP/6-31G(d) of theory level using the CHELPG method (Breneman & Wiberg, 1990). The residual charges for united-atoms are summation of the carbon and bonded hydrogen

Table 2. Force field parameters for molecular force field models

atoms. All the residual charges are summarized in figure 4.

C1-C1 0.1450 3.9800 C2-C2 0.2150 4.0800 P-P 0.2150 4.2950 O2-O2 (OH-OH) 0.0957 3.0332 OW-OW 0.1500 3.1500 C-C 0.0950 3.8800 N-N (NR-NR) 0.1450 3.6950

Bonds l0 (Ǻ) kl (kcal· mol-1· Ǻ-2)

Bond angles θ0 (deg) kθ (kcal· mol-1· rad-2) C2-C1-C2 111.0 115.3 C2-C1-P 111.4 115.0 C1-C2-C1 111.0 115.3 C1-P-O2 111.13 145.0 C1-P-OH 108.2 145.0 O2-P-O2 119.9 280.0 O2-P-OH 108.2 90.0 OH-P-OH 105.2 100.3 P-OH-HO 113.8 93.3

C1-C2 1.526 620 C1-P 1.795 700 P-O2 1.480 1050 P-OH 1.610 460 OH-HO 0.960 1106

Fig. 4. Residual atomic and united-atomic charges for molecular force field models

Though the all-atom force field models are considered more delicate than the united-atom force field models owing to their complexity and the more force field parameters used to model the same molecular system, the united-atom force field models consist of less force field parameters, thus pertain less parameterization errors. Therefore, it cannot be concluded whether one type of force field models are better than the other type or not. The force field models applied in this study reproduce reliable molecular configurations since the geometrical parameters are based on density functional theory at the B3LYP/6-31+G(d) level of theory. At this level of theory, the calculation errors for most bond lengths and bond angles, respectively, are within 0.02 Å and 0.6o except a few large errors between 0.02~0.04 Å and 0.6~1.0o (El-Azhary & Suter, 1996; Baboul, Curtiss et al., 1999). By combination of the B3LYP/6-31G(d) geometry, CHELPG residual charges, and optimized force field parameters, the force field models used in this study minimize the number of force field parameters and retains the accuracy for the description of molecular structures and hydrogen bonding.

#### **4. Hydrogen bonding in acid-base complexes**

#### **4.1 Intermolecular hydrogen bonding**

Various polymeric phosphonic acids are widely accepted as proton exchange membranes besides the polymeric sulfonic acids. The phosphonic acids possess intermediate acidity, weaker than sulfonic acids, but stronger than carboxylic acids. The first acidic protons of phosphonic acids dissociate almost completely under hydrated state providing sufficient concentration of dissociable protons for proton conduction. Similar to concentrated phosphoric acid, phosphonic acids self-dissociate under anhydrous state and possess plausible proton conductivity even under anhydrous states and high temperature. In addition, the phosphonic acids are more stable than sulfonic acids as the C-P bonds dissociate at a higher temperature than the C-S bonds do; therefore, the PEMFCs based on polymeric phosphonic acids could operate at higher temperature than that based on sulfonic acids. On the other hand, the carboxylic acids possess much weaker acidity; their dissociation degrees are limited to a few percents; and their dissociation temperatures are lower than that of sulfonic acids. Therefore, phosphonic acids are potential candidates for the proton exchange membranes based acid-base complexes.

Dissimilar to the limited types of organic acidic groups that exist, there exist a great number of types of organic bases possible for the acid-base complexes. In this study, the heterocyclic moieties based on the 5-membered ring containing one or two nitrogen atoms in the ring, including the pyrrole, pyrazole, and imidazole, are used as basic components.

In figure 5, it summarizes the optimized structure of butanephosphonic acid (C4PA) **1**, pyrrole **2**, pyrazole **3**, and imidazole **4**, as well as acid-base complexes including C4PApyrrole **5**, C4PA-pyrazole **6**, and C4PA-imidazole **7**. For the acid-base complex **5**, the pyrrole donates its acidic proton to the phosphonic oxygen forming a hydrogen bond with a bond length of 1.960 Å. And the N-H bond length, 1.008 Å in free pyrrole **2**, is stretched to 1.019 Å. The hydrogen bond is relatively weak as the bonding enthalpy is only -5.1 kcal· mol-1. For the acid-base complex **6**, where the base component pyrazole consists of two nitrogen atoms, two hydrogen bonds are formed between the phosphonic group and the pyrazole

Fig. 5. Intermolecular hydrogen bonding in acid-base complexes at B3LYP/6-31G(d), the number in parenthesis are hydrogen bonding enthalpies in kcal· mol-1

concentration of dissociable protons for proton conduction. Similar to concentrated phosphoric acid, phosphonic acids self-dissociate under anhydrous state and possess plausible proton conductivity even under anhydrous states and high temperature. In addition, the phosphonic acids are more stable than sulfonic acids as the C-P bonds dissociate at a higher temperature than the C-S bonds do; therefore, the PEMFCs based on polymeric phosphonic acids could operate at higher temperature than that based on sulfonic acids. On the other hand, the carboxylic acids possess much weaker acidity; their dissociation degrees are limited to a few percents; and their dissociation temperatures are lower than that of sulfonic acids. Therefore, phosphonic acids are potential candidates for

Dissimilar to the limited types of organic acidic groups that exist, there exist a great number of types of organic bases possible for the acid-base complexes. In this study, the heterocyclic moieties based on the 5-membered ring containing one or two nitrogen atoms in the ring,

In figure 5, it summarizes the optimized structure of butanephosphonic acid (C4PA) **1**, pyrrole **2**, pyrazole **3**, and imidazole **4**, as well as acid-base complexes including C4PApyrrole **5**, C4PA-pyrazole **6**, and C4PA-imidazole **7**. For the acid-base complex **5**, the pyrrole donates its acidic proton to the phosphonic oxygen forming a hydrogen bond with a bond length of 1.960 Å. And the N-H bond length, 1.008 Å in free pyrrole **2**, is stretched to 1.019 Å. The hydrogen bond is relatively weak as the bonding enthalpy is only -5.1 kcal· mol-1. For the acid-base complex **6**, where the base component pyrazole consists of two nitrogen atoms, two hydrogen bonds are formed between the phosphonic group and the pyrazole

Fig. 5. Intermolecular hydrogen bonding in acid-base complexes at B3LYP/6-31G(d), the

number in parenthesis are hydrogen bonding enthalpies in kcal· mol-1

including the pyrrole, pyrazole, and imidazole, are used as basic components.

the proton exchange membranes based acid-base complexes.

with an overall bonding enthalpy of -16.0 kcal· mol-1. The hydrogen bond lengths, 1.796 and 1.715 Å, respectively, are also shorter than that in **5**. For acid-base complex **7**, the two nitrogen atoms containing in imidazole are separated by a methylidyne spacer, and can not form two hydrogen bonds simultaneously with one phosphonic group because of configuration restriction. However, the hydrogen bond length in **7** is only 1.710 Å, significantly shorter than that in **5**. The hydrogen bonding enthalpy in **7** is -12.0 kcal· mol-1 and is between that of **5** and **6**. The detailed hydrogen bonding parameters are summarized in table 3.

In order to evaluate the hydrogen bonding characteristics among poly(vinylphosphonic acid) and heterocyclic compounds, we also study acid-base complexes based on heptane-1,4 diphosphonic acid and heterocyclic compounds. From figure 6, it is revealed that the second phosphonic group separated by two methylene spacers does not have significant impact on the hydrogen bonding characteristics. The hydrogen bond lengths are similar to the corresponding system with single phosphonic group, and the hydrogen bonding enthalpies are almost doubled as the second phosphonic group is incorporated into the same molecules.

Fig. 6. Intermolecular hydrogen bonding of the acid-base complexes based on heptane-1,4 diphosphonic acid at B3LYP/6-31G(d), the number in parenthesis are hydrogen bonding enthalpies in kcal· mol-1

In summary, hydrogen bonds are formed between phosphonic groups and heterocyclic compounds: pyrrole could form one hydrogen bond; pyrazole could form two hydrogen bonds simultaneously with the same phosphonic group; and imidazole could form one hydrogen bond with the same phosphonic group despite the existence of two nitrogen atoms in the ring. In diphosphonic acids, the second phosphonic group separated by two methylene spacers does not have significant impact on hydrogen bonding characteristics.


Table 3. Hydrogen bonding parameters for the acid-base complexes (dN-H and dO-H are hydrogen bond lengths in Å, NHO is bond angle in degree, υ is vibrational wave number of the hydrogen bonded proton in cm-1, and ∆H is bonding enthalpy in kcal· mol-1).

#### **4.2 Intramolecular hydrogen bonding**

There are two processes to incorporate the acidic and basic groups into the same proton exchange membrane: by blending an acidic polymer and a basic polymer together, or by incorporating both acidic and basic groups into the same polymer. Since molecular level of mixing of two polymers, which is essential for the formation of hydrogen-bond network, is difficult, it is difficult to carry out the blending process. On the other hand, it is versatile to incorporate both acidic and basic groups into the same polymer either by copolymerization or by grafting a second functional group into a polymer.

In this section, the intramolecular hydrogen bonding between a phosphonic group and heterocyclyl group in the same molecule is discussed. In figure 7, it summarizes the optimized structure of acidic and basic bifunctional molecules including 3-(2'-pyrrolyl) butanephosphonic acid **12**, 3-(5'-pyrazolyl)-butanephosphonic acid **13**, and 3-(4'-imidazolyl) butanephosphonic acid **14**. The hydrogen bond length is 1.999 Å in **12** indicating weak hydrogen bonding. For **13**, only one hydrogen bond is formed with a bond length of 1.977 Å despite the possibility forming two hydrogen bonds. The strongest hydrogen bond is formed in **14** with a hydrogen bond length of 1.736 Å. In addition, **12**, **13**, and **14** are also intentionally optimized to structures without any hydrogen bonding as **12'**, **13'**, and **14'**. In these structures, the phosphonic group and heterocyclyl group are in opposite position, the intramolecular

7 1.710 1.014 174.4 2934 -12.0 9 1.021 1.903 169.7 3479 -14.8

10 1.031 1.813 152.3 3309 -32.2

11 1.721 1.011 173.1 2992 -21.8

12 1.017 1.999 153.1 3533 -1.3 13 1.019 1.977 137.8 3509 -0.2

14 1.736 1.006 163.6 3087 -3.1 15 1.017 2.046 157.4 3530 -2.3

16 1.017 2.051 151.2 3541 -4.7

17 1.805 0.997 158.7 3241 -12.0

There are two processes to incorporate the acidic and basic groups into the same proton exchange membrane: by blending an acidic polymer and a basic polymer together, or by incorporating both acidic and basic groups into the same polymer. Since molecular level of mixing of two polymers, which is essential for the formation of hydrogen-bond network, is difficult, it is difficult to carry out the blending process. On the other hand, it is versatile to incorporate both acidic and basic groups into the same polymer either by copolymerization

In this section, the intramolecular hydrogen bonding between a phosphonic group and heterocyclyl group in the same molecule is discussed. In figure 7, it summarizes the optimized structure of acidic and basic bifunctional molecules including 3-(2'-pyrrolyl) butanephosphonic acid **12**, 3-(5'-pyrazolyl)-butanephosphonic acid **13**, and 3-(4'-imidazolyl) butanephosphonic acid **14**. The hydrogen bond length is 1.999 Å in **12** indicating weak hydrogen bonding. For **13**, only one hydrogen bond is formed with a bond length of 1.977 Å despite the possibility forming two hydrogen bonds. The strongest hydrogen bond is formed in **14** with a hydrogen bond length of 1.736 Å. In addition, **12**, **13**, and **14** are also intentionally optimized to structures without any hydrogen bonding as **12'**, **13'**, and **14'**. In these structures, the phosphonic group and heterocyclyl group are in opposite position, the intramolecular

Systems dN-H dO-H NHO υ ∆H 5 1.019 1.960 144.1 3500 -5.1 6 1.030 1.796 152.6 3327 -16.0

1.715 1.014 170.3 2945

1.023 1.916 151.7 3427

 1.777 1.005 165.4 3117 1.031 1.789 153.2 3313 1.707 1.015 170.3 2913

1.707 1.016 169.8 2919

2.588 0.977 160.4 3671

1.016 2.133 133.1 3551

1.019 1.927 147.7 3507

1.709 1.012 163.1 2984

Table 3. Hydrogen bonding parameters for the acid-base complexes (dN-H and dO-H are hydrogen bond lengths in Å, NHO is bond angle in degree, υ is vibrational wave number

of the hydrogen bonded proton in cm-1, and ∆H is bonding enthalpy in kcal· mol-1).

**4.2 Intramolecular hydrogen bonding** 

or by grafting a second functional group into a polymer.

interaction between the phosphonic group and heterocyclyl group is negligible. By comparison of the two sets of structures, **12** and **12'**, **13** and **13'**, and **14** and **14'**, the intramolecular hydrogen bonding enthalpies are evaluated as the enthalpy difference of the two sets of structures. The calculated hydrogen bonding enthalpies for **12**, **13**, and **14** are only - 1.3, -0.2, and -3.1 kcal· mol-1, respectively, indicating marginal intramolecular hydrogen bonding. Therefore, it is concluded that the intramolecular hydrogen bonding between a phosphonic group and a heterocyclyl group separated by two methylene spacers are hindered by configuration restriction. Since these model molecules could represent copolymers of the vinylphosphonic acid and vinyl heterocycles, heterocyclyl grafted poly(vinyl phosphonic acid) or phosphonic group grafted poly(vinyl heterocycles), and polymers grafted with heterocyclyl and phosphonic group, it could be further concluded that intramolecular hydrogen bonding in bi- and multifunctional polymers is hindered by configuration restriction and intramolecular hydrogen-bond network is unlikely to be formed.

Fig. 7. Intramolecular hydrogen bonding in acid-base complexes based on bifunctional molecules optimized at B3LYP/6-31G(d), the number in parenthesis are hydrogen bonding enthalpies in kcal· mol-1

In figure 8, it summarizes the structures of some multifunctional molecules grafted with two phosphonic groups and two heterocyclyl groups, including the 2,6-di(2'-pyrryl)-heptane-1,4 -diphosphonic acid **15**, 2,6-di (5'-pyrazyl)-heptane-1,4-diphosphonic acid **16**, and 2,6-di(4' imidazyl)-heptane-1,4-diphosphonic acid **17**. These multifunctional molecules are optimized to form the maximum intramolecular hydrogen bonds between the phosphonic group and heterocyclyl group. The intramolecular hydrogen bonds formed in the multifunctional molecules possess similar bond lengths to that in the bifunctional molecules **12**, **13**, and **14**, indicating that grafting more functional groups to the molecule exerts little impact on the intramolecular hydrogen bonding characteristics. In addition, these multifunctional

Fig. 8. Intramolecular hydrogen bonding of acid-base complexes based on multifunctional molecules optimized at B3LYP/6-31G(d), the number in parenthesis are hydrogen bonding enthalpies in kcal· mol-1

In figure 8, it summarizes the structures of some multifunctional molecules grafted with two phosphonic groups and two heterocyclyl groups, including the 2,6-di(2'-pyrryl)-heptane-1,4 -diphosphonic acid **15**, 2,6-di (5'-pyrazyl)-heptane-1,4-diphosphonic acid **16**, and 2,6-di(4' imidazyl)-heptane-1,4-diphosphonic acid **17**. These multifunctional molecules are optimized to form the maximum intramolecular hydrogen bonds between the phosphonic group and heterocyclyl group. The intramolecular hydrogen bonds formed in the multifunctional molecules possess similar bond lengths to that in the bifunctional molecules **12**, **13**, and **14**, indicating that grafting more functional groups to the molecule exerts little impact on the intramolecular hydrogen bonding characteristics. In addition, these multifunctional

Fig. 8. Intramolecular hydrogen bonding of acid-base complexes based on multifunctional molecules optimized at B3LYP/6-31G(d), the number in parenthesis are hydrogen bonding

enthalpies in kcal· mol-1

molecules are intentionally optimized to structures without or only one pair of acidic and basic groups in intramolecular hydrogen bonding configurations shown as **15'**, **16'**, and **17'** in figure 8. The hydrogen bonding enthalpies, enthalpy differences between **15** and **15'**, **16** and **16'**, and **17** and **17'**, are only -2.3, -4.7, and -12.0 kcal·mol-1 for the multifunctional molecules, respectively. Therefore, the intramolecular hydrogen bonding is weak in the multifunctional molecules.

From the density functional theory calculations, it is concluded that the intramolecular hydrogen bonds in the bi- and multifunctional molecules are relatively weak compared to the intermolecular hydrogen bonds in the acid-base complexes with the same acidic and basic functional groups. The advantages to apply bi- and multifunctional molecules for proton exchange membranes is not the formation of intramolecular hydrogen bonds, but is the formation of intermolecular hydrogen bonds as the molecular level of mixing of acidic and basic groups is automatically achieved in these systems.

#### **5. Hydrogen-bond network in acid-base complexes**

#### **5.1 The pristine and hydrated PVPA**

The PVPA is a highly hydrophilic polymeric material dissolvable in water. Since the phosphonic group is an acid with intermediate strength, the first acidic protons are almost fully dissociated and the second acidic protons are usually not dissociated under hydrated state. For the simplicity of molecular force field model, all the phosphonic groups are supposed to be singly ionized resulting the poly(vinyl hydrogen phosphonate anion) under hydrated state. In order to compare the structural differences between pristine and hydrated states, both pristine and hydrated PVPA are studied using molecular dynamics simulations.

In figure 9, it shows the typical structural snapshots of the pristine and hydrated PVPA revealed by molecular dynamics simulations. For the pristine PVPA, the oligomers are well ordered and densely packed with an equilibrium density at about 1.67 g· cm-3. In addition, most of the molecular backbones are stretched in the same direction because of the squeezing from the neighboring PVPA oligomers. The phosphonic groups are interconnected in a hydrogen-bond network by donating the acidic protons HO to the acidic oxygen atoms O2 or OH as shown by dotted lines in figure 9a. Furthermore, the radial distribution functions shown in figure 10a reveal that the hydrogen bonds formed between O2 and HO are much stronger than that formed between OH and HO in terms of peak heights at 7.54 and 1.18 for O2-HO and OH-HO, respectively. The peak positions for the radial distribution functions of O2-HO at 1.425 Å is slightly shorter than that of OH-HO at 1.525 Å. These characteristics are attributed to the strong steric repulsion from OH compared to that from O2. And proton could approach to O2 from all directions but to OH from only half of the directions. In fact, the hydrogen bonds between OH and HO are very weak as a peak height of 1.18 is only slightly higher than the average density of unity.

If water molecules are blended into PVPA, the phosphonic groups are well hydrated forming hydrogen phosphonate anions and separated by coulombic repulsion. Since the squeezing effect from neighboring PVPA oligomers is absent in hydrated states, the molecular backbones become less stretched compared to that in pristine PVPA and trend to distort randomly (Fig. 9b). The water molecules have also great impact on the hydrogen bonding characteristics of PVPA. The strongest hydrogen bonds are formed between O2 and HO in pristine PVPA; however, the O2-HO hydrogen bonds are greatly weakened and become insignificant in hydrated PVPA. Actually, all the hydrogen bonds formed between two hydrogen phosphonate anions become insignificant compared to that between hydrogen phosphonate anion and water in terms of peak heights of the corresponding radial distribution functions. In hydrated PVPA, the strongest hydrogen bonds are formed between O2 and HW; and the next strongest hydrogen bonds are formed between water molecules (OW and HW). However, the hydrogen bonds between OW and HW are quite weak compared to that between O2 and HW as the peak height of radial distribution function for OW-HW is only 1.21, about one fourth of that for O2-HW. Therefore, it could be concluded that both hydrogen bonding between PVPA oligomers and between water molecules are greatly weakened during hydration.

Fig. 9. Molecular dynamics simulation snapshots of (a) the pristine and (b) hydrated PVPA. Color code: red for O, white for H, green for P, cyan for C, and the red dotted lines for hydrogen bonds.

Fig. 10. Radial distribution functions showing hydrogen bonding in (a) pristine and (b) hydrated PVPA

HO in pristine PVPA; however, the O2-HO hydrogen bonds are greatly weakened and become insignificant in hydrated PVPA. Actually, all the hydrogen bonds formed between two hydrogen phosphonate anions become insignificant compared to that between hydrogen phosphonate anion and water in terms of peak heights of the corresponding radial distribution functions. In hydrated PVPA, the strongest hydrogen bonds are formed between O2 and HW; and the next strongest hydrogen bonds are formed between water molecules (OW and HW). However, the hydrogen bonds between OW and HW are quite weak compared to that between O2 and HW as the peak height of radial distribution function for OW-HW is only 1.21, about one fourth of that for O2-HW. Therefore, it could be concluded that both hydrogen bonding between PVPA oligomers and between water

molecules are greatly weakened during hydration.

hydrogen bonds.

hydrated PVPA

(a) (b)

(a) (b)

Fig. 10. Radial distribution functions showing hydrogen bonding in (a) pristine and (b)

Fig. 9. Molecular dynamics simulation snapshots of (a) the pristine and (b) hydrated PVPA. Color code: red for O, white for H, green for P, cyan for C, and the red dotted lines for

The limited hydrogen bonding between water molecules are attributed to the low hydration degree of molecular simulation system as most of the water molecules and hydronium cations are attracted by the hydrogen phosphonate anions, and a few water molecules or hydronium cations are free to form water clusters via intermolecular hydrogen bonding. These characteristics are significantly different from that of hydrated poly(perfluorosulfonic acid) where water molecules, hydronium cations and sulfonate anions aggregate forming hydrophilic clusters and channels; and molecular backbones aggregate forming hydrophobic subphase. The hydrophilic clusters and channels provide uninterrupted proton transport channels facilitating the fast macroscopic transport of protons and the hydrophobic subphase provides the membrane with integrity and mechanical properties. The lack of phase segregation in the hydrated PVPA compared to the hydrated poly(perfluorosulfonic acid) are structural bases explaining the great differences between these two types of proton exchange materials.

The hydrogen bond length also shows significant changes during hydration. For example, the hydrogen bond length for O2 and HW is about 1.425 Å, same as that of O2 and HO in pristine PVPA; on the other hand, the hydrogen bond length for OH and HW is about 1.775 Å compared to 1.525 Å for OH and HO in pristine PVPA. Furthermore, the second coordination layer of O2-HW shows high order as the peak height for the second coordination layer at 2.975 Å is 2.35, almost half of the first peak. Ordered second coordination layer also exists for OW-HW, OW-HO, and OH-HW, but absents for O2-HO and OH-HO. Therefore, it is concluded that two layers of water molecules and hydronium cations surround the hydrogen phosphonate anion and the other hydrogen phosphonate anions are excluded from the first two coordination layers.

#### **5.2 The pyrrole solvated PVPA**

In order to increase the operational temperature of proton exchange membrane fuel cells, heterocyclic compounds with high boiling points and low saturated vapor pressures, such as pyrrole, pyrazole, and imidazole, are often employed as protic solvent to substitute water in the solvation of proton exchange membranes. These heterocyclic compounds possess high boiling points and the fuel cells are operational at temperature well above the boiling point of water. On the other hand, water boils at 100oC and begins to evaporate significantly even below 100oC, thus limiting the fuel cell operational temperature to about 80oC.

The microstructures of PVPA solvated by the pyrrole, at solvation degrees of one eighth and three eighths, are shown in figure 11. From the microstructures, it is shown that the PVPA oligomers are well solvated by pyrrole molecules and the pyrrole molecules are evenly dispersed in the PVPA oligomers forming a continuous phase without phase segregation. And hydrogen bonds are formed between different phosphonic groups and between pyrrole and phosphonic group.

For the pyrrole solvated PVPA at solvation degree of one eighth, the strongest hydrogen bonds are formed between acidic oxygen O2 and acidic hydrogen HO of the phosphonic groups in terms of peak height of the corresponding radial distribution functions (fig. 12). For example, the peak height for O2-HO is 9.67, about six times of the second highest peak for OH-HO at 1.54. The second strongest hydrogen bonds are formed between the OH and HO, still between two phosphonic groups. In addition, hydrogen bonds are also formed between a phosphonic group and a pyrrole molecule, where both oxygen O2 and OH accept acidic hydrogen from pyrrole. However, the hydrogen bonds formed between phosphonic group and pyrrole are much weaker than that formed between two phosphonic groups as the peak heights for O2-HN and OH-HN are only 1.10 and 0.84, respectively. The hydrogen bond lengths also show systematic changes. For O2-HO and OH-HO, the hydrogen bond lengths are at 1.375 and 1.525 Å, respectively. These hydrogen bond lengths are significantly closer than the corresponding hydrogen bond lengths in the pristine and hydrated PVPA. The change in hydrogen bond lengths is attributed to the exclusion effect, the interaction between solvent molecules or between PVPA oligomers are more powerful than that between solvent and PVPA, resulting the exclusion of PVPA from the solvent.

Fig. 11. Molecular dynamics simulation snapshots of the pyrrole solvated PVPA at solvation degree of (a) one eighth and (b) three eighths. Color code: red for O, white for H, green for P, cyan for C, blue for N, and the red dotted lines for hydrogen bonds.

Fig. 12. Radial distribution functions showing hydrogen bonding in the pyrrole solvated PVPA at solvation degree of (a) one eighth and (b) three eighths.

For the pyrrole solvated PVPA at solvation degree of three eighths, the overall microstructures do not show significant changes compared to that at solvation degree of one eighth except the relative hydrogen bond strengths and lengths. For example, the peak height for O2-HO increases from 9.67 to 10.30, about 6.5 % increase, as the solvation degree increases from one eighth to three eighths. The exclusion effect reinforces the hydrogen bonds between phosphonic groups in the pyrrole solvated PVPA at high solvation degree compared to that at low solvation degree. In addition, the peak height for O2-HN at 2.43 is now higher than that for OH-HO at 1.51, while at solvation degree of one eighth the peak height for OH-HO is higher than that for O2-HN, despite that the peak position for O2-HN at 1.975 Å is still farther than that for OH-HO at 1.525 Å.

Special attention should be given to the hydrogen bonding characteristics of pyrrole nitrogen. Though amine nitrogen accepts proton in hydrogen bonding, the pyrrole nitrogen does not accept acidic protons from phosphonic groups nor from other pyrrole molecules since the pyrrole nitrogen is in sp2 hybrid forming a planar structure. The planar structure restrains the approaching of proton in hydrogen bonding.

#### **5.3 The pyrazole solvated PVPA**

352 Molecular Interactions

between a phosphonic group and a pyrrole molecule, where both oxygen O2 and OH accept acidic hydrogen from pyrrole. However, the hydrogen bonds formed between phosphonic group and pyrrole are much weaker than that formed between two phosphonic groups as the peak heights for O2-HN and OH-HN are only 1.10 and 0.84, respectively. The hydrogen bond lengths also show systematic changes. For O2-HO and OH-HO, the hydrogen bond lengths are at 1.375 and 1.525 Å, respectively. These hydrogen bond lengths are significantly closer than the corresponding hydrogen bond lengths in the pristine and hydrated PVPA. The change in hydrogen bond lengths is attributed to the exclusion effect, the interaction between solvent molecules or between PVPA oligomers are more powerful than that

between solvent and PVPA, resulting the exclusion of PVPA from the solvent.

(a) (b)

P, cyan for C, blue for N, and the red dotted lines for hydrogen bonds.

(a) (b)

PVPA at solvation degree of (a) one eighth and (b) three eighths.

Fig. 12. Radial distribution functions showing hydrogen bonding in the pyrrole solvated

Fig. 11. Molecular dynamics simulation snapshots of the pyrrole solvated PVPA at solvation degree of (a) one eighth and (b) three eighths. Color code: red for O, white for H, green for

Although pyrrole is a heterocyclic compound containing one nitrogen heteroatom in the ring, the nitrogen does not accept proton in hydrogen bonding because of its planar structure. In order to relay protons in proton hopping, the compound should not only donate but also accept proton in hydrogen bonding. Therefore, another nitrogen heteroatom which could accept proton in hydrogen bonding has to be included in the heterocyclic compound. Pyrazole and imidazole are two 5-membered heterocyclic compounds containing two nitrogen atoms in the ring; one nitrogen atom accepts proton and the other nitrogen atom donates its proton in hydrogen bonding. In this and the following section, the pyrazole solvated PVPA and imidazole solvated PVPA will be simulated.

Fig. 13. Molecular dynamics simulation snapshots for the pyrazole solvated PVPA at solvation degree of (a) one eighth and (b) three eighths. Color code: red for O, white for H, green for P, cyan for C, blue for N, and the red dotted lines for hydrogen bonds.

Fig. 14. Radial distribution functions showing hydrogen bonding in the pyrazole solvated PVPA at solvation degree of (a) one eighth and (b) three eighths.

In figure 13, the molecular dynamics simulation snapshots of the pyrazole solvated PVPA are shown at two levels of solvation degrees, and it could be seen that the pyrazole molecules are dispersed evenly in the PVPA for both systems.

From the radial distribution functions (fig. 14), it is concluded that the strongest hydrogen bonds are formed between the acidic oxygen O2 and acidic hydrogen HO of phosphonic groups with a hydrogen bond length of 1.425 Å. The next two types of hydrogen bonds are formed between an acidic proton HO of a phosphonic group and an imine nitrogen NR of a pyrazole molecule, and between an acidic oxygen OH and an acidic hydrogen HO of phosphonic groups. However, the hydrogen bond length between NR and HO is farther than that between OH and HO owing to the steric repulsion between a phosphonic group and a pyrazole ring. The hydrogen bonds are also formed between acidic oxygen O2 and amine hydrogen HN. The hydrogen bond lengths are divided into two groups: the O2-HN and NR-HO, and the O2-HO and OH-HO. For the pyrazole solvated PVPA at solvation degree of three eighths, the exclusion effect which is similar to that of pyrrole solvated PVPA reinforces the hydrogen bonding between phosphonic groups which are excluded by the pyrazole molecules. The increase in hydration degree also enhances the hydrogen bonding between O2 and HN.

Differing from the pyrrole solvated PVPA, there are pervasive hydrogen-bond network formed in pyrazole solvated PVPA: a phosphonic group not only accepts but also donates proton simultaneously; and a pyrazole molecule also forms two hydrogen bonds simultaneously by donating its amine hydrogen HN and accepting a phosphonic acidic proton by the imine nitrogen NR. Therefore, uninterrupted proton transport channels are developed among phosphonic groups and pyrazole molecules in terms of pervasive hydrogen-bond network and proton hopping is facilitated.

#### **5.4 The imidazole solvated PVPA**

354 Molecular Interactions

(a) (b)

PVPA at solvation degree of (a) one eighth and (b) three eighths.

molecules are dispersed evenly in the PVPA for both systems.

hydrogen-bond network and proton hopping is facilitated.

bonding between O2 and HN.

Fig. 14. Radial distribution functions showing hydrogen bonding in the pyrazole solvated

In figure 13, the molecular dynamics simulation snapshots of the pyrazole solvated PVPA are shown at two levels of solvation degrees, and it could be seen that the pyrazole

From the radial distribution functions (fig. 14), it is concluded that the strongest hydrogen bonds are formed between the acidic oxygen O2 and acidic hydrogen HO of phosphonic groups with a hydrogen bond length of 1.425 Å. The next two types of hydrogen bonds are formed between an acidic proton HO of a phosphonic group and an imine nitrogen NR of a pyrazole molecule, and between an acidic oxygen OH and an acidic hydrogen HO of phosphonic groups. However, the hydrogen bond length between NR and HO is farther than that between OH and HO owing to the steric repulsion between a phosphonic group and a pyrazole ring. The hydrogen bonds are also formed between acidic oxygen O2 and amine hydrogen HN. The hydrogen bond lengths are divided into two groups: the O2-HN and NR-HO, and the O2-HO and OH-HO. For the pyrazole solvated PVPA at solvation degree of three eighths, the exclusion effect which is similar to that of pyrrole solvated PVPA reinforces the hydrogen bonding between phosphonic groups which are excluded by the pyrazole molecules. The increase in hydration degree also enhances the hydrogen

Differing from the pyrrole solvated PVPA, there are pervasive hydrogen-bond network formed in pyrazole solvated PVPA: a phosphonic group not only accepts but also donates proton simultaneously; and a pyrazole molecule also forms two hydrogen bonds simultaneously by donating its amine hydrogen HN and accepting a phosphonic acidic proton by the imine nitrogen NR. Therefore, uninterrupted proton transport channels are developed among phosphonic groups and pyrazole molecules in terms of pervasive Imidazole is an isomer of pyrazole containing two nitrogen atoms in the 5-membered heterocyclic ring. Imidazole is more stable than pyrazole owing to the absent of N-N bond in the ring. In addition, the boiling point of imidazole at 256oC is also much higher than that of pyrazole at 187oC. Since the two nitrogen atoms in pyrazole are in neighboring position, two hydrogen bonds are favored between two pyrazole molecules thus restricting the formation of intermolecular hydrogen bonds with other molecules. On the other hand, the two nitrogen atoms in imidazole are separated by a methylene spacer, the formation of two intermolecular hydrogen bonds with other two molecules simultaneously is favored. The structural difference also causes the differences in relative acidity and basicity of the molecules. For pyrazole, the pKa is 14.0 and the pKb is 2.5, or the basicity is stronger than its acidity. On the other hand, the imidazole possesses a much balanced pKa and pKb at 6.95 and 7.05, respectively, and is almost neutral. Therefore, it is expectable that the hydrogen bonding characteristics in imidazole solvated PVPA significantly differ from that in pyrazole solvated PVPA.

The microstructure of imidazole solvated PVPA is similar to that of pyrazole solvated PVPA. In addition, the hydrogen bonds formed between different phosphonic groups, O2- HO and OH-HO, are similar to that of pyrazole solvated PVPA. Furthermore, the hydrogen bonds formed between the O2 and HN are also similar to that of pyrazole solvated PVPA. The significant difference is the hydrogen bonding characteristics between HO and NR as the imine nitrogen NR in the imidazole rarely accepts acidic proton in hydrogen bonding compared to the powerful proton accepting ability of the pyrazole imine nitrogen NR. This characteristic is important for the proton transport as the hydrogen-bond network formed in imidazole solvated PVPA is less balanced than that formed in the pyrazole solvated PVPA.

Fig. 15. Molecular dynamics simulation snapshots for the imidazole solvated PVPA at solvation degree of (a) one eighth and (b) three eighths. Color code: red for O, white for H, green for P, cyan for C, blue for N, and the red dotted lines for hydrogen bonds.

Fig. 16. Radial distribution functions showing hydrogen bonding in the imidazole solvated PVPA at solvation degree of (a) one eighth and (b) three eighths

#### **6. Conclusions**

Density functional theory calculations and molecular dynamics simulations are applied to the study of proton hopping characteristics in acid-base complexes in terms of hydrogen bonding and hydrogen-bond network. The density functional theory calculations are applied to the study of intra- and intermolecular hydrogen bonding characteristics of acidbase complexes composed of phosphonic acids and heterocyclic compounds, as well as acidic and basic bi- and multifunctional molecules consisting of phosphonic groups and heterocyclyl groups. It is concluded that the intermolecular hydrogen bonds in acid-base complexes are powerful as the entire configuration space could be searched in hydrogen bonding without restriction. On the other hand, the intramolecular hydrogen bonds in the bi- and multifunctional molecules are weak compared to the intermolecular hydrogen bonds because of configurational restriction in hydrogen bonding.

The molecular dynamics simulations based on hybrid force field models, united-atom models for the PVPA and PVHPA and all-atom rigid body models for solvents, are applied to the study of microstructures and hydrogen-bond networks of acid-base complexes including pristine and hydrated PVPA, and pyrrole, pyrazole, or imidazole solvated PVPA. It is revealed that the PVPA oligomers are densely packed with high density in pristine PVPA; while the PVPA oligomers are well hydrated or solvated without phase segregation in the hydrated or solvated PVPA.

It is further concluded that the pervasive hydrogen-bond networks vary significantly for different simulation systems. Firstly, the hydrogen bonds, O2-HO and OH-HO, existing in pristine PVPA, are greatly weakened as the PVPA is hydrated or solvated. In hydrated PVPA, strong hydrogen bonds are formed between phosphonate anion and water molecule even up to well ordered second coordination layer. Thirdly, in heterocycle solvated PVPA, exclusion effect, the exclusion of phosphonic groups from solvent, is recognized. The exclusion of phosphonic groups reinforces the hydrogen-bonding networks. Fourthly, the proton transport channels differ in different heterocycle solvated PVPA systems: the proton transport channels formed in pyrrole solvated PVPA are interrupted, that formed in pyrazole or imidazole solvated PVPA are uninterrupted, and that formed in pyrazole solvated PVPA are more balanced than that formed in imidazole solvated PVPA.

From these researches, it is proposed that special attention should be given to multifunctional polymers or copolymers with both acidic and basic functional groups attached on the same polymeric molecule in future study.

## **7. Acknowledgements**

The authors acknowledge the support from the Chinese National Science Foundation (Nos. 20873081, 21073118), and the Nano project (No. 0952nm01300) of Shanghai Municipal Science & Technology Commission.

#### **8. References**

356 Molecular Interactions

(a) (b)

PVPA at solvation degree of (a) one eighth and (b) three eighths

because of configurational restriction in hydrogen bonding.

in the hydrated or solvated PVPA.

**6. Conclusions** 

Fig. 16. Radial distribution functions showing hydrogen bonding in the imidazole solvated

Density functional theory calculations and molecular dynamics simulations are applied to the study of proton hopping characteristics in acid-base complexes in terms of hydrogen bonding and hydrogen-bond network. The density functional theory calculations are applied to the study of intra- and intermolecular hydrogen bonding characteristics of acidbase complexes composed of phosphonic acids and heterocyclic compounds, as well as acidic and basic bi- and multifunctional molecules consisting of phosphonic groups and heterocyclyl groups. It is concluded that the intermolecular hydrogen bonds in acid-base complexes are powerful as the entire configuration space could be searched in hydrogen bonding without restriction. On the other hand, the intramolecular hydrogen bonds in the bi- and multifunctional molecules are weak compared to the intermolecular hydrogen bonds

The molecular dynamics simulations based on hybrid force field models, united-atom models for the PVPA and PVHPA and all-atom rigid body models for solvents, are applied to the study of microstructures and hydrogen-bond networks of acid-base complexes including pristine and hydrated PVPA, and pyrrole, pyrazole, or imidazole solvated PVPA. It is revealed that the PVPA oligomers are densely packed with high density in pristine PVPA; while the PVPA oligomers are well hydrated or solvated without phase segregation

It is further concluded that the pervasive hydrogen-bond networks vary significantly for different simulation systems. Firstly, the hydrogen bonds, O2-HO and OH-HO, existing in pristine PVPA, are greatly weakened as the PVPA is hydrated or solvated. In hydrated PVPA, strong hydrogen bonds are formed between phosphonate anion and water molecule


Curtiss, L. A., Raghavachari, K. et al. (1993). Gaussian-2 theory using reduced Møller–

Curtiss, L. A., Raghavachari, K. et al. (1997). Assessment of Gaussian-2 and density

Curtiss, L. A., Redfern, P. C. et al. (1998). Assessment of Gaussian-2 and density functional

El-Azhary, A. A. & Suter, H. U. (1996). Comparison between optimized geometries and

Elstner, M., Hobza, P. et al. (2001). Hydrogen bonding and stacking interactions of nucleic

Friesner, R. A., Murphy, R. B. et al. (1999). Correlated ab initio electronic structure

Frisch, M. J., Trucks, G. W. et al. (2003). GAUSSIAN 03, Revision C.2. Wallingford, CT,

Hoover, W. G. (1985). Canonical dynamics: Equilibrium phase-space distributions. *Phys.* 

Jorgensen, W. L., Chandrasekhar, J. et al. (1983). Comparison of simple potential functions

Joswig, J.-O. & Seifert, G. (2009). Aspects of the Proton Transfer in Liquid Phosphonic Acid.

Kreuer, K. D., Fuchs, A. et al. (1998). Imidazole and pyrazole-based proton conducting

Kreuer, K. D., Rabenau, A. et al. (1982). Vehicle mechanism; a new model for the

Lee, C., Yang, W. et al. (1988). Development of the Colle-Salvetti correlation-energy formula

Mayo, S. L., Olafson, B. D. et al. (1990). DREIDING: A genericforce field for molecular

Merrill, G. N. & Kass, S. R. (1996). Calculated gas-phase acidities using density functional

Münch, W., Kreuer, K. D. et al. (2001). The diffusion mechanism of an excess proton in

Nosé, S. (1984). A unified formulation of the constant temperature molecular dynamics

Paddison, S. J., Kreuer, K. D. et al. (2006). About the choice of the protogenic group in

imidazole molecule chains: first results of an ab initio molecular dynamics study.

polymer electrolyte membranes: Ab initio modelling of sulfonic acid, phosphonic acid, and imidazole functionalized alkanes. *Phys. Chem. Chem. Phys.* 8(39): 4530-

into a functional of the electron density. *Phys. Rev. B* 37(2): 785-789.

interpretation of the conductivity of fast proton conductors. *Angew. Chem. Int. Ed.* 

calculations for large molecules. *J. Phys. Chem. A* 103(13): 1913-1928.

for simulating liquid water. *J. Chem. Phys.* 79(2): 926-935.

polymers and liquids. *Electrochim. Acta* 43(10-11): 1281-1288.

functional theories for the computation of enthalpies of formation. *J. Chem. Phys.* 

theories for the computation of ionization potentials and electron affinities. J.

vibrational frequencies calculated by the DFT methods. *J. Phys. Chem.* 100(37):

acid base pairs: A density-functional-theory based treatment. *J. Chem. Phys.* 114(12):

Plesset orders. *J. Chem. Phys.* 98(2): 1293-1298.

106(3): 1063-1079.

15056-15063.

5149-5155.

Gaussian, Inc.

*Engl.* 21(3): 208.

4542.

*Rev. A* 31(3): 1695-1697.

*J. Phys. Chem. B* 113(25): 8475-8480.

simulations. *J. Phys. Chem.* 94(26): 8897-8909.

*Solid State Ionics* 145(1-4): 437-443.

methods. *J. Chem. Phys.* 81(1): 511-519.

theory: Is it reliable? *J. Phys. Chem.* 100(44): 17465-17471.

Chem. Phys. 109(1): 42-45.


Zhu, S., Yan, L. et al. (2011). Molecular dynamics simulation of microscopic structure and hydrogen bond network of the pristine and phosphoric acid doped polybenzimidazole. *Polymer* 52(3): 881-892.
