Structural and Mechanical Properties of Ionic Di-block Copolymers via a Molecular Dynamics Approach

Polymerized ionic copolymers have recently evolved as a new class of materials to overcome the limited range of mechanical properties of ionic homopolymers. In this paper, we investigate the structural and mechanical properties of charged ionic homopolymers and di-block copolymers, while using coarse-grained molecular dynamics simulation. Tensile and compressive deformation are applied to the homopolymers and copolymers in the glassy state. The effect of charge ratio and loading direction on the stress-strain behavior are studied. It is found that the electrostatic interactions among charged pairs play major roles, as evidenced by increased Young’s modulus and yield strength with charge ratio. Increased charge ratio lead to enhanced stress contribution from both bonding and pairwise (Van der Waals + coulombic) interaction. The increase in the gyration of the radius is observed with increasing charge ratio in homopolymers, yet a reversed tendency is observed in copolymers. Introduced charge pairs leads to an increased randomness in the segmental orientation in copolymers.


Introduction
Polymerized ionic copolymers have recently evolved as a new class of materials to overcome the limited range of mechanical properties of ionic homopolymers, in order to achieve a combined high modulus and conductivity. Well defined nanostructures of long-range order from body-centered-cubic spheres to lamellae can be achieved through the self-assembly of incompatible polymers that are covalently attached to one another. Tunable morphology and domain size can be obtained by varying the chain length, interactions between block polymers, etc. For example, in the case of AB di-block copolymer, with reducing f A , the fraction of A block, lamellar, cylindrical, gyroidal, and spherical ordered phases can gradually form from the initial disordered phase.
Block copolymers (BCP) in solid state can form robust films to accelerate the transport of ions and solvents within the self-assembled nanostructure channels, offering unique physiochemical advantages, such as the high electromechanical stability for specific ions within these channels. As the ionic BCPs combine the properties of ionic liquids and BCPs, there is a particular interest to apply them as electrolyte separator in solid-state lithium ion batteries. In the ionic BCPs, one block serves as the solid polymer electrolyte for the ion-conductive pathway, while the other block provides mechanical support [1,2]. Another important application is in the field of ionic electroactive polymer actuators. Upon the application of an electric field, the positive counterions move towards the negative electrode and negatively charged monomers fixed to the polymer backbones experience attractive force from the positive electrode. Previous studies found that this type of electroactive actuators demonstrate q √ 4πσε 0 .
With q = 1.6e −19 C, σ = 0.5 nm, and /k B = 298 , the estimated q * = 10.6. The number of counterions added to the system also varies in order to keep electrically neutral. A two-step equilibration strategy [19,27,28] is employed to obtain the nanostructured BCPs. First, the polymers and counterions are randomly placed in the simulation box, with an initial number density set as ρ * = 0.8. An energy minimization step is performed to avoid monomer and counterions overlap. Subsequently, a soft-core potential, i.e. dissipative particle dynamics (DPD) type potential, is applied to quickly form nanostructured morphology. This type of soft-core DPD force field has been widely used to study colloidal particles and polymeric fluids [19,29]. This model allows for us to obtain ordered phase with much less computational time when compared to the Kremer-Grest (KG) model, which is usually employed to obtain refined structures. In the DPD model, the soft-core potentials are in the form, as follows: where a ij is the pairwise interaction parameter between polymer beads i and j. r represents distance of a pair of polymer beads and r c is the cut-off distance set as 1.68σ. The interaction strength between all components is set to be a ij = 15.4 k B T except that the interaction between different monomer beads is set as a AB = 50.0 k B T. NPT ensemble is employed where the Noose-Hoover thermostat and Berendsen barostat are used to maintain a constant pressure of P * = 4.5 and temperature of T * = 1.0. The temperature friction coefficient is set as 0.5. The equilibration time is set as 5e5 time steps with a time step size of dt * = 0.02. The choice of soft-core model parameters is to ensure the match between the soft-core and hard-core (KG) model in predicting the densities and other structural characteristics.
Interested readers are referred to [19], where detailed discussion on choosing the model parameters can be found. After the ordered phase di-BCP is obtained, the Kremer-Grest (KG) bead spring model with hard-core potential is applied for a second step equilibration in order to obtain refined structure. The KG model includes a shifted Lennard-Jones (LJ) potential that is used to describe the non-bonded interaction, as follows, r ≤ r c 0 , r > r c (2) where σ is the unit length, is the unit energy of this model, and S(r) is the shifting function to ramp the energy smoothly to zero at the cut-off distance [30]. The interaction strength between all of the components is set as = 1.0 kT, except that the strength between A-B pair of beads is set as = 0.1 kT, and r c = 2.5 σ. Adjacent polymer beads on a chain are bonded via a finite extensible nonlinear elastic (FENE) potential, as follows, where K = 30 σ −2 and the maximum possible length of the bond is R 0 = 1.5σ. The pairwise coulombic interaction is also turned on in the hard-potential model to incorporate the electrostatic interaction between charged polymer beads and counterions into the total potential. Coulombic interaction can be expressed, as follows: where ε 0 is the permittivity of free space, ε r is relative permittivity (dielectric constant), and q i , q j are charge value of components i and j. Dielectric constant, ε r , is set to be 10 to simulate the effects of moderately polarized polymeric materials. The cut-off distance for LJ interaction is set as 2.5. The long-range columbic interaction is solved while using the PPPM technique. Since we are not restricted to a particular type of polymers, reduced unit will be used throughout the paper, so that r * = r/σ, T * = k B T/ , and t * = t/τ, where τ = mσ 2 / 1/2 . The two-step strategy is also applied to generate the initial structure of HPs with the interaction strength a and between all of the components set to be same. After the equilibrated structure has been obtained using the KG model, HPs and BCPs are both quenched down to T * = 0.3 at a cooling rate of ∆T * /∆t * = 1e −7 . At T * = 0.3, which is below the glass transformation temperature, the deformation tests are conducted. Uniaxial tensile and compressive loadings are conducted by changing the box length under a constant engineering strain rate of . ε * = 0.002, whereas the lateral pressures are kept as zero. The stress-strain behaviors of di-BCPs and HPs of varying charge ratio are examined up to a strain of ε = 0.1. Tension and compression simulations are both performed along two coordinate directions for BCPs to account for the effect of anisotropy. It is worth mentioning that the KG model is used for all of the subsequent CGMD simulations for understanding the structural and mechanical properties. At each condition, a minimum of three tests have been conducted to ensure repeatability.

