Ultrafast molecular dynamics approach to quantify structural and transport properties of ion exchange polymer: a case study on perfluorinated sulfonic acid polymer

Abstract

A computationally efficient molecular dynamics (MD) simulation approach for evaluating the transport and structural properties of ion exchange polymers (IEPs) is proposed. Prediction of transport and structural properties of IEPs using MD simulation is beneficial in understanding structure–property relations and to design advanced tailor-made variants of such polymers. The IEP is a complex network of polymer chains with ionic end groups. Hence, computational robustness plays a key role, especially in large simulation cells, in avoiding iterative and often time-consuming process to arrive at definitive solutions in terms of physical properties. A novel and robust approach is presented in general and evaluated for perfluorosulfonic acid (PFSA) polymer structure as a case study. While prior researches have analysed transport and structural properties of such polymers using MD simulation in detail, there is a lack of information on the model standard and equilibration protocol. To this end, the present article compares the proposed algorithm to conventional approaches for structure equilibration and demonstrate that the variation in diffusion coefficients (water and hydronium ions) reduces as the number of chains increases, with significantly reduced errors observed in 14 and 16 chains models, even at elevated hydration. The proposed method to achieve equilibration is ∼200% more efficient than conventional annealing and ∼600% more efficient than the lean method.

---

1. Introduction

Ion exchange polymers (IEPs) find applications in energy devices,^1,2 water desalination units,³ biomedical applications, artificial muscles,^4,5 etc. Understanding structural and transport properties is essential for optimizing the performance of IEPs. The structure–property relationship in IEPs is a fundamental concept that connects the molecular structure of these polymers to their physical properties. Furthermore, improving their stability and durability throughout their operational lifetime is crucial. The complex molecular structure of IEPs consists of a polymeric backbone with a side-chain functional group that attracts and exchanges ions.^6–8 Numerous experimental studies are accessible to comprehend the intricate molecular mechanisms that regulate the functioning of IEPs. However, the limitations of experimental methods have led to the introduction of molecular dynamics (MD) simulation methods, which help in understanding the properties at the atomic level. The essential aspect of MD simulation is its computational efficiency. Studying a large-scale simulation cell requires a lot of time and processing resources. Thus, a reliable computing approach is necessary to overcome this challenge.^9–15

One of the most well-known IEPs is the PFSA polymer, whose commercial name is Nafion®. It was first discovered by DuPont in 1960. Nafion® has gained significant recognition among researchers, emerging as a material of choice for various applications, especially in fuel cells. Its outstanding characteristics, including strong ionic conductivity, excellent chemical and thermal stability, and hydrophilicity, are responsible for its widespread application. These characteristics make Nafion® an excellent choice for use as a proton exchange membrane (PEM) in fuel cells, where it facilitates efficient ion transport while ensuring durability and long-term performance.^16–26 PEMFCs have garnered significant interest because of their notable environmental sustainability, rapid initiation at low temperatures (about 80 °C), and prompt adaptability to power requirements. Consequently, it has become an appealing alternative for portable electronic devices and automotive uses.^27–35

The MD simulation approach encompasses fundamental steps such as modeling molecular structure, selecting appropriate force fields, and implementing equilibration methods to achieve an equilibration state. The equilibration state of a simulation cell holds significant importance in MD simulations. Achieving equilibration to get target density is a computationally intensive step. Most of the literature follows the conventional annealing method to get an equilibrated state and the target density of the simulation cell.^36–42 The annealing method involves the sequential implementation of processes corresponding to the NVT (canonical ensemble) and NPT (isothermal–isobaric ensemble) ensembles within a temperature range of 300 K to 1000 K.^36,43–51 Devanathan et al.^43,45 employed the annealing process to attain an equilibrated state. In their study, the hydrated Nafion® model structure is minimized first using the steepest descent technique, followed by annealing in four sequential phases. The simulation employs an NPT ensemble at a temperature of 300 K. Subsequently, NVT ensembles are utilized within the temperature range of 300 K to 600 K, and then from 600 K to 300 K. These iterative methods are employed until the desired density is achieved, followed by a final equilibration stage using NPT at 300 K. Feng et al.³⁶ equilibrated the structure utilizing the annealing process, with twenty annealing cycles conducted within a temperature range of 300 K to 1000 K. This approach was employed to get the desired density. Yunqi Li et al.⁵² performed GPU-accelerated MD simulations to explore the assembled structures of Nafion® at varying water contents, using an anisotropic coarse-grained model with the Gay–Berne potential. The initial simulation configurations were equilibrated from 3600 K to 298 K over 18.5 ns with a time step of 0.5 fs, using annealing method. Cui et al.⁵³ used atomistic MD simulations to investigate water dynamics in a polyamide (PA) matrix, the bulk phase of reverse osmosis membranes. The PA matrix was adjusted to target densities using a 21-step compression and relaxation scheme under NVT and NPT ensembles. Cha's work⁵⁴ follows the annealing method within a temperature range of 300 to 550 K. The process iterated until the simulated density reached the targeted value of 1.7 g cm⁻³. Jang et al.⁵⁵ initially progressively extended the structure at a temperature range of 300 K to 600 K by 50% of its original volume. At 300 K, the structure is compressed to its original volume. In their study, the target density is achieved with five times repetition of the NPT ensemble. Some studies follow a variant but conventional approach to achieve equilibration, referred to in this article as the ‘lean method’. The lean method encompasses two steps, namely the NVT and NPT ensembles. Specifically, in the NPT ensemble of this method, the simulation system is allowed to run for a longer duration.^16,44 For example, in their MD simulation approach, Savage et al.¹⁶ initiated the equilibration procedure with the NPT ensemble at a temperature of 400 K. Subsequently, the structure was cooled down to 300 K while maintaining the NPT ensemble. Finally, the system is subjected to a final NVT ensemble simulation. Nevertheless, conventionally, the annealing method is widely employed in most literature studies, and the lean method is explored in a few previous studies to achieve an equilibrated state. One of the significant disadvantages of these methods is they are computationally expensive. A robust computational algorithm is essential to minimize the simulation time for large structures. This allows for quicker analysis and optimization of complex systems. Therefore, the primary objective of this article is to ascertain an ultrafast computationally efficient method necessary for achieving an equilibrated state. The proposed method offers significant benefits for achieving an equilibrated state across various polymer configurations and compositions.

