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

Abstract

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.

---

1. 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₂). The key component of PEFC is the polymer electrolyte membrane (PEM), and the most widely used PEM is of perfluorosulfonic acid (PFSA) type such as DuPont's Nafion®. In addition to its high proton conductivity, Nafion has chemical, thermal and mechanical stabilities.¹ Nafion polymer consists of a polytetrafluoroethylene backbone (or Teflon®) and a side chain terminated with a sulfonic acid group. The structures formed by hydrated Nafion have been investigated by various experimental methods such as X-ray scattering,^2–6 neutron scattering,^7–9 transmission electron microscopy (TEM),¹⁰ atomic force microscopy (AFM),¹¹ infrared spectroscopy (IR),¹² and so on. As a result, it has been known that the hydrated Nafion 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–16 Radical reaction analysis using electron spin resonance has also been conducted.^17–19

To elucidate the hydrated structure of Nafion, many molecular calculations have been carried out to date. Voth et al. analyzed the proton transport mechanism in Nafion at atomistic level using the self-consistent multistate empirical valence bond (SCMS-EVB) approach.^20–22 Choe et al. simulated proton transport using the first principle molecular dynamics.^23–25 Dupuis et al. performed a lot of molecular dynamics (MD) simulations for proton hopping and hydration of Nafion.^26–28 Kawakami and Shigemoto also used MD simulation to verify the diffusion mechanism of proton.²⁹ In addition, using density functional theory (DFT) calculation, water cluster interacting with sulfonic acid group,³⁰ Nafion/Pt interface,³¹ and degradation of Nafion³² have been also reported.

Though Nafion has excellent conductivity as mentioned above, other types of PFSA which have different side chain length from Nafion have been studied.^33–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 Nafion 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.⁴⁰

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 first way is all-atom MD simulations using highly parallelized programs and huge computational resources. There have been various excellent programs such as MODYLAS,⁴¹ GROMACS,⁴² NAMD,⁴³ and LAMMPS.⁴⁴ 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.⁴⁵ Knox and Voth also carried out atomistic MD simulations to probe morphological models of Nafion.⁴⁶ Komarov et al. performed MD simulations with cell size of up to 36 nm including 4 million atoms.⁴⁷

The second way is coarse-grained MD (CGMD) and related methods such as self-consistent mean field method (SC-MFT) and dissipative particle dynamics (DPD).^48–51 These methods are very useful because large scale behaviors in molecular level are obtainable at reasonable computational costs. Many coarse-grained simulations have been performed for PEM. Voth et al. reported a mesoscale study of the proton transport using smoothed particle hydrodynamics (SPH).^52,53 Wescott et al.⁵⁴ and Galperin et al.⁵⁵ used the polymer SC-MFT theory to reproduce the phase structure. Eikerling et al. performed a CGMD-based study of assembly of ionomer.^56,57 DPD has been frequently used for researches of PEM. Morohoshi et al. used DPD to verify the structural dependences of physical properties of PEM^58,59 and investigated the gas permeation by combining dynamic Monte Carlo (MC).⁶⁰ Neimark et al. modeled the proton dissociation and conductivity.^61–63 Ref. 64–66 utilized DPD as well. Hereafter, we focus on DPD simulations.

The DPD method was proposed by Hoogerbrugge and Koelman,^66,67 and was extended to polymer systems by Groot et al.^48–51 However, the evaluation of a set of interaction parameters (χ) among particles have still been of difficult issue. There have been a couple of major routes of χ parameter predictions. The first route is based on the solubility parameter (SP) values. The SP value method was devised by Hildebrand,⁶⁸ and various empirical estimation models such as atomic group contribution models^69–72 were developed. Note that the SP values could be obtained from molecular simulations for single molecules as well.⁷¹ The second route is based on the evaluation of interactions between heterogeneous molecules, due to the contact energy between segments⁷² and the difference in cohesive energies of the aggregated model.⁷³ 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 field (FF) molecular calculation and reproduced the network structure of the water clusters generated by Nafion.⁷⁴ However, since the interaction between Nafion 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 Nafion, based on ab initio molecular orbital (MO) calculations.⁷⁵ Although ab initio evaluations of χ 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–82 This procedure could be considered as an extension of Fan's method⁷² based on the Flory–Huggins theory.^83,84

In this paper, we report a series of DPD-based investigations on the morphology of hydrated Nafion 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–88} since such an effective index should relate with the conductivity. Nafion and SPEEK were compared. The rest of this paper is configured as follows. In Section 2, the χ parameter evaluations and models used for DPD simulations are described in detail. In Section 3 the simulated results are presented and discussed.

2. Simulation details

2.1. χ Parameter evaluation

In the Flory–Huggins lattice theory, the free energy change (ΔG) for a binary system is expressed as^83,84where φ_i and x_i are the volume fraction and the chain length (i = 1, 2 for the two components), respectively. The first and the second terms on the right side of this equation describe the entropy changes, and the third term provides the enthalpy change. The χ parameter is defined aswhere Z is the coordination number of the model lattice and the contact energy ΔE₁₂ 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 ΔE₁₂ 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 ΔE₁₂, based on MC simulations with classical FF set.⁷²