Morphologies and Number Distribution
The simulation boxes of HPs and di-BCPs equilibrated at T * = 1.0, under the anisotropic NPT ensembles are demonstrated in Figures 1 and 2 (top panel), respectively. At first, the pressures in the x, y, and z directions are adjusted in a coupled manner. The anisotropic nature of the system leads to a non-zero averaged pressure in three directions. Therefore, a further equilibration under NPT ensemble is carried out where the pressures in the x, y, and z directions are separately adjusted. The cubic simulation boxes of di-BCPs change their shape due to the anisotropic nature of the systems. The lamellar phase increases in length in the perpendicular to the interface direction during equilibration, accompanied by the shrinkage along the lateral side, similarly to what has been observed by Ryu et al [27]. The ratio of box size in the three directions is approximately 5:3:3 at T * = 1.0 (Figure 2), where the elongated direction is the direction perpendicular to the lamella layer. In contrast, the simulation box of HPs remains approximately cubic during the equilibration (Figure 1). The ratio of box size in the three directions does not significantly change with increasing charge ratio. The counterions locate mainly in the A block adjacent to charged monomers, as seen in Figure 2 (bottom panel). The domain width of A block increases gradually with charge ratio due to the increasing number of counterions that reside in A block.
It is worth mentioning that the system is simulated while using the hard-core potential for sufficiently long time until no significant change in morphology or instability is observed. During the later stage of equilibration, the radius of gyration R g stabilizes with small fluctuation and no tendency of R g increase or decrease with time is observed. The structural and mechanical properties are analyzed after the equilibration and quench process. However, given the limitation in the time scale of MD simulation, it is still possible that the system reaches a local equilibrium state instead of a full equilibrium state. The equilibrium morphology of the charged BCPs has been an active research topic in recent years. It has been theoretically predicted that the morphology diagram is highly asymmetric in charged copolymers. The counter-ions tend to be located in the domains containing the charged ions [31]. Sing and Pryamitsyn et al. discussed that nanostructures can be tuned through the highly asymmetric charge cohesion effect [22,23]. It has been demonstrated that different phases from bcc spheres to hcp cylinders and lamellae can be formed in a charged BCP of the A block fraction f A = 0.25. This is drastically different from the prediction of phase diagram of neutral BCPs. The change in the phase diagram of the charged polymers from that of the uncharged one strongly depends on the strength of the electrostatic cohesion that is indicated by the parameter Γ = e 2 /(8πε 0 ε r kTσ). The parameter Γ expresses the magnitude of the Coulombic potential between two charges of strength e at their closest distance, in units of thermal energy k B T. The phase diagram has not been found to alter much at a low value of Γ.
cubic simulation boxes of di-BCPs change their shape due to the anisotropic nature of the systems. The lamellar phase increases in length in the perpendicular to the interface direction during equilibration, accompanied by the shrinkage along the lateral side, similarly to what has been observed by Ryu et al [27]. The ratio of box size in the three directions is approximately 5:3:3 at * = 1.0 (Figure 2), where the elongated direction is the direction perpendicular to the lamella layer. In contrast, the simulation box of HPs remains approximately cubic during the equilibration (Figure 1). The ratio of box size in the three directions does not significantly change with increasing charge ratio. The counterions locate mainly in the A block adjacent to charged monomers, as seen in Figure  2 (bottom panel). The domain width of A block increases gradually with charge ratio due to the increasing number of counterions that reside in A block.  observed by Ryu et al [27]. The ratio of box size in the three directions is approximately 5:3:3 at * = 1.0 (Figure 2), where the elongated direction is the direction perpendicular to the lamella layer. In contrast, the simulation box of HPs remains approximately cubic during the equilibration (Figure 1). The ratio of box size in the three directions does not significantly change with increasing charge ratio. The counterions locate mainly in the A block adjacent to charged monomers, as seen in Figure  2 (bottom panel). The domain width of A block increases gradually with charge ratio due to the increasing number of counterions that reside in A block. It is worth mentioning that the system is simulated while using the hard-core potential for sufficiently long time until no significant change in morphology or instability is observed. During the later stage of equilibration, the radius of gyration stabilizes with small fluctuation and no tendency of increase or decrease with time is observed. The structural and mechanical properties are analyzed after the equilibration and quench process. However, given the limitation in the time scale of MD simulation, it is still possible that the system reaches a local equilibrium state instead of a full equilibrium state. The equilibrium morphology of the charged BCPs has been an active research topic in recent years. It has been theoretically predicted that the morphology diagram is highly asymmetric in charged copolymers. The counter-ions tend to be located in the domains containing the charged ions [31]. Sing and Pryamitsyn et al. discussed that nanostructures can be tuned through the highly asymmetric charge cohesion effect [22,23]. It has been demonstrated that different phases from bcc spheres to hcp cylinders and lamellae can be formed in a charged BCP of the A block fraction = 0.25. This is drastically different from the prediction of phase diagram of neutral BCPs. The change in the phase diagram of the charged polymers from that of the uncharged one strongly depends on the strength of the electrostatic cohesion that is indicated by the parameter . The parameter Γ expresses the magnitude of the Coulombic potential between The parameter Γ is estimated approximately as 5 when r equals 10. In this case, the phase diagram is shifted vertically instead of leftwards, as predicted by Sing et al. [22]. Therefore, only the segregation strength dependence is varied, but volume fraction dependence is kept in the similar range as that in the neutral copolymers. This possibly explains that the lamellar microstructure remains stable, even with increasing charge density at f A = 0.5 in our study.