On the other hand, many investigations have explored different morphological models of Nafion® chains to better understand the structure and transport properties of hydrated PFSA membranes. In the molecular model of Nafion®, each Nafion® chain, in general, is composed of ten Nafion® monomer units.^{43,45,47,49,54,56,57} Devanathan et al.^43,45 studied the Nafion® four-chain model to reveal the nanostructure of hydrated Nafion® membranes. In their study, the diffusion coefficients are computed as a function of hydration level and the mean residence time of H₂O molecules and H₃O⁺. These diffusion coefficients agreed well with experimental NMR⁵⁸ and QENS⁵⁹ findings. Venkatnathan et al.⁴⁹ also demonstrated the Nafion® four-chain model to analyze the structural and dynamic properties, including density, radial distribution functions, coordination numbers, mean square deviations, and diffusion coefficients, and conclude hydronium ions change the hydration structure around sulfonate groups. Cha⁵⁴ highlighted the Nafion® four-chain model and concluded shorter side chains caused more hydroniums to enter the conductive state, resulting in higher ion conductivity. Zhang et al.⁵⁷ revealed that increasing temperature and hydration levels enhanced the diffusion of water molecules and hydronium ions. For the Nafion® eight-chain model, Mabuchi and Tokumasu⁶⁰ studied reactive MD simulation to analyze the correlation between proton transport and water clustering in polymer electrolyte membranes. The cluster analysis results provide a consistent perspective on the correlation between the connection and confinement of water clusters and the transport of protons. Ban et al.⁶¹ used the Nafion® sixteen-chain model to analyze the properties of hydrated Nafion® membranes, explicitly focusing on gas adsorption, diffusion, and penetration. The experimental findings confirm the molecular model of Nafion® in terms of material, shape, free volume, and water diffusivity. Flottat et al.⁶² simulated the twenty-five-chain model of Nafion® and Aquivion. They found the influence of the chemical composition and hydration level on the nanostructure and dynamics of confined water and concluded that water volume fraction directly influences structural and dynamical properties. Cui et al.⁶³ calculated the free energy associated with hydrogen bond creation with information-theoretic methodology. This analysis reveals that the hydrogen bond network production process in Nafion® becomes thermodynamically favorable with an increase in water content.

The major takeaway from the above literature is that each of these studies follows different structural morphology (i.e., the number of polymer chains is not constant), which perhaps depends on the availability of computational resources. In other words, the adequacy of a number of chains for predicting the structural and transport properties is not explicitly clear from the literature. The computational resources needed to assess the membrane's structural and transport properties are directly proportional to the number of Nafion® chain models employed. However, it would be very convenient if the estimated properties are morphologically and computationally independent. This motivates one to look for a morphological threshold beyond which the structural and transport properties significantly do not vary. Therefore, the objective is also to ascertain the minimum number of Nafion® chains required to appropriately evaluate the structural and transport properties of a PFSA membrane. The structural properties include radial distribution functions (RDF), such as sulphur–water molecules, sulphur–hydronium ions, and coordination numbers, which are used to determine the conformation of the structure. The mean square displacement (MSD) and diffusion coefficient are also employed to study transport properties.

2. Computational methodology

The computational methodology comprehensively describes the morphological model, force field parameters, and equilibration approach.

2.1 Atomic model and simulation box

Fig. S1(a) in the SI, illustrates the molecular arrangement of Nafion®. The present investigation considers ten monomer units in one Nafion® chain.⁴⁹ Therefore, the value of ‘n’ in Fig. S1(a) corresponds to 10. The individual chain consists of a total of 682 atoms. Ten hydronium ions (H₃O⁺) are added to a solitary chain to preserve charge neutrality. Water molecules solvate the Nafion® polymer chain and the hydronium ions. This configuration is used to build further chain models consisting of 2, 4, 6, 8, 10, 12, 14, and 16 chains, with hydronium ions and varying water molecules, as shown in Table 1. Initial configuration is generated by Avogadro software.⁶⁴ Next, these molecules are randomly placed in a cubic box using the Packmol software.⁶⁵ The water content (λ) decides the quantity of water molecules within each cell, where λ is the ratio of molecules of water and hydronium ions to the sulfonic groups.⁵⁴ The molecular composition of the hydronium ion-molecules, water molecules, and Nafion® chain is demonstrated in Fig. S1(b).

Table 1 Details of molecules of water (H₂O) and hydronium ions (H₃O⁺) in the simulation cell

No. of Nafion® chain	No. of hydronium ion	No. of water molecules
		λ = 3.5	λ = 7	λ = 10	λ = 13	λ = 16
2	20	50	120	180	240	300
4	40	100	240	360	480	600
6	60	150	360	540	720	900
8	80	200	480	720	960	1200
10	100	250	600	900	1200	1500
14	140	300	720	1080	1440	1800
16	160	350	840	1260	1680	2100

---

2.2 Force field

The present investigation uses all-atom force fields to conduct classical MD simulations on Nafion®, water, and hydronium ions. Simulations are performed with the modified Dreiding force field.^47,56 The modified Dreiding force field signifies a notable progression in the realm of ab initio force fields, as it allows for precise and simultaneous estimation of various properties in both the gas phase (such as structural, conformational, and vibrational properties) and the condensed phase for a diverse array of molecules. Additionally, this force field represents a pioneering advancement in integrating factors pertaining to organic and inorganic materials.

A modified Dreiding force field's functional form can be categorized into bonded and non-bonded interaction terms. The bonded interaction is composed of a bond-stretching term (E_bond), the angle-bending term (E_angle), and torsion angle term (E_torsion) and non-bond interactions, including the van der Waals (E_vdW) term and the electrostatic term (E_{electrostatic}) as explained in eqn (1). Interatomic interactions pertain to the connections between pairs of atoms that are positioned at a distance of two or more atoms in between or are part of separate molecules.

E_total = E_bond + E_angle + E_torsion + E_vdW + E_{electrostatic} (1)

The eqn (1) is expressed in the functional form as

