Theoretical analyses on water cluster structures in polymer electrolyte membrane by using dissipative particle dynamics simulations with fragment molecular orbital based effective parameters

The mesoscopic structures of polymer electrolyte membrane (PEM) affect the performances of fuel cells. Nafion® with the Teflon® backbone has been the most widely used of all PEMs, but sulfonated poly-ether ether-ketone (SPEEK) having an aromatic backbone has drawn interest as an alternative to Nafion. In the present study, a series of dissipative particle dynamics (DPD) simulations were performed to compare Nafion and SPEEK. These PEM polymers were modeled by connected particles corresponding to the hydrophobic backbone and the hydrophilic moiety of sulfonic acid group. The water particle interacting with Nafion particles was prepared as well. The crucial interaction parameters among DPD particles were evaluated by a series of calculations based on the fragment molecular orbital (FMO) method in a non-empirical way (Okuwaki et al., J. Phys. Chem. B, 2018, 122, 338–347). Through the DPD simulations, the water and hydrophilic particles aggregated, forming cluster networks surrounded by the hydrophobic phase. The structural features of formed water clusters were investigated in detail. Furthermore, the differences in percolation behaviors between Nafion and SPEEK revealed much better connectivity among water clusters by Nafion. The present FMO-DPD simulation results were in good agreement with available experimental data.