Radial Dstribution Function
Radial distribution function g(r) is the function used to describe the molecular-level structure of the system. The g(r) for pairs of beads is given, as follows, where N j r ± ∆r 2 is the number of charged/neutral beads residing in the spherical shell of radius r of bead j, Ω r ± ∆r 2 is the volume of the spherical shell of thickness ∆r, ρ is the averaged density of charged/neutral beads, and N c is the total number of bead j. In Figure 3, the radial distribution function of neutral monomer-neutral monomer, charged monomer-charged monomer, charged monomer-counterion, and counterion-counterion pair are computed for both HPs and di-BCPs with varying charge ratio, denoted as g nm−nm , g cm−cm , g cm−ci , and g ci−ci , respectively. Only the neutral monomers A block in the di-BCPs are counted in the calculation. In g nm−nm , the primary peak locates at around r * = 1.0. Its intensity decreases with charge ratio most apparently at f c = 50% and it shifts slightly to the higher radius. Instead, the peak locating at r * = 1.0 for charged pairs (g cm−cm , g ci−ci ) become stronger. Interestingly, the primary peak of g cm−ci decreases with charge ratio, possibly due to the increased electrostatic screening effect. The charge ratio dependent radial distribution functions g(r) are very similar between HPs and BCPs in terms of the peak location and relative strength. However, g(r) does not converge to unity with increasing r in BCPs as in HPs mainly due to the microphase separation and the charged atoms located only in the A block. varying charge ratio, denoted as , , , and , respectively. Only the neutral monomers A block in the di-BCPs are counted in the calculation. In , the primary peak locates at around * = 1.0. Its intensity decreases with charge ratio most apparently at = 50% and it shifts slightly to the higher radius. Instead, the peak locating at * = 1.0 for charged pairs ( , ) become stronger. Interestingly, the primary peak of decreases with charge ratio, possibly due to the increased electrostatic screening effect. The charge ratio dependent radial distribution functions ( ) are very similar between HPs and BCPs in terms of the peak location and relative strength. However, ( ) does not converge to unity with increasing in BCPs as in HPs mainly due to the microphase separation and the charged atoms located only in the A block.  Figure 4 demonstrates density versus temperature relation of both HPs and di-BCPs with varying charge ratio. Glass transition temperature ( ) can be obtained by the change in the thermal expansion coefficient [32][33][34].