where K_b, K_θ, and V_n are the bond, angle, and torsional force constants, respectively, r₀ and θ₀ are the equilibrium bond lengths and bond angles, respectively, n is the multiplicity, and φ is the phase angle for torsional parameters. In van der Waals interaction ε is van der Waals well depth (kcal mol⁻¹), σ represents the point at which the potential energy between two particles becomes zero, whereas r denotes the distance separating these interacting particles, and in coulombic interaction, q is the partial charge on atoms as described in eqn (2). The Nafion® chain's partial charges are derived from ref. 54. The partial charges utilised in each of these simulations are displayed in Table S1. All interaction parameters-related information leveraged for the Nafion® membrane is available in SI (Tables S2–S5). Three sites, the F₃C model is used for water molecules,^43,54,66 and the classical hydronium model⁴⁰ is used for hydronium molecules. For the Lennard-Jones interactions (van der Waals), the Lorentz–Berthelot mixing rule⁶⁷ is adopted to describe how various atom types interact with one another eqn (3). The particle–particle particle–mesh method (PPPM)⁶⁸ technique is utilised to handle all coulombic interaction, with a specified cut-off distance of 15 Å.

2.3 Computational procedure

An initial equilibration procedure is conducted on the simulation cell to evaluate the structural and transport properties to get the desired density. After achieving equilibration, the subsequent phase is the production run. In production run stage, to ensure statistical convergence, five independent (unique) MD simulation is performed. Ensemble averaging is carried out over these five independent runs to obtain statistically significant result. This approach helps minimize fluctuations and improve the accuracy of computed properties like RDF, MSD and diffusion coefficient. The trajectories obtained from the production run stage are utilised to analyze structural and transport properties during this stage. The method of achieving an equilibrated state is computationally expensive. As we discussed in the Introduction section, the annealing method and the lean method are the two methods that are frequently employed in past studies to reach an equilibrated state. This study proposes an alternative approach for attaining equilibration. The detailed procedure of each method is outlined below. In the next section, the three methods are compared in terms of structural confirmation, target density, and computing efficiency.

The annealing method involves iteratively applying the NVT (canonical ensemble) and NPT (isothermal–isobaric ensemble) ensemble cyclically to achieve the desired density, as demonstrated in Fig. 1. The number of cycles depends on the simulation box size or the system's initial configuration. In this method, the system undergoes temperature cycling (heating and cooling) along with pressure adjustments to reach the target density. This process can be time-consuming as the system requires sufficient relaxation time for molecular rearrangements, especially when starting far from equilibrium. The simulation follows as (i) minimization → (ii) MD simulation for 100 ps at 300 K in NVT ensemble → (iii) MD simulation for 100 ps at 300 K in NPT ensemble → (iv) the structure is heated from 300 K to 700 K, for 100 ps in the NVT ensemble → (v) then the structure is cooled from 700 K to 300 K, for 100 ps in the NVT ensemble → (vi) final MD simulation for 100 ps at 300 K in NPT ensemble → continuing the iterative process outlined in steps (iv) to (vi) until the target density is achieved. The lean method involves conducting the NVT ensemble for a short duration, followed by the NPT ensemble for a more extended period until the desired density is achieved, as demonstrated in Fig. 1. The simulation follows as (i) minimization → (ii) MD simulation for 100 ps at 300 K in NVT ensemble → (iii) MD simulation for 100 ps at 300 K in NPT ensemble → (iv) MD simulation for 12 ns at 300 K in NPT ensemble to reach target density.

Fig. 1 Combined graphical representation of annealing, lean, and proposed method. (Note that gradient blue arrow indicates time duration of individual methods.)

However, the proposed method reduces overall simulation time by introducing two initial steps: reducing (deform) the box size (deformation magnitude is adaptive and determined based on the initial system size and the target density) and performing energy minimization. First, the gradual deformation of the simulation cell helps bypass much of the ‘relaxation’ phase required in the annealing method, where molecular motion and configuration changes would otherwise take longer to reach equilibrium. Following this, energy minimization is performed, allowing the system to find a more stable, lower-energy configuration that is closer to equilibrium. This approach helps avoid configurations that are too far from the desired state, reducing the need for long equilibration to explore unfavourable regions. Subsequently, the remaining three steps are carried out, after which the entire sequence of five steps is repeated until the system reaches the target density. The step-by-step procedure of the proposed method is briefed in Fig. 1. The initial cells with a density of 0.145 g cm⁻³ and box size (100 × 100 × 100 ų) are first optimized using the steepest descent minimization algorithm in LAMMPS, with the minimize command parameters set as: etol = 10 × 10⁻¹⁰ (stopping tolerance for energy, unitless), ftol = 10 × 10⁻¹⁰ (stopping tolerance for force, in force units), maxiter = 10 000 (maximum number of minimization iterations), and maxeval = 10 000 (maximum number of force/energy evaluations). Following the optimization process, further simulation operations are executed → (i) MD simulation for 100 ps at 300 K in NVT ensemble → (ii) MD simulation for 150 ps at 300 K in NPT ensemble → (iii) gradually reduce the box size in three steps, starting with 80 × 80 × 80 ų in first step → (iv) minimization → (v) the structure is heated from 300 K to 700 K, for 100 ps in the NVT ensemble → (vi) then the structure is cooled from 700 K to 300 K, for 100 ps in the NVT ensemble → (vii) MD simulation for 100 ps at 300 K in NPT ensemble → repeat the process two more times, using a box size of 62 × 62 × 62 ų in the second step and 51 × 51 × 51 ų in the third step, to progressively reach the target density. The final configuration has utilized as the initial configuration for equilibrating procedures. The equilibration process is followed by 500 ps NPT simulation. The final step is the production run for 1 ns in NPT simulation at 300 K. The trajectory obtained from the production run stage is utilized to calculate several parameters, including density, RDF, and MSD.

All MD simulation cases in this work are run in LAMMPS (Large-scale Molecular Massively Parallel Simulator) code from Plimpton at Sandia.⁶⁹ The simulation employs three-dimensional periodic boundary conditions with one femtosecond (fs) time step. The Verlet algorithm is employed to numerically integrate the equations of motion.⁵⁸ Nose–Hoover thermostats is used to maintain the desired temperature, and Nose–Hoover barostats^70,71 is used to maintain the atmospheric pressure during simulation. This work utilises a computational resource consisting of an Intel® Xeon® W-2255 CPU @ 3.70 GHz with 10 cores and 20 logical processors, along with 64 GB of RAM. For the Nafion® ten-chain model, at λ = 16 hydration level, the simulation box has 11 720 atoms, 10 110 bonds, 15 200 angles, and 19 100 dihedrals. The system takes approximately 14 hours to complete a simulation from the initial to the production run stage.