Introduction
Polymer electrolyte fuel cells (PEFC)also called polymer exchange membrane fuel cell (PEMFC)are attractive candidates for the use in vehicles because of no emission of carbon dioxide (CO 2 ). The key component of PEFC is the polymer electrolyte membrane (PEM), and the most widely used PEM is of peruorosulfonic acid (PFSA) type such as DuPont's Naon®. In addition to its high proton conductivity, Naon has chemical, thermal and mechanical stabilities. 1 Naon polymer consists of a polytetrauoroethylene backbone (or Teon®) and a side chain terminated with a sulfonic acid group. The structures formed by hydrated Naon have been investigated by various experimental methods such as X-ray scattering, 2-6 neutron scattering, 7-9 transmission electron microscopy (TEM), 10 atomic force microscopy (AFM), 11 infrared spectroscopy (IR), 12 and so on. As a result, it has been known that the hydrated Naon membrane provides nanophase-segregated structures consisting of hydrophobic phase including main chains and hydrophilic phase containing a sulfonic acid group and water. Then, the formed water cluster-networks are related with the crucial proton conductivity. [2][3][4][5][6][7][8][9][10][11][12][13][14][15][16] Radical reaction analysis using electron spin resonance has also been conducted. [17][18][19] To elucidate the hydrated structure of Naon, many molecular calculations have been carried out to date. Voth et al. analyzed the proton transport mechanism in Naon at atomistic level using the self-consistent multistate empirical valence bond (SCMS-EVB) approach. [20][21][22] Choe et al. simulated proton transport using the rst principle molecular dynamics. [23][24][25] Dupuis et al. performed a lot of molecular dynamics (MD) simulations for proton hopping and hydration of Naon. [26][27][28] Kawakami and Shigemoto also used MD simulation to verify the diffusion mechanism of proton. 29 In addition, using density functional theory (DFT) calculation, water cluster interacting with sulfonic acid group, 30 Naon/Pt interface, 31 and degradation of Naon 32 have been also reported.
Though Naon has excellent conductivity as mentioned above, other types of PFSA which have different side chain length from Naon have been studied. [33][34][35] However, there is still room for improvement on chemical durability, gas permeability, high production cost, and so on. Alternatively, the sulfonated poly ether-ether-ketone (SPEEK) having aromatic hydrocarbon has attracted attention as a promising alternative to Naon and related materials. 36,37 SPEEK demonstrates high thermal stability, mechanical properties and high cost effectiveness, but the problems of chemical stability and poor conductivity have been pointed out. 38,39 Against such a situation, Miyake et al. recently developed a new PEM composed of poly-phenylene with high chemical stability, and further improvements of the hydrocarbon PEM have been desired. 40 In order to accelerate researches to optimize such material performances of PEM, theoretical studies based on molecular simulations (that deal with large sizes and time scales) should be useful. There are mainly two ways to perform the large-scale simulations as follows.
The rst way is all-atom MD simulations using highly parallelized programs and huge computational resources. There have been various excellent programs such as MODYLAS, 41 GROMACS, 42 NAMD, 43 and LAMMPS. 44 In fact, extensive MD simulations for PEM have been reported as briefed below. Okazaki et al. investigated the morphology of the PEM using large-scale MD simulations. 45 Knox and Voth also carried out atomistic MD simulations to probe morphological models of Naon. 46 Komarov et al. performed MD simulations with cell size of up to 36 nm including 4 million atoms. 47 The second way is coarse-grained MD (CGMD) and related methods such as self-consistent mean eld method (SC-MFT) and dissipative particle dynamics (DPD). [48][49][50][51] These methods are very useful because large scale behaviors in molecular level are obtainable at reasonable computational costs. Many coarsegrained simulations have been performed for PEM. Voth et al. reported a mesoscale study of the proton transport using smoothed particle hydrodynamics (SPH). 52 The DPD method was proposed by Hoogerbrugge and Koelman, 66,67 and was extended to polymer systems by Groot et al. [48][49][50][51] However, the evaluation of a set of interaction parameters (c) among particles have still been of difficult issue. There have been a couple of major routes of c parameter predictions. The rst route is based on the solubility parameter (SP) values. The SP value method was devised by Hildebrand, 68 and various empirical estimation models such as atomic group contribution models [69][70][71][72] were developed. Note that the SP values could be obtained from molecular simulations for single molecules as well. 71 The second route is based on the evaluation of interactions between heterogeneous molecules, due to the contact energy between segments 72 and the difference in cohesive energies of the aggregated model. 73 This second route has advantage to directly determine the interactions between the contacting particles in a molecular level approach. Yamamoto and Hyodo performed novel DPD simulation using the parameter obtained from a classical force eld (FF) molecular calculation and reproduced the network structure of the water clusters generated by Naon. 74 However, since the interaction between Naon and water should involve both polarization and charge transfer, simulation results with FF might could have some limitations. For that reason, the parameters converted from experimental values were used for many of the simulations currently.
Certainly, simulations with non-empirical parameters are desirable to develop new type materials for which experimental data are hardly obtainable. In 2016, Sepehr et al. evaluated the effective interaction parameters of DPD for Naon, based on ab initio molecular orbital (MO) calculations. 75 Although ab initio evaluations of c parameters for DPD are desirable, its applicability could be limited due to the enlarged computational costs when the molecular sizes of segment pairs grows. Therefore, we have developed a new approach 76,77 to calculate such effective parameters based on the fragment molecular orbital (FMO) method. [78][79][80][81][82] This procedure could be considered as an extension of Fan's method 72 based on the Flory-Huggins theory. 83,84 In this paper, we report a series of DPD-based investigations on the morphology of hydrated Naon and SPEEK. The crucial interaction parameters among DPD particles were evaluated through our new protocol with FMO. 76,77 Furthermore, the network connectivity of water clusters was evaluated through the percolation analysis of formed water clusters in PEM, 15,58,[85][86][87][88] since such an effective index should relate with the conductivity. Naon and SPEEK were compared. The rest of this paper is congured as follows. In Section 2, the c parameter evaluations and models used for DPD simulations are described in detail. In Section 3 the simulated results are presented and discussed.