Recently, we have developed a new procedure^76,77 to estimate the χ 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 χ 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 Nafion. According to Yamamoto's previous study,⁷⁴ the basic unit of Nafion chain is divided into three segments of the same size (A: –CF₂–CF₂–CF₂–CF₂–, B: –O–CF₂–C(CF₃)FO–, and C: –CF₂–CF₂–SO₃H). The termini were capped with F for segments A and C and with CF₃ 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 Nafion. The termini were capped with CH₃. 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.

Fig. 1 Molecular structure and segment models of Nafion monomer.

Fig. 2 Molecular structure and segment models of SPEEK monomer.

Fig. 3 Various water structures used for parameter evaluations. (a) Tetramer. (b) Linear dimer. (c) Bifurcated dimer. (d) Cyclic dimer. (e) Monomer. Blue dots indicate hydrogen bonds between water molecules.

The geometries of segments were optimized at the dispersion-corrected B97D (ref. 89) DFT level with the 6-31G(d′,p′) basis set^90,91 by using the GAUSSIAN09 program.⁹² Note that the orbital exponents of polarization functions of 6-31G(d′,p′) were optimized for the respective elements. The generation of geometrical configurations (each number was typically 2000) was performed with the J-OCTA program.⁹³ A number of FMO calculations for all the possible combinations among segments were carried out at FMO-MP2/6-31G(d′) level, where the ABINIT-MP program⁷⁹ was used in parallel executions on several in-house servers with Intel's Xeon processors.

2.2. DPD simulation

DPD is a sort of soft 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–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 motionandwhere, 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 asThe first three forces of the right-hand-side are considered within a certain cutoff radius r_c. The conservative force F^C_ij is a soft repulsion action as followswhere a_ij is the maximum repulsion force between particles i and j, and related definitions 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 χ parameter set. When the reduced density ρ is assumed to be 3, a linear relation with χ_ij is usually set as^48,49

a_ij = a_ii + 3.27χ_ij, (9)

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 modified Verlet algorithm⁴⁸ with the empirical factor λ = 0.65 and time step Δt = 0.05. The cubic cell system was defined as a dimension of L = 30 DPD length unit, corresponding to 23.5 nm in real unit (DPD particle diameter is 0.71 nm).⁵¹ Each simulation involved about 81 000 beads, and the beads were initially packed randomly at the standard mean density ρ = 3.⁴⁸ 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 Nafion. In Fig. 4, the structure (b) corresponds to Nafion 117 (EW = 1100).³³ 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.

Fig. 4 DPD particle models for Nafion polymer and SPEEK. Particle size is approximately 0.71 nm. Particles with Green, light green and yellow colors indicate segments A, B and C, respectively, shown in Fig. 1 for Nafion and in Fig. 2 for SPEEK. Symbolic structure (b) corresponds to Nafion 117 (EW = 1100),³³ and structures (a) and (c) are defined as low and high EW models, respectively.

3. 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 Nafion, 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 Nafion are shown in Table 2. The E_ij values among the hydrophobic segments (A–A, A–B, B–B) are approximately −1 kcal mol⁻¹, whereas the interaction energies between the hydrophilic segments (C–C, C–W, W–W) are about −10 kcal mol⁻¹. On the other hand, the S_ij values between hydrophobic segments are about 0.9, and those between hydrophilic segments are less 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 χ obtained from the results listed in Table 2 using eqn (4), where the values of Yamamoto and Hyodo's preceding study⁷⁴ are listed for comparison. The χ value using FMO is small for A–B (−0.17) and C–W (−4.1) pairs, while it exceeds 20 for A–W and B–W pair. Although the tendency of these parameters is consistent with the results of ref. 74, the absolute values of our χ related to water particles is larger, suggesting the contributions of polarization and charge transfer interactions incorporated by the FMO calculations.

Table 1 Interaction energies (E_ij), and scaling factors (S_ij) of segment C of Nafion shown in Fig. 1 and various water models (W) at 350 K

	Tetramer	Dimer			Monomer
		Bifurcated	Cyclic	Linear
Eij	−10.18	−12.93	−10.88	−15.58	−10.33
S12 (C–W)	0.24	0.23	0.25	0.20	0.20
S21 (W–C)	0.45	0.27	0.44	0.30	0.52

---

Table 2 Interaction energies (E_ij), scaling factors (S_ij), and coordination numbers (Z_ij) for each segment of Nafion and water at 350 K. Note that label W corresponds to linear shape of water dimer for C–W pair and tetramer for the others

		Segment 1
		A			B			C			W
		Eij	Sij	Zij	Eij	Sij	Zij	Eij	Sij	Zij	Eij	Sij	Zij
Seg. 2	A	−0.76	0.89	7.2	−0.87	0.91	8.0	−0.82	0.83	8.1	−1.32	0.84	8.1
	B	−0.87	0.90	7.0	−0.95	0.89	8.0	−0.91	0.84	8.2	−1.38	0.83	7.3
	C	−0.82	0.92	7.6	−0.91	0.90	8.7	−4.26	0.46	8.6	−15.58	0.30	7.4
	W	−1.32	0.83	13.7	−1.38	0.83	12.3	−15.58	0.20	14.2	−9.87	0.51	10.6

---

Table 3 Calculated χ parameters and corresponding repulsion parameters a_ij for each pair of hydrated Nafion

	This work		Previous work^74
	χ	aij	χ