Examining configurations derived from MD simulations plays a vital role in understanding the microstructure of hydrated membranes. One method to explore the hydrated membrane's microstructure is RDF. The RDF is a mathematical representation of the likelihood of locating an atom P at a specific distance r from a reference atom A.⁷² This probability can be determined using eqn (4).

In eqn (4) V denotes the overall system volume, N_P denotes the total number of P atoms in the simulation box, n_P is the number of P atoms located at r distance from particle A in a dr-thick shell.

To quantitatively evaluate the correlation between water molecules and hydronium ions with sulfonic acid groups, the coordination number of these molecules around sulphur atoms calculated using the following eqn (5).⁷²

The transportation of water and hydronium ions plays a vital role in the functioning of a polyelectrolyte membrane. Estimating self-diffusion coefficients in classical MD simulations involves an analysis of MSD.^43,54,66 The mathematical description of MSD follows the eqn (6).

In eqn (6), c_i(t) signifies the positional coordinates of atom i at time t, whereas N denotes the number of atoms undergoing free diffusion. As per Einstein's diffusion law, the self-diffusion coefficient D^48,73 can be determined using eqn (7):

3. Results and discussion

The computational efficiency of conventional and proposed algorithms for MD simulation is tested for Nafion® ten-chain model at λ = 10 and λ = 16 hydration levels. Fig. 2(a) shows the computational time for conventional and proposed methods require to reach equilibration. The simulation box size is 115 × 115 × 115 ų, which is the same for both hydration levels, resulting in a lower initial density for low hydration levels compared to high ones. The initial density of the Nafion® ten-chain model is 0.1450 g cm⁻³ (λ = 10) and 0.1568 g cm⁻³ (λ = 16). The annealing method achieves target density in 3.6 ns and 3.9 ns for λ = 16 and λ = 10, respectively. This method takes more time to achieve equilibration for lower hydration levels due to the reduced number of atoms within the same-sized box. The lean method, on the other hand, takes 6 ns and 12 ns for λ = 16 and λ = 10, respectively. On the other hand, the proposed equilibration method took only 1.8 ns and 2 ns for λ = 16 and λ = 10, respectively.

Fig. 2 (a) Computational time required to achieve equilibrated state for all three methods at hydration level (λ = 10 and 16), density vs. time for (b) annealing (c) lean (d) proposed method (for λ = 10 and 16) at 300 K.

Fig. 2(b–d) show the relationship between density variation as a function of computational time (in ‘ns’). Density is computed once the conventional and proposed methods reach their equilibrated states. The density exhibits fluctuations around the average value of 1.65 g cm⁻³ (for λ = 10) and 1.58 g cm⁻³ (for λ = 16) for all three cases. In the annealing method, Fig. 2(b) illustrates the equilibration of density starting from 2.55 ns (λ = 16) and 2.8 ns (λ = 10) and final equilibration achieved at 3.9 ns (λ = 10) and 3.6 ns (λ = 16). Fig. 2(c) demonstrates the equilibration of density starting from 0.5 ns (λ = 16) and 9.2 ns (λ = 10) and final equilibration achieved at 12 ns (λ = 10) and 6 ns (λ = 16). Finally, in the proposed method, Fig. 2(d) showcases the equilibration of density starting from 1 ns (λ = 10) and 0.9 ns (λ = 16) and final equilibration achieved at 2 ns (λ = 10) and 1.8 ns (λ = 16). The RDF of S–O_w and S–O_h for the Nafion® ten-chain model at λ = 10 and λ = 16 for conventional and proposed methods is calculated. The intensity peak appears at 4 Å in all three cases. The results indicate a high degree of similarity in the structural conformation among all three cases. The detailed description is given in SI Section S2 and RDF plots are given in Fig. S2 and S3. Furthermore, to ensure the robustness of the proposed method, key thermodynamic properties including total energy, pressure, and temperature are continuously monitored throughout the simulation. This allows to confirm a smooth and stable evolution of the system, free from any nonphysical behaviour. The proposed method is adaptable and can be extended to other ion-exchange polymers, including structurally complex systems such as branched and cross-linked systems.

In addition to this the structural and transport properties of Nafion® chain models are investigated, with chain lengths of 2, 4, 6, 8, 10, 12, 14, and 16, considering varying hydration levels (λ = 3.5, 7, 10, 13, and 16). Our analysis aims to comprehensively assess the characteristics of these chain models and evaluate their respective properties. The local interactions, as characterized by RDFs, involving the sulfonate sulphur and water oxygen (S–O_w), as well as the sulfonate sulphur and hydronium oxygen (S–O_h), provide information about the sulphur immediate surroundings. Finally, transport properties such as the MSD and diffusion coefficient of water and hydronium are analysed. The water and proton diffusion coefficients are subsequently determined based on the linear regime of the MSD curves.

Fig. 3 shows the density values corresponding to various hydration levels as reported in the literature and by our model. The validation of the equilibration of the simulation cell is determined by whether the system has achieved the target density or not. As discussed in the previous section, the proposed method took considerably less time to reach equilibration; hence, it is used to equilibrate all Nafion® chain models, henceforth. As the hydration level (λ) increases, the membrane experiences a greater degree of swelling with a reduction in density—the density values observed during the equilibration process range from 1.57 to 1.73 g cm⁻³. Density is evaluated at various hydration levels, including 3.5, 7, 10, 13, and 16. The density values corresponding to different hydration levels are presented in Table 2. The density of our simulation box closely aligns with the reference values found in existing literature.

Fig. 3 The density of the simulated system and the relevant experimental data at different hydration levels, our model (square), Venkatnathan et al. (circle),⁴⁹ Tse et al. (triangle),⁴⁴ Ban et al. (upside-down triangle).⁶¹

Table 2 Density value at different hydration levels (λ)

Hydration level (λ)	Density (g cm^−3)
3.5	1.71
7	1.68
10	1.65
13	1.62
16	1.58

---