c Parameter evaluation
In the Flory-Huggins lattice theory, the free energy change (DG) for a binary system is expressed as 83,84 where 4 i and x i are the volume fraction and the chain length (i ¼ 1, 2 for the two components), respectively. The rst and the second terms on the right side of this equation describe the entropy changes, and the third term provides the enthalpy change. The c parameter is dened as where Z is the coordination number of the model lattice and the contact energy DE 12 is given by the following equation, E ij in this equation is the average interaction energy between the segments i and j in the two components, and DE 12 corresponds to the energy gain per segment due to the mixing. These relations imply a scale down of problem from mesoscale to nanoscale. Fan et al. proposed the procedure to calculate Z and DE 12 , based on MC simulations with classical FF set. 72 Recently, we have developed a new procedure 76,77 to estimate the c parameter set through a series of FMO calculations, 80,81 where the electronic effects of polarization and charge transfer were incorporated in the energy evaluation and also the molecular anisotropy was taken into account. Consequently, the c parameter could be calculated using the following equation where Z ij is the coordination number of segment j around segment i, and S ij corresponds to the scaling factor associating with anisotropy (refer to ref. 78 and 79 for details). Fig. 1 shows the structure of Naon. According to Yamamoto's previous study, 74 the basic unit of Naon chain is divided into three segments of the same size (A: The termini were capped with F for segments A and C and with CF 3 for segment B. The structure of SPEEK was shown in Fig. 2. The chain was divided into three segments A, B, and C similarly to Naon. The termini were capped with CH 3 . Furthermore, various conformations were considered for water molecules. The water particle (W) is typically modeled by a water tetramer with cyclic hydrogen bonding for DPD simulations due to the segment size problem, as was done in ref. 74. However, such a model might fail to interact with outer particles because of its internal hydrogen bonding. Therefore, we employed three kinds of dimers (with shapes of linear, cyclic, and bifurcated types) and even a monomer as the candidates interacting with sulfonic side chain (C). All the water structures employed in the parameter evaluation are shown in Fig. 3.
The geometries of segments were optimized at the dispersion-corrected B97D (ref. 89) DFT level with the 6-31G(d 0 ,p 0 ) basis set 90,91 by using the GAUSSIAN09 program. 92 Note that the orbital exponents of polarization functions of 6-31G(d 0 ,p 0 ) were optimized for the respective elements. The generation of geometrical congurations (each number was typically 2000) was performed with the J-OCTA program. 93 A number of FMO calculations for all the possible combinations among segments were carried out at FMO-MP2/6-31G(d 0 ) level, where the ABINIT-MP program 79 was used in parallel executions on several in-house servers with Intel's Xeon processors.

DPD simulation
DPD is a sort of so particle dynamics with conservative, dissipative and random forces. 66,67 The fundamental DPD scheme was extended to polymer systems by Groot et al. by introducing a bead-spring type particle model. [48][49][50] The outline of Groot's DPD model is described below. The time evolution of the given system is simulated by solving the Newton equation of motion and where, r i , v i , and m i are the position, velocity, and mass of the particle i, respectively. For convenience, the masses and diameter of the particles are scaled relative to 1. The force f i contains four parts as The rst three forces of the right-hand-side are considered within a certain cutoff radius r c . The conservative force F C ij is a so repulsion action as follows where a ij is the maximum repulsion force between particles i and j, and related denitions are r ij ¼ r j À r i , r ij ¼ |r ij |, and n ij ¼ r ij /|r ij |. The repulsion parameters between different type particles correspond to the mutual solubility provided by the c parameter set. When the reduced density r is assumed to be 3, a linear relation with c ij is usually set as 48,49 In eqn (7), the dissipative force F D ij represents the hydrodynamic drags, and the random force F R ij incorporates thermal noises of the Gaussian statistics. The fourth force F S ij is a harmonic force calculated for particles directly connected with spring bonds. The detailed functional forms of these potentials are described in eqn (S1)-(S8) in the ESI, † and the numerical values of associated parameters are compiled in Tables S1 and S2. † In the present study, the time evolution was calculated by the modied Verlet algorithm 48 with the empirical factor l ¼ 0.65  This journal is © The Royal Society of Chemistry 2018 and time step Dt ¼ 0.05. The cubic cell system was dened as a dimension of L ¼ 30 DPD length unit, corresponding to 23.5 nm in real unit (DPD particle diameter is 0.71 nm). 51 Each simulation involved about 81 000 beads, and the beads were initially packed randomly at the standard mean density r ¼ 3. 48 The cut-off radius r c was set to 1. A default value of 4 (ref. 51) was used for the spring constant of F S ij . The polymer models used are shown in Fig. 4. Three types of symbolic structures which have different equivalent weight (EW) were prepared for Naon. In Fig. 4, the structure (b) corresponds to Naon 117 (EW ¼ 1100). 33 For SPEEK, one hydrophilic segment was placed for every four particles. The time integration was performed for 10 000 steps corresponding to 500 DPD-time unit (t). The simulations were carried out for water contents of 10-30 vol% with 2 vol% intervals. The trajectory data were saved every 100 steps (5 DPD-time unit) for subsequent analyses. All the DPD simulations were performed with the COGNAC program 93,94 on in-house servers.