Glass Transition Temperature
of both charged HPs and BCPs gradually increase with charge ratio. These results are in agreement with our previous study in the ionic homopolymer system, where an increase in the glass transition temperature of charged HPs was found as compared to their neutral counterpart [32,33]. At the same charge ratio, the of BCPs is lower than that of HPs due to the lower of B block. As the number of charged ions in HPs is same as that in BCPs when , is half of , , the of 12.5% charged HPs is close to that of 25% charged BCPs, and the of 25% charged HPs is close to that of 50% charged BCPs. Both blocks of the lamellar phase are in their glassy state at * = 0.3.  Figure 4 demonstrates density versus temperature relation of both HPs and di-BCPs with varying charge ratio. Glass transition temperature (T g ) can be obtained by the change in the thermal expansion coefficient [32][33][34]. T g of both charged HPs and BCPs gradually increase with charge ratio. These results are in agreement with our previous study in the ionic homopolymer system, where an increase in the glass transition temperature of charged HPs was found as compared to their neutral counterpart [32,33]. At the same charge ratio, the T g of BCPs is lower than that of HPs due to the lower T g of B block. As the number of charged ions in HPs is same as that in BCPs when f c,A is half of f c,HP , the T g of 12.5% charged HPs is close to that of 25% charged BCPs, and the T g of 25% charged HPs is close to that of 50% charged BCPs. Both blocks of the lamellar phase are in their glassy state at T * = 0.3.   Figure 5 demonstrates the stress-strain behaviors of ionic HPs and BCPs under tensile and compressive deformation at * = 0.3. For the di-BCPs, perpendicular and parallel directions will be used in the following context to refer to the directions perpendicular and parallel to the interfaces. The BCPs deformed along the two different directions do not demonstrate major change in behaviors, and therefore only those in parallel direction are shown in Figure 5. The stress is computed via virial expression.

Mechanical Properties of Copolymers
where N is the number of atoms in the simulation box and N' includes the number of ghost atoms outside the central box. V is the system volume, and are the I and J component of position and force f of atom k, respectively.
The linear elastic regime appears within a small strain range, and then the system enters the nonlinear elastic and plastic regime. A strong dependence on the charge ratio can be observed for both HPs and BCPs. Introducing charged beads also slightly extends the elastic deformation regime.  Figure 5 demonstrates the stress-strain behaviors of ionic HPs and BCPs under tensile and compressive deformation at T * = 0.3. For the di-BCPs, perpendicular and parallel directions will be used in the following context to refer to the directions perpendicular and parallel to the interfaces. The BCPs deformed along the two different directions do not demonstrate major change in behaviors, and therefore only those in parallel direction are shown in Figure 5. The stress σ is computed via virial expression.