The snapshots of the Nafion® ten-chain model at various hydration levels (λ = 3.5, 7, 10, 13, and 16) obtained from the production run step NPT simulation for 1 ns at 300 K temperature using visualization tool OVITO⁷⁴ is shown in Fig. S4. The images reveal that the hydrophilic regions comprise molecules of water, hydronium ions, and sulfonate molecules, while the hydrophobic section contains the Nafion® polymer backbones. A series of RDFs are computed to analyze the relationships between various atomic species inside the hydrated Nafion® systems to comprehensively examine the interplay between sulfonate molecules, hydronium ions, and water molecules.

The RDFs of S–O_w (pair correlation function between sulfonate molecules and molecules of water) for 2, 4, 6, 8, 10, 12, 14, 16 chain model with hydration level (λ = 3.5, 7, 10, 13, and 16) at 300 K is plotted in Fig. 4(a), respectively. The maximum peak occurred at about 4 Å; this indicates that a hydration shell surrounds SO₃⁻ at a distance of 4 Å from the sulphur atom. This is similar to prior investigations.^45,48,49 An increase in water content is associated with a drop in peak intensity, indicating a reduction in the interaction between molecules of water and side chains as hydration levels rise. Observing a notable decrease in peak intensity is intriguing as the water content increases from 3.5 to 7. Nevertheless, as the water content increases from 7 to 16, there is a modest reduction in the peak intensity.

Fig. 4 (a) RDFs of S–O_w (oxygen of water) (i–viii) for various Nafion® chain models 2, 4, 6, 8, 10, 12, 14, 16 respectively for hydration level (λ = 3.5 to 16) at 300 K. The maximum peak observed at 4 Å, indicating the most probable distance between the sulphur atom in SO₃⁻ group and O_w. (b) RDFs of S–O_h (oxygen of hydronium) (i–viii) for various Nafion® chain models 2, 4, 6, 8, 10, 12, 14, 16 respectively for hydration level (λ = 3.5 to 16) at 300 K. The maximum peak observed at 4 Å, indicating the most probable distance between sulphur atoms in SO₃⁻ group and O_h.

The RDFs of S–O_h (pair correlation function between sulfonate molecules and hydronium ion) for 2, 4, 6, 8, 10, 12, 14, and 16 chain model with hydration level (λ = 3.5, 7, 10, 13, and 16) at 300 K are shown in Fig. 4(b), respectively. The peak positions exhibit a high degree of similarity to those observed in S–O_w radial distribution functions, occurring at approximately 4 Å.^45,48,49 This observation provides evidence that the sulfonate groups have an affinity for molecules of water and hydronium ions. Nevertheless, the magnitudes of the S–O_h radial distribution functions are significantly higher than the S–O_w RDFs; it is explained by the powerful interaction between the negatively charged sulfonate molecules and the positively charged hydronium ions.

The RDF for sulphur–water and sulphur–hydronium interactions shows inconsistencies for the Nafion® two to six-chain model. However, for the Nafion® six to sixteen-chain model, the RDF exhibits well-defined and consistent profile. This suggests that a higher number of Nafion® chain model allows for a better representation of the structural arrangement and interactions, leading to more accurate and reliable RDF data.

Fig. 5(a) presents the RDFs for water–water interactions at varying hydration levels (λ = 3.5, 7, 10, 13, 16) for the Nafion® 10 chain model at 300 K. The first peak observed at 3 Å, corresponding to the O_w–O_w distance in the first solvation shell. As the hydration level increases, water molecules aggregate into clusters, exhibiting bulk water behaviour, which promotes the formation of ion transport channels. Similarly, Fig. 5(b) shows the RDFs between hydronium ion pairs (O_h–O_h) at 300 K for different hydration levels (λ = 3.5 to 16). The first peak in an RDF plot represents the most probable distance between two hydronium ions. In this case, the strong first peak (highest in λ = 3.5) suggests a significant amount of close-range interaction between hydronium ions at this hydration level. The second peak represents the next coordination shell, which shows hydronium ions that are farther apart than those found in the first coordination shell, indicating medium-range interactions. The distance between the first and second peaks can give insights into the structural organization of ions at intermediate distances. At lower hydration levels (λ = 3.5), the hydronium ions are more clustered, with a prominent first RDF peak indicating strong ionic interactions due to limited water molecules available for solvation. At higher hydration levels increased (λ = 7 to 16) the system becomes more hydrated, causing the ions to spread out, leading to lower RDF peaks and weaker interactions as the ions are better solvated by water molecules. The result is agreed from the previous study.^45,49

Fig. 5 (a) RDFs of O_w–O_w (oxygen of water), the first peak observed at 3 Å and (b) RDFs of O_h–O_h (oxygen of hydronium), across hydration levels (λ = 3.5 to 16) for the Nafion® 10 chain model at 300 K.

Further, coordination numbers (CNs) are a quantitative assessment of the correlation between various molecules calculated by integrating the area under the RDF.⁴⁹ Table 3 displays the CN of molecules of water and hydronium ions surrounding the sulphur atoms for the Nafion® ten-chain model. Similar outcomes are noted in previous studies.^45,48 The CN of water molecules and hydronium ions for the remaining Nafion® chain models are included in SI in Tables S6 and S7, respectively. The observation reveals the CN of water molecules increases as the water content (λ) increases, whereas the CN of hydronium ions decreases with the increase in water content. Water molecule's solvation of sulfonate molecules and hydronium ions is the underlying cause of this phenomenon. The affinity between the sulfonate molecules and the hydronium ions appears to be reduced when more water molecules compete to interact with the sulfonate molecules. The collective CNs of molecules of water and hydronium ions exhibit an upward trend as the water content increases at a constant temperature. As the water content increases, there is a corresponding increase in the size of water clusters within the membrane, leading to a linked transmission channel. This transmission channel helps in the transportation of hydronium ions.

Table 3 Coordination number of Nafion® ten-chain model for water molecules using a cut-off of 4.725 Å and hydronium ion using a cut-off of 4.275 Å at 300 K around the sulphur atoms

Hydration level (λ)	H2O	H3O^+
3.5	3.31	1.92
7	5.70	0.950
10	6.61	0.648
13	7.09	0.488
16	7.26	0.406

---