Results and discussion
3.1. Results of the parameter evaluation E ij and S ij between various water models and hydrophilic segment C are shown in Table 1. It is found that the E ij value varies greatly depending on the water structures. The interaction energy with the linear dimer is the strongest one, as expected from its open shape for external interactions. This linear shape of water dimer was thus actually used in evaluating the interaction energies with Naon, where the coordination number in eqn (2) was still evaluated with the tetramer. The E ij , S ij and Z ij between each segment pair at 350 K of Naon are shown in Table 2. The E ij values among the hydrophobic segments (A-A, A-B, B-B) are approximately À1 kcal mol À1 , whereas the interaction energies between the hydrophilic segments (C-C, C-W, W-W) are about À10 kcal mol À1 . On the other hand, the S ij values between hydrophobic segments are about 0.9, and those between hydrophilic segments are less  Fig. 1 for Nafion and in Fig. 2 for SPEEK. Symbolic structure (b) corresponds to Nafion 117 (EW ¼ 1100), 33 and structures (a) and (c) are defined as low and high EW models, respectively. Table 1 Interaction energies (E ij ), and scaling factors (S ij ) of segment C of Nafion shown in Fig. 1  than or equal to 0.5. These facts imply that hydrophobic segments are isotropic but hydrophilic segments are anisotropic. Table 3 shows the values of c obtained from the results listed in Table 2  The energy values with the water dimer of linear form were used for parameter evaluation for SPEEK as well as Naon. E ij , S ij and Z ij between each segment pair at 350 K of SPEEK are shown in Table 4. The interactions between hydrophobic particles are about 1.5 times larger than that of Naon. The dispersion interaction caused by the p electron of benzene rings may be responsible for this enhancement. In addition, segment C having sulfonic acid has a strong interaction with other main chain segments, and also interacts with water particles greatly. The value of c obtained from the above results is shown in Table  5. The c values among polymer segments (A-B, A-C, and B-C) were À0.75, 4.94, and 3.53 respectively, suggesting high affinity for each other. In addition, the c of C-W is À3.78, which indicates strong interaction. Thus, the tendency of the parameters predicted from the structure of SPEEK is reproduced. The nal a ij values of Naon and SPEEK for DPD are found in Table S2. † is visible that small domains of water aggregation develop into large water clusters surrounded by hydrophilic particle of Naon. Fig. 6 shows the time-dependent water density distribution (for more than 40% of water particles) of the same simulation case. The complicated network structures of water clusters are being formed from the initial uniform situation with respect to the time evolution. This progress is consistent with the visualized results of Fig. 5.