A–B	−0.17	24.44	0.02
A–C	7.51	49.54	3.11
B–C	7.36	49.05	1.37
A–W	25.37	107.91	5.79
B–W	27.86	116.05	4.90
C–W	−4.1	11.60	−2.79

---

The energy values with the water dimer of linear form were used for parameter evaluation for SPEEK as well as Nafion. 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 Nafion. The dispersion interaction caused by the π 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 χ obtained from the above results is shown in Table 5. The χ 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 χ 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 final a_ij values of Nafion and SPEEK for DPD are found in Table S2.†

Table 4 Interaction energies (E_ij), scaling factors (S_ij), and coordination numbers (Z_ij) for each segment of SPEEK and water at 350 K. Note that label W corresponds to linear shape of water dimer for C–W pair and tetramer for the others

		Segment 1
		A			B			C			W
		Eij	Sij	Zij	Eij	Sij	Zij	Eij	Sij	Zij	Eij	Sij	Zij
Seg. 2	A	−1.27	0.82	8.2	−1.27	0.89	11.0	−2.38	0.76	9.6	−0.98	0.88	9.2
	B	−1.27	0.85	7.1	−1.28	0.89	9.2	−2.28	0.85	7.9	−1.25	0.87	8.3
	C	−2.38	0.77	7.5	−2.28	0.78	11.0	−6.38	0.55	8.3	−8.13	0.76	8.2
	W	−0.98	0.84	17.6	−1.25	0.88	24.5	−8.13	0.37	18.1	−9.89	0.50	10.6

---

Table 5 Calculated χ parameters and corresponding repulsion parameters a_ij for each pair of hydrated SPEEK

Pair	χ	aij
A–B	−0.75	22.55
A–C	4.94	41.14
B–C	3.53	36.54
A–W	28.09	116.80
B–W	19.49	88.69
C–W	−3.78	12.65

---

3.2. Results of DPD simulation

Fig. 5 illustrates the time-dependent morphologies of the 20 vol% hydrated Nafion with symbolic structure (b) of Fig. 4. It is visible that small domains of water aggregation develop into large water clusters surrounded by hydrophilic particle of Nafion. 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.

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.

Fig. 6 Time-dependent water density distributions of the case of 20 vol% water content for Nafion with symbolic structure (b) shown in Fig. 4. The regions in which more than 40% of water particles exist are displayed with yellow color.