Fig. 6(a) and (b) show the MSDs of water molecules and hydronium ions computed for several Nafion® chain models during the production run stage of NPT simulation of 1 ns, considering different hydration levels at a 300 K temperature. The final MSD curves represent the average of five independent simulations. The diffusion coefficients are determined by analysing the slope of the MSD graphs for hydronium ions and water molecules, as presented in Tables S8 and S9, respectively. The molecular movement of water is significantly greater than that of hydronium ions. This is the reason molecules of water exhibit a considerably greater MSD compared to hydronium ions. The water diffusion coefficient exhibits strong concordance with experimental measurements of the Nafion® membrane using nuclear magnetic resonance (NMR)⁵⁸ and QENS.⁵⁹

Fig. 6 (a) Mean square displacement (i–viii) of water molecules for various Nafion® chain models 2, 4, 6, 8, 10, 12, 14, 16 respectively for hydration level (λ = 3.5 to 16) at 300 K. (b) Mean square displacement (i–viii) of hydronium ion for various Nafion® chain model 2, 4, 6, 8, 10, 12, 14, 16 respectively for hydration level (λ = 3.5 to 16) at 300 K.

The diffusion coefficient of water (D_w) of a ten-chain model at 300 K for λ = 7 is 0.28782 × 10⁻⁵ cm² s⁻¹, λ = 10 is 0.64182 × 10⁻⁵ cm² s⁻¹, and λ = 13 is 0.77155 × 10⁻⁵ cm² s⁻¹. These values agree with experimental NMR⁵⁸ values for a similar water content at λ = 6 is 0.37297 × 10⁻⁵ cm² s⁻¹, λ = 9 is 0.45135 × 10⁻⁵ cm² s⁻¹ and λ = 14 is 0.58108 × 10⁻⁵ cm² s⁻¹. Our values also agree with simulation results,⁴⁸ which reports water diffusion coefficient at λ = 9 is 0.29245 × 10⁻⁵ cm² s⁻¹, λ = 12 is 0.49057 × 10⁻⁵ cm² s⁻¹ and λ = 15 is 0.66981 × 10⁻⁵ cm² s⁻¹. The error bar plots for the water diffusion and hydronium diffusion coefficients are presented in Fig. S5 and S6, respectively. The Nafion® ten-chain model and experimental values of D_w are shown and compared in Fig. 7(a). Furthermore, the diffusion coefficient value for the Nafion® four-chain model is consistent with the findings from previous investigation.⁴³

Fig. 7 (a) Diffusion coefficient of water molecules and (b) diffusion coefficient of hydronium ion in relation to hydration level as determined by our simulation (rhombus), an NMR⁵⁸ experiment (circle), QENS⁵⁹ experiment (triangle), Devanathan et al. (square),⁴³ Tsung et al. (hexagon)⁴⁸ respectively.

The hydronium ion diffusion coefficient (D_h) of a ten-chain model at 300 K for λ = 7 is 0.04591 × 10⁻⁵ cm² s⁻¹, λ = 10 is 0.10066 × 10⁻⁵ cm² s⁻¹, λ = 13 is 0.21443 × 10⁻⁵ cm² s⁻¹ and λ = 16 is 0.3067 × 10⁻⁵ cm² s⁻¹. Our values agree with simulation results for Nafion® two hundred-chain model, which reported hydronium ion diffusion coefficient at λ = 6 is 0.01765 × 10⁻⁵ cm² s⁻¹, λ = 11 is 0.12941 × 10⁻⁵ cm² s⁻¹, and λ = 15 is 0.0.19118 × 10⁻⁵ cm² s⁻¹. These values are less as compared to experimental NMR values for a similar water content at λ = 6 is 0.12647 × 10⁻⁵ cm² s⁻¹, λ = 10 is 0.33529 × 10⁻⁵ cm² s⁻¹, and λ = 11 is 0.39412 × 10⁻⁵ cm² s⁻¹. Fig. 7(b) shows the diffusion coefficient of hydronium ion for the Nafion® ten-chain model of different studies and experimental studies such as NMR⁵⁸ & QENS⁵⁹ to compare the results. The diffusion coefficient of the hydronium ion is lower than the experimental values obtained for Nafion® membranes. The diffusion coefficient obtained in this work exclusively accounts for the straightforward diffusion of hydronium ions without incorporating the proton hopping phenomenon elucidated by the Grotthuss mechanism. However, our result agrees with previous simulation studies.^45,48

The movement of water molecules is six times higher than the hydronium ions at low water content (λ = 3.5) at 300 K. This is due to the intense electrostatic interaction between the sulfonic acid group (SO₃⁻) and hydronium ions. But, at higher water content (λ = 16) the movement of water molecules is only three times quicker than the movement of hydronium ions. This phenomenon demonstrates the water molecules enhanced ability to dissolve and surround the hydronium ions as the water content increases.

Fig. 8 illustrates the hydronium ion's diffusion coefficient (D_h) with varying hydration levels (λ = 3.5 to 16) for different Nafion® chain models 2 to 16, with corresponding error bars indicating the deviation from the mean value. As the number of chains increases, the error associated with the diffusion coefficient calculation progressively decreases. This trend suggests that the higher Nafion® chain model provides more consistent and reliable results, reducing the uncertainty in the diffusion coefficient. Specifically, beyond the Nafion® 14 chain model, the error margin decreases significantly, indicating that higher Nafion® chain models are more suitable for accurately determining transport properties.

Fig. 8 The diffusion coefficient (D_h) of hydronium ions for various Nafion® chain models (2 to 16) at different hydration levels (λ = 3.5, 7, 10, 13, and 16). The error bars represent the deviation across five independent simulations.

4. Conclusion

An ultrafast, computationally efficient MD simulation is performed on an IEP-based PFSA polymer known as Nafion®. A robust proposed approach is compared with the conventional approach in terms of computational efficiency. Nafion® ten-chain model is used to achieve equilibration for two hydration cases (λ = 10 and λ = 16). The proposed method speeds up achieving the target density by introducing two key steps: gradual simulation box deformation and energy minimization. These steps significantly reduce the prolonged relaxation times typically required in the annealing method, which relies on temperature cycling and pressure adjustments for equilibrium. Deforming the box size bypasses much of the initial relaxation, while energy minimization ensures a more stable, lower-energy configuration. This approach minimizes the need for longer equilibration, ensuring a quicker convergence to the desired density without compromising the system's accuracy or molecular configuration. The RDF analysis of sulphur–water (S–O_w) and sulphur–hydronium (S–O_h) reveals a consistent peak at 4 Å across all cases. This finding shows that the structural conformation is constant in all instances. Notably, the proposed method is six times faster than the lean method and two times faster than the annealing method. Thus, the proposed method is an ultrafast approach to achieve an equilibrated structure.