Results of DPD simulation
Next, the structural changes due to the differences in EW of Naon (corresponding to the symbolic structures of (a), (b) and (c) in Fig. 4) are investigated. The results of 20 vol% case are shown in Fig. 7. From this gure, one can see the behavior that the size of water cluster (lower row) enlarges as the chain (or length of segment A of Naon) elongates. Fig. 8 illustrates the dependence of morphology on water contents of 10, 20 and 30 vol% for the Naon, and Fig. 9 does the corresponding results for SPEEK. Comparison between Fig. 8 and 9 indicate that the water clusters formed by Naon connect each other more easily than those formed by SPEEK. Especially, the water clusters look still rather segmented at 20 vol% for SPEEK, although the mutual connections almost complete for Naon. A probable reason for this difference is the fact that there are more potential sites for hydrogen bond with water in SPEEK than in Naon. In short, Naon could provide better connectivity of water clusters than does SPEEK. 86

Analysis of Naon cluster structure
Based on the qualitative discussion on the mesoscopic structures of water clusters in the previous section, comparative discussion with available experimental data is made for Naon. Fig. 10 illustrates the radial distribution functions (RDFs) of water particles for Naon with symbolic structure (b) (Naon 117), where the three cases of water contents of 10, 20 and  30 vol% are plotted. The Fourier transformation of RDF could provide the small-angle scattering pattern, 74 as plotted in Fig. 11. The peak derived from the water-network (or so-called ionomer peak) is observed around q ¼ 0.2. When the water content increments, this peak shis to lower angles and its intensity increases. This tendency was in good agreement with the experimental data of small angle neutron scattering (SANS) (refer to Fig. 6 of ref. 8).
According to the precedent study 74 by Yamamoto and Hyodo, the diameter of water cluster was measured at the point where g(r) became 1 aer the rst peak (see again Fig. 10). The separation between the rst and second peaks was taken as the cluster space. Fig. 12 plots the dependences of cluster diameter and space on water contents of 10, 20 and 30 vol%. The incremental trend is observed for both diameter and space against the water content. For Naon 117 (symbolic structure (b) in Fig. 4), as the water content increases for 10-30 vol%, the cluster diameter changes to 3.3-5.3 nm and its spacing changes to 4.4-6.3 nm. The corresponding experimental values of the cluster diameter and distance are 4 nm and 5 nm, respectively, 2,16 being in good agreement with the present simulation results. The elongation of chain length (or the decreasing of EW) leads to the increase of both cluster size and space, and this trend is consistent with the result reported by Morohoshi et al. 60

Evaluation of diffusion coefficient
Diffusions of protons in the mesoscopic water network has been considered as a crucial indicator of performance for the PEM. 95 Fig. 5 Time-dependent (with symbol "t") morphologies of the case of 20 vol% water content for Nafion with symbolic structure (b) shown in Fig. 4. Hydrophobic, hydrophilic, and water particles are depicted with green, yellow and blue colors, respectively. The diffusion coefficients (D) were thus evaluated from the DPD results, based on the mean square displacement (MSD) of t ¼ 150-450 by using the following relation h(r(t) À r(0)) 2 i ¼ 6Dt.
(10) Fig. 13 plots the dependences of diffusion coefficient on water contents for Naon and SPEEK. For Naon, all three cases show almost linear relation against the water content. Although a linear relation is found for SPEEK as well, its slope is rather small, reecting the difference in connectivity of water clusters from Naon. In other words, Naon should provide better connections for proton diffusions. The hydration level l, which is the number of water molecules per sulfonate part, for water contents of 10, 20 and 30 vol% are 2.5, 6 and 10, respectively, for Naon 117 (structure (b)). According to experiments of Zawodzinski et al., 95 the diffusion coefficients of that region are roughly 0.8 Â 10 À6 to 4.5 Â 10 À6 cm 2 s À1 . The present DPDbased results have a corresponding range of 0.6 Â 10 À6 to 3.2 Â 10 À6 cm 2 s À1 , providing reasonable agreement with the experimental data.