Next, the structural changes due to the differences in EW of Nafion (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 figure, one can see the behavior that the size of water cluster (lower row) enlarges as the chain (or length of segment A of Nafion) elongates. Fig. 8 illustrates the dependence of morphology on water contents of 10, 20 and 30 vol% for the Nafion, and Fig. 9 does the corresponding results for SPEEK. Comparison between Fig. 8 and 9 indicate that the water clusters formed by Nafion 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 Nafion. A probable reason for this difference is the fact that there are more potential sites for hydrogen bond with water in SPEEK than in Nafion. In short, Nafion could provide better connectivity of water clusters than does SPEEK.⁸⁶

Fig. 7 Structural changes due to the differences in EW of the case of 20 vol% water content for Nafion. Labels (a), (b) and (c) correspond to symbolic structures shown in Fig. 5. Upper row shows morphologies, and lower rows does water density distributions.

Fig. 8 Dependence of morphology on water contents for Nafion (t = 500) with symbolic structure (b) in Fig. 4. Morphologies of 10, 20, and 30 vol% are depicted in upper row. Water density distributions are given drawn in lower row.

Fig. 9 Dependence of morphology on water contents for SPEEK (t = 500). Morphologies of 10, 20, and 30 vol% are depicted in upper row. Water density distributions are given drawn in lower row.

3.3. Analysis of Nafion 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 Nafion. Fig. 10 illustrates the radial distribution functions (RDFs) of water particles for Nafion with symbolic structure (b) (Nafion 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,⁷⁴ 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 shifts 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).

Fig. 10 Radial distribution functions of water particles for Nafion (t = 500). Solid, dashed and dotted lines plot the results of 10, 20, and 30 vol% water contents, respectively.

Fig. 11 Calculated small-angle scattering patterns for Nafion (t = 500). Symbols of triangle, circle and square represent the results of 10, 20 and 30 vol% water contents, respectively.

According to the precedent study⁷⁴ by Yamamoto and Hyodo, the diameter of water cluster was measured at the point where g(r) became 1 after the first peak (see again Fig. 10). The separation between the first 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 Nafion 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.⁶⁰

Fig. 12 Dependences of cluster diameter and space on water contents of 10, 20 and 30 vol% for Nafion. Symbols of triangle, circle and square represent the results of symbolic structures (a), (b) and (c) in Fig. 4.

3.4. Evaluation of diffusion coefficient

Diffusions of protons in the mesoscopic water network has been considered as a crucial indicator of performance for the PEM.⁹⁵ 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

〈(r(t) − r(0))²〉 = 6Dt. (10)

Fig. 13 plots the dependences of diffusion coefficient on water contents for Nafion and SPEEK. For Nafion, 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, reflecting the difference in connectivity of water clusters from Nafion. In other words, Nafion should provide better connections for proton diffusions. The hydration level λ, 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 Nafion 117 (structure (b)). According to experiments of Zawodzinski et al.,⁹⁵ the diffusion coefficients of that region are roughly 0.8 × 10⁻⁶ to 4.5 × 10⁻⁶ cm² s⁻¹. The present DPD-based results have a corresponding range of 0.6 × 10⁻⁶ to 3.2 × 10⁻⁶ cm² s⁻¹, providing reasonable agreement with the experimental data.

Fig. 13 Dependences of diffusion coefficient on water contents for Nafion and SPEEK. Symbols of triangle, circle, and square represent the results of symbolic structures (a), (b), and (c) of Nafion in Fig. 4, respectively. Symbol of diamond corresponds to SPEEK.

3.5. Percolation analysis

So far, several differences due to the connectivity of water clusters between Nafion and SPEEK have been discussed. To directly evaluate the connectivity value, the percolation analysis was performed. The size of water cluster was defined as

R(i,j) ≤ R_C, (11)

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 followswhere 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 first 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 Nafion 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 Nafion 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 λ (or hydration level), the connectivity became slightly higher with an increase in EW. This tendency was in good agreement with Fontanella's experimental observation.⁸⁸ 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.⁸⁶ Since the connection structure of water clusters should be different between Nafion 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 first peak position of the RDF). The corresponding results are plotted in Fig. 15. For Nafion, the plots of connectivity shift toward lower water content, and a degeneracy of (a) and (b) is resolved. The correspondence to the experimental results⁸⁸ is improved by a certain extent. The results for SPEEK are lower shifted, but the final 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 Nafion (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 Nafion and SPEEK.⁸⁶ Certainly, Nafion is better as the PEM material than SPEEK.

Fig. 14 Water connectivity plots for Nafion and SPEEK with respect to water contents using the transient structure of every 100 steps. Symbols of triangle, circle, and square represent the results of symbolic structures (a), (b), and (c) of Nafion in Fig. 4, respectively. Symbol of diamond corresponds to SPEEK.

Fig. 15 Water connectivity plots for Nafion and SPEEK with respect to water contents using the average structure of every 100 steps. Symbols of triangle, circle, and square represent the results of symbolic structures (a), (b), and (c) of Nafion in Fig. 4, respectively. Symbol of diamond corresponds to SPEEK.

4. Conclusion

In the present study, we investigated the mesoscopic structures of Nafion and SPEEK by using DPD simulations. A crucial set of effective interaction (χ) parameters has been determined with a series of non-empirical FMO calculations for all the possible particle–particle interactions.⁷⁶ The results of DPD simulations were analyzed in several ways, illuminating the differences in mesoscopic structures of water clusters between Nafion 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 Nafion is better that by SPEEK. This difference could be considered as a reason why Nafion has been the first 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–98 have been reported. Deng et al. have made a DPD study of sulfonic gemini surfactants.⁹⁹ 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.¹⁰⁰ 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 χ parameters through the FMO calculations has been packaged as a workflow system (FCEWS, FMO-based Chi-parameter Evaluation System).⁷⁷

Conflicts of interest

The authors declare no competing financial interest.

Acknowledgements

The authors would thank to Dr Yuto Komeiji (AIST Japan) for fruitful comments on the manuscript. This work was supported in part by Ministry of Education, Culture, Sports, Science and Technology (MEXT) as a social and scientific priority issue #6 (Accelerated Development of Innovative Clean Energy Systems) to be tackled by using post-K computer.

References

1. K. A. Mauritz and R. B. Moore, State of Understanding of Nafion, Chem. Rev., 2004, 104, 4535–4586.
2. T. D. Gierke, G. E. Munn and F. C. Wilson, The morphology in nafion perfluorinated membrane products, as determined by wide- and small-angle x-ray studies, J. Polym. Sci. Polym. Phys. Ed., 1981, 19, 1687–1704.
3. H.-G. Haubold, T. Vad, H. Jungbluth and P. Hiller, Nano structure of NAFION: a SAXS study, Electrochim. Acta, 2001, 46, 1559–1563.
4. A. Saccà, A. Carbone, R. Pedicini, G. Portale, L. D'Ilario, A. Longo, A. Martorana and E. Passalacqua, Structural and electrochemical investigation on re-cast Nafion membranes for polymer electrolyte fuel cells (PEFCs) application, J. Memb. Sci., 2006, 278, 105–113.
5. V. Barbi, S. S. Funari, R. Gehrke, N. Scharnagl and N. Stribeck, Nanostructure of Nafion membrane material as a function of mechanical load studied by SAXS, Polymer (Guildf)., 2003, 44, 4853–4861.
6. M. Fujimura, T. Hashimoto and H. Kawai, Small-angle x-ray scattering study of perfluorinated ionomer membranes. 2. Models for ionic scattering maximum, Macromolecules, 1982, 15, 136–144.
7. J. A. Dura, V. S. Murthi, M. Hartman, S. K. Satija and C. F. Majkrzak, Multilamellar Interface Structures in Nafion, Macromolecules, 2009, 42, 4769–4774.
8. F. Xu, O. Diat, G. Gebel and A. Morin, Determination of Transverse Water Concentration Profile Through MEA in a Fuel Cell Using Neutron Scattering, J. Electrochem. Soc., 2007, 154, B1389.
9. G. Gebel, Structural evolution of water swollen perfluorosulfonated ionomers from dry membrane to solution, Polymer (Guildf)., 2000, 41, 5829–5838.
10. Z. Porat, J. R. Fryer, M. Huxham and I. Rubinstein, Electron Microscopy Investigation of the Microstructure of Nafion Films, J. Phys. Chem., 1995, 99, 4667–4671.
11. P. J. James, J. A. Elliott, T. J. McMaster, J. M. Newton, A. M. S. Elliott, S. Hanna and M. J. Miles, Hydration of Nafion studied by AFM and X-ray scattering, J. Mater. Sci., 2000, 35, 5111–5119.
12. R. Buzzoni, S. Bordiga, G. Ricchiardi, G. Spoto and A. Zecchina, Interaction of H₂O, CH₃OH, (CH₃)₂O, CH₃CN, and Pyridine with the Superacid Perfluorosulfonic Membrane Nafion: An IR and Raman Study, J. Phys. Chem., 1995, 99, 11937–11951.
13. K. a. Mauritz and C. E. Rogers, A water sorption isotherm model for ionomer membranes with cluster morphologies, Macromolecules, 1985, 18, 483–491.
14. M. Eikerling, A. A. Kornyshev and U. Stimming, Electrophysical Properties of Polymer Electrolyte Membranes: A Random Network Model, J. Phys. Chem. B, 1997, 101, 10807–10820.
15. P. A. Cirkel and T. Okada, A Comparison of Mechanical and Electrical Percolation during the Gelling of Nafion Solutions, Macromolecules, 2000, 33, 4921–4925.
16. W. Y. Hsu and T. D. Gierke, Ion transport and clustering in nafion perfluorinated membranes, J. Memb. Sci., 1983, 13, 307–326.
17. R. Uegaki, Y. Akiyama, S. Tojo, Y. Honda and S. Nishijima, Radical-induced degradation mechanism of perfluorinated polymer electrolyte membrane, J. Power Sources, 2011, 196, 9856–9861.
18. M. K. Kadirov, A. Bosnjakovic and S. Schlick, Membrane-Derived Fluorinated Radicals Detected by Electron Spin Resonance in UV-Irradiated Nafion and Dow Ionomers: Effect of Counterions and H 2 O 2, J. Phys. Chem. B, 2005, 109, 7664–7670.
19. A. Panchenko, H. Dilger, E. Möller, T. Sixt and E. Roduner, In situ EPR investigation of polymer electrolyte membrane degradation in fuel cell applications, J. Power Sources, 2004, 127, 325–330.
20. M. K. Petersen and G. A. Voth, Characterization of the Solvation and Transport of the Hydrated Proton in the Perfluorosulfonic Acid Membrane Nafion, J. Phys. Chem. B, 2006, 110, 18594–18600.
21. M. K. Petersen, F. Wang, N. P. Blake, H. Metiu and G. A. Voth, Excess Proton Solvation and Delocalization in a Hydrophilic Pocket of the Proton Conducting Polymer Membrane Nafion, J. Phys. Chem. B, 2005, 109, 3727–3730.
22. S. Feng and G. A. Voth, Proton Solvation and Transport in Hydrated Nafion, J. Phys. Chem. B, 2011, 115, 5903–5912.
23. Y.-K. Choe, E. Tsuchida and T. Ikeshoji, First-principles molecular dynamics study on aqueous sulfuric acid solutions, J. Chem. Phys., 2007, 126, 154510.
24. Y.-K. Choe, E. Tsuchida, T. Ikeshoji, S. Yamakawa and S. Hyodo, Nature of Water Transport and Electro-Osmosis in Nafion: Insights from First-Principles Molecular Dynamics Simulations under an Electric Field, J. Phys. Chem. B, 2008, 112, 11586–11594.
25. Y.-K. Choe, E. Tsuchida, T. Ikeshoji, S. Yamakawa and S.-A. Hyodo, Nature of proton dynamics in a polymer electrolyte membrane, nafion: a first-principles molecular dynamics study, Phys. Chem. Chem. Phys., 2009, 11, 3892–3899.
26. R. Devanathan, A. Venkatnathan, R. Rousseau, M. Dupuis, T. Frigato, W. Gu and V. Helms, Atomistic simulation of water percolation and proton hopping in Nafion fuel cell membrane, J. Phys. Chem. B, 2010, 114, 13681–13690.
27. R. Devanathan, A. Venkatnathan and M. Dupuis, Atomistic Simulation of Nafion Membrane: I. Effect of Hydration on Membrane Nanostructure, J. Phys. Chem. B, 2007, 111, 8069–8079.
28. A. Venkatnathan, R. Devanathan and M. Dupuis, Atomistic Simulations of Hydrated Nafion and Temperature Effects on Hydronium Ion Mobility, J. Phys. Chem. B, 2007, 111, 7234–7244.
29. T. Kawakami and I. Shigemoto, Molecular dynamics studies on the structures of polymer electrolyte membranes and diffusion mechanism of protons and small molecules, Polymer (Guildf)., 2014, 55, 6309–6319.
30. J. Kabrane and A. J. a Aquino, Electronic Structure and Vibrational Mode Study of Nafion Membrane Interfacial Water Interactions, J. Phys. Chem. A, 2015, 119, 1754–1764.
31. R. Pacheco and A. M. Rocco, in ECS Transactions, ECS, 2009, vol. 25, pp. 1383–1391.
32. T. H. Yu, Y. Sha, W.-G. Liu, B. V Merinov, P. Shirvanian and W. A. Goddard, Mechanism for Degradation of Nafion in PEM Fuel Cells from Quantum Mechanics Calculations, J. Am. Chem. Soc., 2011, 133, 19857–19863.
33. K. D. Kreuer, M. Schuster, B. Obliers, O. Diat, U. Traub, A. Fuchs, U. Klock, S. J. Paddison and J. Maier, Short-side-chain proton conducting perfluorosulfonic acid ionomers: Why they perform better in PEM fuel cells, J. Power Sources, 2008, 178, 499–509.
34. D. Brandell, J. Karo, A. Liivat and J. O. Thomas, Molecular dynamics studies of the Nafion®, Dow® and Aciplex® fuel-cell polymer membrane systems, J. Mol. Model., 2007, 13, 1039–1046.
35. R. Devanathan and M. Dupuis, Insight from molecular modelling: does the polymer side chain length matter for transport properties of perfluorosulfonic acid membranes?, Phys. Chem. Chem. Phys., 2012, 14, 11281.
36. J. Pozuelo, E. Riande, E. Saiz and V. Compañ, Molecular Dynamics Simulations of Proton Conduction in Sulfonated Poly(phenyl sulfone)s, Macromolecules, 2006, 39, 8862–8866.
37. S. Roy, D. Markova, A. Kumar, M. Klapper and F. Müller-Plathe, Morphology of Phosphonic Acid-Functionalized Block Copolymers Studied by Dissipative Particle Dynamics, Macromolecules, 2009, 42, 841–848.
38. S. J. Paddison, The modeling of molecular structure and ion transport in sulfonic acid based ionomer membranes, J. New Mater. Mater. Electrochem. Syst., 2001, 4, 197–207.
39. M. Rikukawa, Possibility and Problems for Polymer Electrolytes Based on Hydrocarbon Polymers, Membrane, 2003, 28, 14–20.
40. J. Miyake, R. Taki, T. Mochizuki, R. Shimizu, R. Akiyama, M. Uchida and K. Miyatake, Design of flexible polyphenylene proton-conducting membrane for next-generation fuel cells, Sci. Adv., 2017, 3, eaao0476.
41. Y. Andoh, S. Okazaki and R. Ueoka, Molecular dynamics study of lipid bilayers modeling the plasma membranes of normal murine thymocytes and leukemic GRSL cells, Biochim. Biophys. Acta - Biomembr., 2013, 1828, 1259–1270.
42. M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess and E. Lindah, Gromacs: High performance molecular simulations through multi-level parallelism from laptops to supercomputers, SoftwareX, 2015, 1–2, 19–25.
43. J. C. Phillips, R. Braun, W. Wang, J. Gumbart, E. Tajkhorshid, E. Villa, C. Chipot, R. D. Skeel, L. Kalé and K. Schulten, Scalable molecular dynamics with NAMD, J. Comput. Chem., 2005, 26, 1781–1802.
44. S. Plimpton, Fast Parallel Algorithms for Short-Range Molecular Dynamics, J. Comput. Phys., 1995, 117, 1–19.
45. A.-T. Kuo, W. Shinoda and S. Okazaki, Molecular Dynamics Study of the Morphology of Hydrated Perfluorosulfonic Acid Polymer Membranes, J. Phys. Chem. C, 2016, 120, 25832–25842.
46. C. K. Knox and G. A. Voth, Probing Selected Morphological Models of Hydrated Nafion Using Large-Scale Molecular Dynamics Simulations, J. Phys. Chem. B, 2010, 114, 3205–3218.
47. P. V. Komarov, P. G. Khalatur and A. R. Khokhlov, Large-scale atomistic and quantum-mechanical simulations of a Nafion membrane: Morphology, proton solvation and charge transport, Beilstein J. Nanotechnol., 2013, 4, 567–587.
48. R. D. Groot and P. B. Warren, Dissipative particle dynamics: Bridging the gap between atomistic and mesoscopic simulation, J. Chem. Phys., 1997, 107, 4423–4435.
49. R. D. Groot and T. J. Madden, Dynamic simulation of diblock copolymer microphase separation, J. Chem. Phys., 1998, 108, 8713–8724.
50. C. M. Wijmans, B. Smit and R. D. Groot, Phase behavior of monomeric mixtures and polymer solutions with soft interaction potentials, J. Chem. Phys., 2001, 114, 7644–7654.
51. R. D. Groot and K. L. Rabone, Mesoscopic Simulation of Cell Membrane Damage, Morphology Change and Rupture by Nonionic Surfactants, Biophys. J., 2001, 81, 725–736.
52. S. Liu, J. Savage and G. A. Voth, Mesoscale study of proton transport in proton exchange membranes: Role of morphology, J. Phys. Chem. C, 2015, 119, 1753–1762.
53. R. Jorn and G. A. Voth, Mesoscale simulation of proton transport in proton exchange membranes, J. Phys. Chem. C, 2012, 116, 10476–10489.
54. J. T. Wescott, Y. Qi, L. Subramanian and T. Weston Capehart, Mesoscale simulation of morphology in hydrated perfluorosulfonic acid membranes, J. Chem. Phys., 2006, 124, 134702.
55. D. Y. Galperin and A. R. Khokhlov, Mesoscopic morphology of proton-conducting polyelectrolyte membranes of nafion?? type: A self-consistent mean field simulation, Macromol. Theory Simulations, 2006, 15, 137–146.
56. K. Malek, M. Eikerling, Q. Wang, Z. Liu, S. Otsuka, K. Akizuki and M. Abe, Nanophase segregation and water dynamics in hydrated Nafion: Molecular modeling and experimental validation, J. Chem. Phys., 2008, 129, 204702.
57. M. Ghelichi, K. Malek and M. H. Eikerling, Ionomer Self-Assembly in Dilute Solution Studied by Coarse-Grained Molecular Dynamics, Macromolecules, 2016, 49, 1479–1489.
58. G. Dorenbos and K. Morohoshi, Percolation thresholds in hydrated amphiphilic polymer membranes, J. Mater. Chem., 2011, 21, 13503.
59. G. Dorenbos and K. Morohoshi, Pore morphologies and diffusion within hydrated polyelectrolyte membranes: Homogeneous vs. heterogeneous and random side chain attachment, J. Chem. Phys., 2013, 138, 064902.
60. G. Dorenbos and K. Morohoshi, Modeling gas permeation through membranes by kinetic Monte Carlo: Applications to H 2, O 2, and N 2 in hydrated Nafion®, J. Chem. Phys., 2011, 134, 044133.
61. M.-T. Lee, A. Vishnyakov and A. V. Neimark, Coarse-grained model of water diffusion and proton conductivity in hydrated polyelectrolyte membrane, J. Chem. Phys., 2016, 144, 014902.
62. M.-T. Lee, A. Vishnyakov and A. V. Neimark, Modeling Proton Dissociation and Transfer Using Dissipative Particle Dynamics Simulation, J. Chem. Theory Comput., 2015, 11, 4395–4403.
63. A. Vishnyakov and A. V. Neimark, Self-Assembly in Nafion Membranes upon Hydration: Water Mobility and Adsorption Isotherms, J. Phys. Chem. B, 2014, 118, 11353–11364.
64. P. V. Komarov, I. N. Veselov, P. P. Chu, P. G. Khalatur and A. R. Khokhlov, Atomistic and mesoscale simulation of polymer electrolyte membranes based on sulfonated poly(ether ether ketone), Chem. Phys. Lett., 2010, 487, 291–296.
65. D. Wu, S. J. Paddison and J. A. Elliott, A comparative study of the hydrated morphologies of perfluorosulfonic acid fuel cell membranes with mesoscopic simulations, Energy Environ. Sci., 2008, 1, 284.
66. P. J. Hoogerbrugge and J. M. V. A. Koelman, Simulating Microscopic Hydrodynamic Phenomena with Dissipative Particle Dynamics, Europhys. Lett., 1992, 19, 155–160.
67. J. M. V. A. Koelman and P. J. Hoogerbrugge, Dynamic Simulations of Hard-Sphere Suspensions Under Steady Shear, EPL, 1993, 21, 363.
68. H. A. Benesi and J. H. Hildebrand, A Spectrophotometric Investigation of the Interaction of Iodine with Aromatic Hydrocarbons, J. Am. Chem. Soc., 1949, 71, 2703–2707.
69. D. W. Van Krevelen and K. Te. Nijenhuis, Properties of polymers: their correlation with chemical structure; their numerical estimation and prediction from additive group contributions, Amsterdam, Elsevier, 4th edn, 1976.
70. M. M. Coleman, C. J. Serman, D. E. Bhagwagar and P. C. Painter, A practical guide to polymer miscibility, Polymer., 1990, 31, 1187–1203.
71. B. Prathab, V. Subramanian and T. M. Aminabhavi, Molecular dynamics simulations to investigate polymer–polymer and polymer–metal oxide interactions, Polymer, 2007, 48, 409–416.
72. C. F. Fan, B. D. Olafson, M. Blanco and S. L. Hsu, Application of molecular simulation to derive phase diagrams of binary mixtures, Macromolecules, 1992, 25, 3667–3676.
73. S. S. Jawalkar, S. G. Adoor, M. Sairam, M. N. Nadagouda and T. M. Aminabhavi, Molecular Modeling on the Binary Blend Compatibility of Poly(vinyl alcohol) and Poly(methyl methacrylate): An Atomistic Simulation and Thermodynamic Approach, J. Phys. Chem. B, 2005, 109, 15611–15620.
74. S. Yamamoto and S. Hyodo, A Computer Simulation Study of the Mesoscopic Structure of the Polyelectrolyte Membrane Nafion, Polym. J., 2003, 35, 519–527.
75. F. Sepehr and S. J. Paddison, Dissipative Particle Dynamics interaction parameters from ab initio calculations, Chem. Phys. Lett., 2016, 645, 20–26.
76. K. Okuwaki, Y. Mochizuki, H. Doi and T. Ozawa, Fragment Molecular Orbital Based Parametrization Procedure for Mesoscopic Structure Prediction of Polymeric Materials, J. Phys. Chem. B, 2018, 122, 338–347.
77. K. Okuwaki, H. Doi and Y. Mochizuki, An automated framework to evaluate effective interactionparameters for dissipative particle dynamics simulations basedon the fragment molecular orbital (FMO) method, J. Comput. Chem. Japan, 2018, 17, 102–109.
78. D. G. Fedorov, T. Nagata and K. Kitaura, Exploring chemistry with the fragment molecular orbital method, Phys. Chem. Chem. Phys., 2012, 14, 7562.
79. S. Tanaka, Y. Mochizuki, Y. Komeiji, Y. Okiyama and K. Fukuzawa, Electron-correlated fragment-molecular-orbital calculations for biomolecular and nano systems, Phys. Chem. Chem. Phys., 2014, 16, 10310–10344.
80. M. S. Gordon, D. G. Fedorov, S. R. Pruitt and L. V. Slipchenko, Fragmentation Methods: A Route to Accurate Calculations on Large Systems, Chem. Rev., 2012, 112, 632–672.
81. D. G. Fedorov and K. Kitaura, The Fragment Molecular Orbial Method: Practical Applications to Large Molecular Systems, CRC Press, Boca Raton, 2009.
82. K. Kitaura, E. Ikeo, T. Asada, T. Nakano and M. Uebayasi, Fragment molecular orbital method: an approximate computational method for large molecules, Chem. Phys. Lett., 1999, 313, 701–706.
83. D. R. Paul, Polymer Blends, Academic Press, New York, 1978.
84. P. J. Flory, Principles of Polymer Chemistry, Comell University Press, NY, 1953.
85. B. R. Matos, E. I. Santiago, J. F. Q. Rey, C. H. Scuracchio, G. L. Mantovani, L. A. Hirano and F. C. Fonseca, dc Proton conductivity at low-frequency in Nafion conductivity spectrum probed by time-resolved SAXS measurements and impedance spectroscopy, J. Polym. Sci. Part B Polym. Phys., 2015, 53, 822–828.
86. X. Wu, X. Wang, G. He and J. Benziger, Differences in water sorption and proton conductivity between Nafion and SPEEK, J. Polym. Sci. Part B Polym. Phys., 2011, 49, 1437–1445.
87. R. Wódzki, A. Narębska and W. K. Nioch, Percolation conductivity in Nafion membranes, J. Appl. Polym. Sci., 1985, 30, 769–780.
88. J. J. Fontanella, C. A. Edmondson, M. C. Wintersgill, Y. Wu and S. G. Greenbaum, High-Pressure Electrical Conductivity and NMR Studies in Variable Equivalent Weight NAFION Membranes, Macromolecules, 1996, 29, 4944–4951.
89. S. Grimme, Semiempirical GGA-type density functional constructed with a long-range dispersion correction, J. Comput. Chem., 2006, 27, 1787–1799.
90. G. A. Petersson, A. Bennett, T. G. Tensfeldt, M. A. Al-Laham, W. A. Shirley and J. Mantzaris, A complete basis set model chemistry. I. The total energies of closed-shell atoms and hydrides of the first-row elements, J. Chem. Phys., 1988, 89, 2193–2218.
91. G. A. Petersson and M. A. Al-Laham, A complete basis set model chemistry. II. Open-shell systems and the total energies of the first-row atoms, J. Chem. Phys., 1991, 94, 6081–6090.
92. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery Jr, J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. Dannenberg, S. Dapprich, A. D. Daniels, O. Farkas, J. B. Foresman, J. V Ortiz, J. Cioslowski and D. J. Fox, Gaussian 09, Revision B.01, Gaussian, Inc., Wallingford CT, 2016.
93. J-OCTA, http://www.j-octa.com/.
94. T. Aoyagi, F. Sawa, T. Shoji, H. Fukunaga, J. Takimoto and M. Doi, A general-purpose coarse-grained molecular dynamics program, Comput. Phys. Commun., 2002, 145, 267–279.
95. T. A. Zawodzinski, M. Neeman, L. O. Sillerud and S. Gottesfeld, Determination of water diffusion coefficients in perfluorosulfonate ionomeric membranes, J. Phys. Chem., 1991, 95, 6040–6044.
96. H. Doi, K. Okuwaki, Y. Mochizuki, T. Ozawa and K. Yasuoka, Dissipative particle dynamics (DPD) simulations with fragment molecular orbital (FMO) based effective parameters for 1-Palmitoyl-2-oleoyl phosphatidyl choline (POPC) membrane, Chem. Phys. Lett., 2017, 684, 427–432.
97. H. Doi, K. Okuwaki, Y. Mochizuki and T. Ozawa, A New Treatment for Water near the Interface between Lipid Membrane and Silica in Dissipative Particle Dynamics Simulation, J. Comput. Chem. Japan, 2017, 16, 28–31.
98. E. Shinsho, K. Okuwaki, H. Doi, Y. Mochizuki, K. Fukuzawa and E. Yonemochi, Formation mechanism of lipid membrane and vesicle using small angle X-ray scattering and dissipative particle dynamics (DPD) method, J. Comput. Chem. Japan DOI:10.2477/jccj.2018-0012.
99. X. Deng, Y. Yang, Y. Ma, X. Sun, G. Zhou, H. Wu and G. Lu, Self-assembled structure of sulfonic gemini surfactant solution, AIP Adv, 2018, 8, 075003.
100. T. Inokuchi, N. Li, K. Morohoshi and N. Arai, Multiscale prediction of functional self-assembled materials using machine learning: high-performance surfactant molecules, Nanoscale, 2018, 10, 16013–16021.

---

Footnotes

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c8ra07428c

‡ Present address: Research Center for Computational Design of Advanced Functional Materials, National Institute of Advanced Industrial Science and Technology (AIST), Central 2, 1-1-1 Umezono, Tsukuba-shi, Ibaraki 305-8568, Japan.

---