In addition, the proposed method is employed to analyze the structural and transport characteristics of PFSA membranes. These properties include density, RDF, CN, MSD, and diffusion coefficient. To evaluate the outcomes of the several Nafion® chain models the RDF and diffusion coefficient properties are analyzed. The observed diffusion coefficient for water closely matched experimental studies using NMR and QENS techniques. Conversely, the diffusion coefficient of hydronium ions is significantly lower, ranging from three to six times smaller than that of water.

Furthermore, the deviation in the diffusion coefficient of hydronium ion across five independent simulations decreases as the number of Nafion® chain models increases. For the Nafion® 14 chains model, the error becomes minimal, suggesting that the system achieves a more stable and reliable representation of ion diffusion beyond this point. This emphasizes the importance of simulating higher Nafion® chain models to enhance the reliability and accuracy of structural and transport properties of PFSA membranes.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included as part of the SI. The SI includes Dreiding force field parameter information, RDF, CN, and diffusion coefficient for water molecules and hydronium ions. See DOI: https://doi.org/10.1039/d5ra05434f.

Acknowledgements

This work is supported by DST-SERB Core Research Grant (File No.: CRG/2022/008594), DST-INSPIRE Faculty Award (File No.: DST/INSPIRE/04/2016/000735), DST-SERB Early Career Research Award (File No.: ECR/2018/002635) Schemes by the Department of Science and Technology, Government of India and CPRI-RSoP scheme of Ministry of Power (MoP) (File No.: RSoP/21-26/GD/17), Government of India. Further, the authors acknowledge the use of PARAM UTKARSH HPC facility of the National Supercomputing Mission project and thank Centre for Cyber Physical Systems (CCPS), National Institute of Technology Karnataka, Surathkal for the financial support for the same.

References