Percolation analysis
So far, several differences due to the connectivity of water clusters between Naon and SPEEK have been discussed. To directly evaluate the connectivity value, the percolation analysis was performed. The size of water cluster was dened as  where R(i,j) is the distance between particle i and j, and R C is a criterion for contact. If eqn (11) is true, these two particles (i and j) belong to the same cluster. The cluster connectivity M is thus calculated as follows where N is the total number of particles in the system, and g(i) is the cluster size to which particle i belongs. The R C was set to 1.1 DPD unit length, corresponding to the spacing of the rst coordination area obtained from the RDF of DPD (see Fig. 10). For the percolation analysis, a series of additional DPD simulations were carried out for the water content range of 10-30 vol% with 2 vol% intervals. Fig. 14 presents the results of connectivity for Naon and SPEEK, where the transient structures of every 100 steps of t ¼ 300-500 were used in the evaluations. It is found that the connectivity rapidly grows over 0.8 at about 20 vol% for Naon with the symbolic structures of (a) and (b) and that there is a delay in grow for the symbolic structure (c). Finally, the connectivity value reach almost one (as full connection) as expected. These critical behaviors qualitatively accord with an experimental threshold of 10 vol%. In addition, for a given value of l (or hydration level), the connectivity became slightly higher with an increase in EW. This tendency was in good agreement with Fontanella's experimental observation. 88 On the other hand, the critical grow of connectivity occurs at about 30 vol% for SPEEK, and this tendency agrees with the experimental result by Wu et al. 86 Since the connection structure of water clusters should be   different between Naon and SPEEK, an alternative scheme was tried for the percolation analysis where the structures were averaged in every 100 steps and the R C value was loosed to 0.8 DPD unit (the rst peak position of the RDF). The corresponding results are plotted in Fig. 15. For Naon, the plots of connectivity shi toward lower water content, and a degeneracy of (a) and (b) is resolved. The correspondence to the experimental results 88 is improved by a certain extent. The results for SPEEK are lower shied, but the nal connectivity is still about 0.9 (see also Fig. S1 † of time-dependent connectivity when necessary). Both relative positioning of hydrophilic/ hydrophobic beads and length of polymer chains should affect the connectivity, as seen in the situation that the results of Naon (c) model having long chain are rather close to those of SPEEK model in Fig. 15. In summary, the percolation analysis has shed on light on the difference in mesoscale connection structures of water clusters formed by Naon and SPEEK. 86 Certainly, Naon is better as the PEM material than SPEEK.

Conclusion
In the present study, we investigated the mesoscopic structures of Naon and SPEEK by using DPD simulations. A crucial set of effective interaction (c) parameters has been determined with a series of non-empirical FMO calculations for all the possible particle-particle interactions. 76 The results of DPD simulations were analyzed in several ways, illuminating the differences in mesoscopic structures of water clusters between Naon and SPEEK. The simulation data were compared with available experimental data, and reasonable agreement have been obtained, indicating the reliability of non-empirical FMO-DPD scheme. Our analysis showed that the connectivity of water clusters formed by Naon is better that by SPEEK. This difference could be considered as a reason why Naon has been the rst choice of PEM materials industrially. Our simulation results have been in accord with other studies. 45,74,86,88 We would believe that the present FMO-DPD simulation scheme has a wide-range applicability to various mesoscopic systems. Several demonstrative applications such as lipid-water system [96][97][98] have been reported. Deng et al. have made a DPD study of sulfonic gemini surfactants. 99 Our present approach may be applicable to such novel materials. Very recently, Inokuchi et al. published an interesting paper in which DPD results of surfactants were analyzed and predicted by machine learning (ML) techniques. 100 The combination of DPD and ML should be suggestive for our future works. Finally, we would note that the sequenced protocols to evaluate a set of c parameters through the FMO calculations has been packaged as a workow system (FCEWS, FMO-based Chi-parameter Evaluation System). 77

Conflicts of interest
The authors declare no competing nancial interest.