Mechanical Properties of Copolymers
where N is the number of atoms in the simulation box and N' includes the number of ghost atoms outside the central box. V is the system volume, r kI and f kJ are the I and J component of position r and force f of atom k, respectively. The linear elastic regime appears within a small strain range, and then the system enters the nonlinear elastic and plastic regime. A strong dependence on the charge ratio can be observed for both HPs and BCPs. Introducing charged beads also slightly extends the elastic deformation regime.
Elastic modulus is obtained by the linear fitting of elastic region (ε < 0.025 ), and yield strength is measured by the 0.2% offset method. Figure 6 plots the extracted values. Under both tension and compression along the two directions, increased mechanical properties (Young's modulus and yield strength) are found with increasing charge ratio in the HPs and BCPs, most obviously from neutral to 12.5% charged in HPs and from neutral to 25% charged in BCPs. In the ionic HPs, the improvement in mechanical properties is not prominent when the charge ratio is further increased from 25% to 50%. As expected, the Young's modulus and yield strength under compression is higher than that under tension. Elastic modulus is obtained by the linear fitting of elastic region ( 0.025), and yield strength is measured by the 0.2% offset method. Figure 6 plots the extracted values. Under both tension and compression along the two directions, increased mechanical properties (Young's modulus and yield strength) are found with increasing charge ratio in the HPs and BCPs, most obviously from neutral to 12.5% charged in HPs and from neutral to 25% charged in BCPs. In the ionic HPs, the improvement in mechanical properties is not prominent when the charge ratio is further increased from 25% to 50%. As expected, the Young's modulus and yield strength under compression is higher than that under tension.  Elastic modulus is obtained by the linear fitting of elastic region ( 0.025), and yield strength is measured by the 0.2% offset method. Figure 6 plots the extracted values. Under both tension and compression along the two directions, increased mechanical properties (Young's modulus and yield strength) are found with increasing charge ratio in the HPs and BCPs, most obviously from neutral to 12.5% charged in HPs and from neutral to 25% charged in BCPs. In the ionic HPs, the improvement in mechanical properties is not prominent when the charge ratio is further increased from 25% to 50%. As expected, the Young's modulus and yield strength under compression is higher than that under tension. Among all of the materials systems studied, the lamellar BCPs with the highest charge ratio ( f c,A = 50%) demonstrates the best mechanical properties, possessing the highest yield stress among all charge ratios and morphologies. The yield strength and Young's modulus reach σ * y = 1.08 and E * = 36.46 under tension, and σ * y =1.14 and E * = 39.59 under compression when deformed in the parallel direction. This leads to the yield strength on the order of 100 MPa and Young's modulus in the order of 1 GPa, assuming that σ = 0.5 nm and /k B = 298 K.
It is worth mentioning that even though no obvious difference is observed when the ionic BCPs are deformed along the perpendicular and parallel directions up to the strain of 0.1, a larger tensile  (Figure 7). When stretched in the perpendicular direction, there is a clear strain hardening, which is hard to observe when stretched in the parallel direction. The anisotropic mechanical properties are very likely related to the chain and bond alignment, as will be demonstrated in the following sections. * = 36.46 under tension, and σ * =1.14 and * = 39.59 under compression when deformed in the parallel direction. This leads to the yield strength on the order of 100 MPa and Young's modulus in the order of 1 GPa, assuming that = 0.5 and ⁄ = 298 K. It is worth mentioning that even though no obvious difference is observed when the ionic BCPs are deformed along the perpendicular and parallel directions up to the strain of 0.1, a larger tensile deformation (ε = 0.3) demonstrates that the loading direction does play a role in the deformation of BCPs (Figure 7). When stretched in the perpendicular direction, there is a clear strain hardening, which is hard to observe when stretched in the parallel direction. The anisotropic mechanical properties are very likely related to the chain and bond alignment, as will be demonstrated in the following sections.

Radius of Gyration during Deformation
To understand the conformation of the polymer molecules of different morphologies during deformation, mean square radius of gyration, 〈 〉 is calculated, as follows, where is the degree of polymerization of the molecule, = ( , , ) is the position of polymer beads, and = ( , , ) is the center of mass, and 〈⋯ 〉 denoting a statistical average over all of the chains. Its perpendicular and parallel components 〈 ∥ 〉 and 〈 〉 are given as follows if the x direction is the parallel direction and the perpendicular direction lies on the y-z plane:

Radius of Gyration during Deformation
To understand the conformation of the polymer molecules of different morphologies during deformation, mean square radius of gyration, R 2 g is calculated, as follows, where N p is the degree of polymerization of the molecule, r i = (x i , y i , z i ) is the position of polymer beads, and r G = (x G , y G , z G ) is the center of mass, and · · · denoting a statistical average over all of the chains. Its perpendicular and parallel components R 2 g and R 2 g⊥ are given as follows if the x direction is the parallel direction and the perpendicular direction lies on the y-z plane: The gyration of radius R g , as well as its components R g and R g⊥ are defined as R g = R 2 g , R g = R 2 g , and R g⊥ = R 2 g⊥ . Here, the calculated R * g of neutral HPs has the value of~3.0, which is a little larger than the estimation from flexible ideal chains given by b 2 N/6 (≈ 2.57), where b is the monomer length and N the degree of polymerization. R * g of ionic HPs gradually increase with charge ratio (Figure 8). Under tensile and compressive deformation, R * g only slightly varies, when its components R * g and R * g⊥ show a near linear relation with applied strain. Ionic BCPs have significantly increased values of R * g as compared with HPs. Its parallel and perpendicular components are no longer comparable but show a major difference in the direction perpendicular and parallel to the interface (Figure 9), which indicates the chains aligning perpendicular to the interface. Interestingly, the introduced charged beads reduce R * g , in contrast to ionic HPs.
Here, the calculated 〈 * 〉 of neutral HPs has the value of ~3.0, which is a little larger than the estimation from flexible ideal chains given by 6 ⁄ (≈ 2.57), where is the monomer length and the degree of polymerization. 〈 * 〉 of ionic HPs gradually increase with charge ratio ( Figure  8). Under tensile and compressive deformation, 〈 * 〉 only slightly varies, when its components 〈 ∥ * 〉 and 〈 * 〉 show a near linear relation with applied strain. Ionic BCPs have significantly increased values of 〈 * 〉 as compared with HPs. Its parallel and perpendicular components are no longer comparable but show a major difference in the direction perpendicular and parallel to the interface (Figure 9), which indicates the chains aligning perpendicular to the interface. Interestingly, the introduced charged beads reduce 〈 * 〉, in contrast to ionic HPs.  When loaded along the perpendicular direction of lamellar interface, the chains tend to be stretched, thus having larger 〈 * 〉. A reduced 〈 * 〉 is expected when compressed along the same direction ( Figure 9). On the contrary when the lamellar domains are stretched and compressed parallel to the interface, a reduced and increased 〈 * 〉 is expected, respectively, as observed from Figure 10. In both cases of BCPs, the components 〈 ∥ * 〉 and 〈 * 〉 show a near linear relation with applied strain, similar as that observed in HPs. estimation from flexible ideal chains given by 6 ⁄ (≈ 2.57), where is the monomer length and the degree of polymerization. 〈 * 〉 of ionic HPs gradually increase with charge ratio ( Figure  8). Under tensile and compressive deformation, 〈 * 〉 only slightly varies, when its components 〈 ∥ * 〉 and 〈 * 〉 show a near linear relation with applied strain. Ionic BCPs have significantly increased values of 〈 * 〉 as compared with HPs. Its parallel and perpendicular components are no longer comparable but show a major difference in the direction perpendicular and parallel to the interface (Figure 9), which indicates the chains aligning perpendicular to the interface. Interestingly, the introduced charged beads reduce 〈 * 〉, in contrast to ionic HPs.  When loaded along the perpendicular direction of lamellar interface, the chains tend to be stretched, thus having larger 〈 * 〉. A reduced 〈 * 〉 is expected when compressed along the same direction ( Figure 9). On the contrary when the lamellar domains are stretched and compressed parallel to the interface, a reduced and increased 〈 * 〉 is expected, respectively, as observed from Figure 10. In both cases of BCPs, the components 〈 ∥ * 〉 and 〈 * 〉 show a near linear relation with applied strain, similar as that observed in HPs. When loaded along the perpendicular direction of lamellar interface, the chains tend to be stretched, thus having larger R * g . A reduced R * g is expected when compressed along the same direction ( Figure 9). On the contrary when the lamellar domains are stretched and compressed parallel to the interface, a reduced and increased R * g is expected, respectively, as observed from Figure 10. In both cases of BCPs, the components R * g and R * g⊥ show a near linear relation with applied strain, similar as that observed in HPs.

Second Order Legendre Polynomials
Bonds that connect the monomer beads can play a major role in the mechanical performance of the polymeric material, second order Legendre polynomials 〈 (cos )〉 is thus calculated, as Figure 10. Variation of radius of gyration R * g and its parallel and perpendicular components R * g and R * g⊥ when deformed in the parallel to interface direction. The parallel direction is chosen as the deformation direction.

Second Order Legendre Polynomials
Bonds that connect the monomer beads can play a major role in the mechanical performance of the polymeric material, second order Legendre polynomials P 2 (cos θ) is thus calculated, as follows, where θ is the angle between the direction of two adjacent monomer beads in the chain and the reference direction, chosen as the deformation direction. The value of P 2 (cos θ ) ranges from −0.5 to 1. When P 2 (cos θ) = 1 and −0.5, it means that the average bond orientation in the system is parallel and perpendicular to the tensile or compressive direction, respectively. As shown in Figure 11, ionic HPs have random distributed bonds initially, which under tension and compression preferentially align along and perpendicular to the loading direction, respectively. In BCPs, the bonds tend to align toward the direction perpendicular to lamellar interface initially, as reflected by a nonzero P 2 (cos θ) , which has a positive value when deformed perpendicular to the interface and negative when deformed in parallel ( Figure 12). P 2 (cos θ) increases under tension and decreases under compression in all phases, which shows that the bonds gradually align to the direction parallel to the tensile direction or perpendicular to the compressive direction, respectively. The absolute change in P 2 (cos θ) value with deformation has been found to be very close between the parallel and perpendicular directions. Yet, the perpendicular loading direction can more easily reorient the chain segments, as observed from the slightly more significant change in the value of P 2 (cos θ) with loading. The change in f c has been observed to shift the values of P 2 (cos θ) closer to zero (Figure 12), which indicates more randomly ordered chain segments associated with the increasing number of counterions, in agreement with the analysis from radius of gyration.

Discussion
Introduced charged beads and counterions have been demonstrated to affect the configurational and mechanical properties of HPs and BCPs. With an increase in charge ratio, the radius of gyration of HPs gradually increase, same as the prediction from previous theoretical study on a single charged polymer chain that shows a stiffened chain due to the repulsive intrachain interaction [35]. In the charged polymer-counterion melt/glass, the types of interaction become more complicated in that the counterions also interact with polymer monomers via Van der Waals force described by the LJ potential. Thus, it is difficult to attribute the swollen polymers chains purely to the electrostatic effect. The charge of both monomers and counterions in the 12.5% charged HPs is set as zero and the system is equilibrated for another 1e7 timesteps until the radius of gyration fluctuates at a constant value to clarify the influence of Van der Waals and electrostatic effects. It was found the 〈 * 〉 increases from around 3.1 to 3.2 when the charge is set to zero, i.e., when only the Van der Waals interaction takes place. This is most likely due to the increased intrachain and interchain attraction through crosslinked charged monomers via counterions, which reduces the gyration of radius. The variation of 〈 〉 is observed to decrease in BCPs with charge ratio, still most likely attributed to the increased intrachain and interchain attraction, which prevents the chains from aligning perpendicular to the interface direction, as evidenced by the increased 〈 〉 with charge ratio (Figure 9).

Discussion
Introduced charged beads and counterions have been demonstrated to affect the configurational and mechanical properties of HPs and BCPs. With an increase in charge ratio, the radius of gyration of HPs gradually increase, same as the prediction from previous theoretical study on a single charged polymer chain that shows a stiffened chain due to the repulsive intrachain interaction [35]. In the charged polymer-counterion melt/glass, the types of interaction become more complicated in that the counterions also interact with polymer monomers via Van der Waals force described by the LJ potential. Thus, it is difficult to attribute the swollen polymers chains purely to the electrostatic effect. The charge of both monomers and counterions in the 12.5% charged HPs is set as zero and the system is equilibrated for another 1e7 timesteps until the radius of gyration fluctuates at a constant value to clarify the influence of Van der Waals and electrostatic effects. It was found the R * g increases from around 3.1 to 3.2 when the charge is set to zero, i.e., when only the Van der Waals interaction takes place. This is most likely due to the increased intrachain and interchain attraction through crosslinked charged monomers via counterions, which reduces the gyration of radius. The variation of R g is observed to decrease in BCPs with charge ratio, still most likely attributed to the increased intrachain and interchain attraction, which prevents the chains from aligning perpendicular to the interface direction, as evidenced by the increased R g⊥ with charge ratio (Figure 9).
The variation of R g and R g ⊥ and P 2 (cos θ) demonstrate a linear dependence with applied strain (Figures 8-10), showing similar dependence at a different charge ratio. To further understand the strengthening origin with increasing charge ratio, energy decomposition into bonded, Van der Waals, and coulombic contribution are conducted and demonstrated in Figure 13. Two neutral systems, 0% charged HP and 0% charged BCP, are compared with two charged systems, 12.5% HP and 25% BCP, during tensile loading. The two charged systems contain the same number of charged ions to enable a direct comparison. The results show that the neutral and charged systems have close values in bonded and Van der Waals energy, respectively. Interestingly, the introduction of charged beads lowers the bonded energy. Stress decomposition reveals that the bonded contribution in all cases is compressive stress, which indicates that bonds are in compression. Therefore, in charged system, the bonds in compressive state are getting closer to the equilibrium length ( Figure 13). Van der Waals energy is also lower in charged systems as compared with that in neutral systems. Coulombic interaction has a lower absolute value than Van der Waals considering only 25% monomers are charged. The Van der Waals energy has a relatively sharper increase at lower strain than at larger strain, which is consistent with the tensile loading of polymeric materials at similar strain range [36], despite some difference in the physical models. Different from Van der Waals energy, the coulombic energy has a close-to-linear relation with applied strain. The total pair energy, i.e., the sum of Van der Waals and coulombic energy, demonstrates the energy that was stored in charged systems is higher than that in neutral systems.
cases is compressive stress, which indicates that bonds are in compression. Therefore, in charged system, the bonds in compressive state are getting closer to the equilibrium length ( Figure 13). Van der Waals energy is also lower in charged systems as compared with that in neutral systems. Coulombic interaction has a lower absolute value than Van der Waals considering only 25% monomers are charged. The Van der Waals energy has a relatively sharper increase at lower strain than at larger strain, which is consistent with the tensile loading of polymeric materials at similar strain range [36], despite some difference in the physical models. Different from Van der Waals energy, the coulombic energy has a close-to-linear relation with applied strain. The total pair energy, i.e., the sum of Van der Waals and coulombic energy, demonstrates the energy that was stored in charged systems is higher than that in neutral systems.  Figure 14. In Equation (6), can be decomposed into bonded and non-bonded (pair) terms, i.e., = , + , , and the corresponding stress contributions are thus obtained. The first term related to kinetic energy is not included in the calculated stress that is shown in Figure 14, and hence only the second term in the curly bracket is considered, Even though the bonded energy plot shows no clear tendency, the stress from bonded contribution has a clear initial increase, followed by slight decrease or flattening. It is worth mentioning that the plots are made relative to the values at zero deformation, and the absolute bonded stress remain at negative values through the loading process. The bonded stress has greater contribution in charged systems than in neutral systems, similar to the case of pair stress. Interestingly, the enhancement of stress from bonded and non-bonded contribution are comparable with increasing charge ratio, and lead to the increased stress values, as demonstrated in the + plot. The increase in bonded and pair stress have also been observed in the 25% charged HP and 50% charged BCP (figures not shown). In addition, it is found that the HPs and BCPs that possess the same number of ions have similar enhancement in mechanical properties ( Figure 14).

The stress decomposition into bonded and non-bonded interactions is further demonstrated in
The final remark is that our study found the anisotropy of charged BCP is not apparent at small deformation. The preliminary study shows that the effect might become apparent at larger deformation ( Figure 7). Therefore, more detailed study of the anisotropy effect will be considered for  Figure 14. In Equation (6), f k can be decomposed into bonded and non-bonded (pair) terms, i.e., f k = f k,bond + f k,pair , and the corresponding stress contributions are thus obtained. The first term related to kinetic energy is not included in the calculated stress that is shown in Figure 14, and hence only the second term in the curly bracket is considered, Polymers 2017, 9, x FOR PEER REVIEW 14 of 16 future work. Of equivalent interest are the viscoelastic behaviors at rubbery state and the effect of different loading mode, e.g., shearing, which will also be considered in our work in the near future. Figure 14. The change of bonded and pair (Van der Waals + coulombic) stress contribution relative to the initial values at zero deformation during tensile loading perpendicular to the interface.

Conclusion
In this study, CGMD simulations are performed to investigate the structural and mechanical behaviors of ionic HPs and di-BCPs while using CGMD simulation. The materials are evaluated in terms of structural, configurational, and mechanical properties, which have demonstrated quite distinct features as a result of the varying charge ratio and loading directions. Some main findings are summarized, as below: 1. A strong dependence of mechanical properties on increasing charge ratio have been observed in ionic HPs up to 25% and in BCPs up to 50% charged A block. Further increasing the charge ratio from 25% to 50% in HPs does not result in the improvement of mechanical properties. 2. Anisotropy in mechanical properties in the lamellar BCPs is not apparent when deformed to a strain level of 0.1 along the directions parallel and perpendicular to the microphase interface. When stretched up to a strain level of 30%, the mechanical responses start to show large discrepancies, which is expected as the weak interface between A and B block orients differently with respect to the two loading directions. 3. Chain segments preferentially align perpendicular to the microphase interface in ionic di-BCPs, in contrast to the random distribution in ionic HPs. The radius of gyration 〈 〉 increases with increasing charge ratio in HPs. However, 〈 〉 shows a reverse tendency with the charge ratio in BCPs. A linear dependence of the parallel and perpendicular components of the radius of gyration and bond orientation on the applied strain is observed under all loading conditions of HPs and BCPs. 4. HPs do not show a significant change in the radius of gyration, whereas di-BCPs have an Even though the bonded energy plot shows no clear tendency, the stress from bonded contribution has a clear initial increase, followed by slight decrease or flattening. It is worth mentioning that the plots are made relative to the values at zero deformation, and the absolute bonded stress remain at negative values through the loading process. The bonded stress has greater contribution in charged systems than in neutral systems, similar to the case of pair stress. Interestingly, the enhancement of stress from bonded and non-bonded contribution are comparable with increasing charge ratio, and lead to the increased stress values, as demonstrated in the σ bond + σ pair plot. The increase in bonded and pair stress have also been observed in the 25% charged HP and 50% charged BCP (figures not shown). In addition, it is found that the HPs and BCPs that possess the same number of ions have similar enhancement in mechanical properties ( Figure 14).
The final remark is that our study found the anisotropy of charged BCP is not apparent at small deformation. The preliminary study shows that the effect might become apparent at larger deformation ( Figure 7). Therefore, more detailed study of the anisotropy effect will be considered for future work. Of equivalent interest are the viscoelastic behaviors at rubbery state and the effect of different loading mode, e.g., shearing, which will also be considered in our work in the near future.

Conclusions
In this study, CGMD simulations are performed to investigate the structural and mechanical behaviors of ionic HPs and di-BCPs while using CGMD simulation. The materials are evaluated in terms of structural, configurational, and mechanical properties, which have demonstrated quite distinct features as a result of the varying charge ratio and loading directions. Some main findings are summarized, as below:

1.
A strong dependence of mechanical properties on increasing charge ratio have been observed in ionic HPs up to 25% and in BCPs up to 50% charged A block. Further increasing the charge ratio from 25% to 50% in HPs does not result in the improvement of mechanical properties.

2.
Anisotropy in mechanical properties in the lamellar BCPs is not apparent when deformed to a strain level of 0.1 along the directions parallel and perpendicular to the microphase interface.
When stretched up to a strain level of 30%, the mechanical responses start to show large discrepancies, which is expected as the weak interface between A and B block orients differently with respect to the two loading directions.

3.
Chain segments preferentially align perpendicular to the microphase interface in ionic di-BCPs, in contrast to the random distribution in ionic HPs. The radius of gyration R g increases with increasing charge ratio in HPs. However, R g shows a reverse tendency with the charge ratio in BCPs. A linear dependence of the parallel and perpendicular components of the radius of gyration and bond orientation on the applied strain is observed under all loading conditions of HPs and BCPs.

4.
HPs do not show a significant change in the radius of gyration, whereas di-BCPs have an increase in R g when stretched perpendicular to the interface or compressed parallel to the interface. BCPs have a decrease in R g when stretched parallel to the interface or compressed perpendicular to the interface.