1. X. Li, H. Zhang, Z. Mai, H. Zhang and I. Vankelecom, Energy Environ. Sci., 2011, 4, 1147–1160.
2. Y. Liu, Z. Cai, L. Tan and L. Li, Energy Environ. Sci., 2012, 5, 9007–9013.
3. S. J. Kim, S. H. Ko, K. H. Kang and J. Han, Nat. Nanotechnol., 2010, 5, 297–301.
4. V. H. Nguyen, J. Kim, R. Tabassian, M. Kotal, K. Jun, J. H. Oh, J. M. Son, M. T. Manzoor, K. J. Kim and I. K. Oh, Advanced Science, 2019, 6, 1–7.
5. E. Polukhov, D. Vallicotti and M. A. Keip, Comput. Methods Appl. Mech. Eng., 2018, 337, 165–197.
6. C. H. Park, C. H. Lee, M. D. Guiver and Y. M. Lee, Prog. Polym. Sci., 2011, 36, 1443–1498.
7. T. Xu, J. Membr. Sci., 2005, 263, 1–29.
8. P. K. Koorata, J. Electrochem. Soc., 2022, 169, 104505.
9. K. A. Mauritz and R. B. Moore, Chem. Rev., 2004, 104, 4535–4586.
10. S. Wang, L. Y. Zhao and Y. Liu, Comput. Methods Appl. Mech. Eng., 2022, 392, 114644.
11. M. Griebel and J. Hamaekers, Comput. Methods Appl. Mech. Eng., 2004, 193, 1773–1788.
12. Q. Berrod, S. Lyonnard, A. Guillermo, J. Ollivier, B. Frick, A. Manseri, B. Améduri and G. Gébel, Macromolecules, 2015, 48, 6166–6176.
13. A. Tarokh, K. Karan and S. Ponnurangam, Macromolecules, 2020, 53, 288–301.
14. U. Shinde and P. K. Koorata, J. Power Sources, 2022, 527, 231179.
15. K. K. Poornesh and C. Cho, Fuel Cells, 2015, 15, 196–203.
16. J. Savage, Y.-L. S. Tse and G. A. Voth, J. Phys. Chem. C, 2014, 118, 17436–17445.
17. T. Shimoaka, C. Wakai, T. Sakabe, S. Yamazaki and T. Hasegawa, Phys. Chem. Chem. Phys., 2015, 17, 8843–8849.
18. K. Zhang, W. Yu, X. Ge, L. Wu and T. Xu, J. Membr. Sci., 2022, 661, 120922.
19. Y. Wang, D. Y. C. Leung, J. Xuan and H. Wang, Renewable Sustainable Energy Rev., 2016, 65, 961–977.
20. A. Parekh, Front. Energy Res., 2022, 10, 1–13.
21. N. Cele and S. S. Ray, Macromol. Mater. Eng., 2009, 294, 719–738.
22. N. J. Economou, J. R. O'Dea, T. B. McConnaughy and S. K. Buratto, RSC Adv., 2013, 3, 19525–19532.
23. H. Gao and K. Lian, RSC Adv., 2014, 4, 33091–33113.
24. E. Planes, J. Peet, J. B. Brubach, L. Porcar, G. De Moor, C. Iojoiu and S. Lyonnard, Energy Adv., 2024, 3, 1594–1611.
25. M. R. Asghar, W. Zhang, H. Su, J. Zhang, H. Liu, L. Xing, X. Yan and Q. Xu, Energy Adv., 2025, 4, 185.
26. P. K. Koorata and S. D. Bhat, Int. J. Hydrogen Energy, 2021, 46, 5570–5579.
27. R. Ding, S. Zhang, Y. Chen, Z. Rui, K. Hua, Y. Wu, X. Li, X. Duan, X. Wang, J. Li and J. Liu, Energy and AI, 2022, 9, 100170.
28. P. Padavu, P. K. Koorata, S. Kattimani and D. N. Gaonkar, Int. J. Hydrogen Energy, 2024, 84, 435–446.
29. Y. Wang, B. Seo, B. Wang, N. Zamel, K. Jiao and X. C. Adroher, Energy and AI, 2020, 1, 100014.
30. V. Kumar, P. K. Koorata, U. Shinde, P. Padavu and S. C. George, Polym. Degrad. Stab., 2022, 205, 110151.
31. Y. Wang, K. S. Chen, J. Mishler, S. C. Cho and X. C. Adroher, Appl. Energy, 2011, 88, 981–1007.
32. U. Shinde, P. K. Koorata and P. Padavu, Int. J. Hydrogen Energy, 2024, 51, 1497–1511.
33. A. P. Sunda and A. Venkatnathan, J. Comput. Chem., 2011, 32, 3319–3328.
34. P. Padavu, P. K. Koorata and S. Kattimani, Int. J. Heat Mass Transfer, 2023, 214, 124460.
35. X. Ren, E. Gobrogge and F. L. Beyer, J. Membr. Sci., 2021, 637, 119645.
36. C. Feng, Y. Li, K. Qu, Z. Zhang and P. He, RSC Adv., 2019, 9, 9594–9603.
37. S. Kundu, L. C. Simon, M. Fowler and S. Grot, Polymer, 2005, 46, 11707–11715.
38. K. Abu-Hakmeh, P. Sood, B. Jae Chun, J. Il Choi and S. S. Jang, Mater. Perform. Charact., 2015, 4, 131–147.
39. A.-T. Kuo, A. Tanaka, J. Irisawa, W. Shinoda and S. Okazaki, J. Phys. Chem. C, 2017, 121, 21374–21382.
40. J. Wang and T. Hou, J. Comput. Chem., 2011, 32, 3505–3519.
41. S. Chaiwongwattana, M. Phonyiem, V. Vchirawongkwin, S. Prueksaaroon and K. Sagarik, J. Comput. Chem., 2012, 33, 175–188.
42. E. Krieger and G. Vriend, J. Comput. Chem., 2015, 36, 996–1007.
43. R. Devanathan, A. Venkatnathan and M. Dupuis, J. Phys. Chem. B, 2007, 111, 13006–13013.
44. Y.-L. S. Tse, A. M. Herring, K. Kim and G. A. Voth, J. Phys. Chem. C, 2013, 117, 8079–8091.
45. R. Devanathan, A. Venkatnathan and M. Dupuis, J. Phys. Chem. B, 2007, 111, 8069–8079.
46. T. Ohkubo, K. Kidena, N. Takimoto and A. Ohira, J. Mol. Model, 2012, 18, 533–540.
47. I. H. Hristov, S. J. Paddison and R. Paul, J. Phys. Chem. B, 2008, 112, 2937–2949.
48. A.-T. Kuo, W. Shinoda and S. Okazaki, J. Phys. Chem. C, 2016, 120, 25832–25842.
49. A. Venkatnathan, R. Devanathan and M. Dupuis, J. Phys. Chem. B, 2007, 111, 7234–7244.
50. S. Akbari, M. T. H. Mosavian, F. Moosavi and A. Ahmadpour, RSC Adv., 2017, 7, 44537–44546.
51. J. Wang, L. Fan, Q. Du and K. Jiao, Energy Adv., 2022, 1, 269–276.
52. W. Chen, F. Cui, L. Liu and Y. Li, J. Phys. Chem. B, 2017, 121, 9718–9724.
53. F. Cui, W. Chen, X. Kong, L. Liu, C. Shi and Y. Li, J. Phys. Chem. B, 2019, 123, 3086–3095.
54. J. Cha, Sci. Rep., 2020, 10, 22014.
55. S. S. Jang, V. Molinero, T. Çaǧın and W. A. Goddard, J. Phys. Chem. B, 2004, 108, 3149–3157.
56. T. Mabuchi and T. Tokumasu, J. Chem. Phys., 2014, 141, 104904.
57. G. Zhang, G. Yang, S. Li, Q. Shen, H. Wang, Z. Li, Y. Zhou and W. Ye, Membranes, 2021, 11, 695.
58. S. Ochi, O. Kamishima, J. Mizusaki and J. Kawamura, Solid State Ionics, 2009, 180, 580–584.
59. J.-C. Perrin, S. Lyonnard and F. Volino, J. Phys. Chem. C, 2007, 111, 3393–3404.
60. T. Mabuchi and T. Tokumasu, J. Phys. Chem. B, 2018, 122, 5922–5932.
61. S. Ban, C. Huang, X.-Z. Yuan and H. Wang, J. Phys. Chem. B, 2011, 115, 11352–11358.
62. T. Flottat, B. Latour, F. Goujon, P. Hauret and P. Malfreyt, Int. J. Hydrogen Energy, 2023, 48, 33283–33296.
63. R. Cui, S. Li, C. Yu and Y. Zhou, Macromolecules, 2023, 56, 1688–1703.
64. M. D. Hanwell, D. E. Curtis, D. C. Lonie, T. Vandermeersch, E. Zurek and G. R. Hutchison, J. Cheminform., 2012, 4, 17.
65. L. Martinez, R. Andrade, E. G. Birgin and J. M. Martínez, J. Comput. Chem., 2009, 30, 2157–2164.
66. S. S. Awulachew and K. N. Nigussa, Chem. Phys. Lett., 2023, 816, 140380.
67. M. K. Petersen, F. Wang, N. P. Blake, H. Metiu and G. A. Voth, J. Phys. Chem. B, 2005, 109, 3727–3730.
68. D. E. Breen, D. H. House and M. J. Wozny, Text. Res. J., 1994, 64, 663–685.
69. A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in ’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott and S. J. Plimpton, Comput. Phys. Commun., 2022, 271, 108171.
70. W. G. Hoover, Phys. Rev. A: At., Mol., Opt. Phys., 1985, 31, 1695–1697.
71. W. G. Hoover, Phys. Rev. A: At., Mol., Opt. Phys., 1986, 34, 2499–2500.
72. Z. Rao, C. Zheng and F. Geng, Comput. Mater. Sci., 2018, 142, 122–128.
73. F. Müller-Plathe, Acta Polym., 1994, 45, 259–293.
74. A. Stukowski, Modell. Simul. Mater. Sci. Eng., 2010, 18, 015012.

---