Dissipative Particle Dynamics Simulations of a Protein-Directed Self-Assembly of Nanoparticles.

Design and fabrication of multifunctional porous structures play key roles in the development of high-performance energy storage devices. Our experiments demonstrated that nanostructured porous components, such as electrodes and interlayers, generated from the protein-directed self-assembly of nanoparticles can significantly improve the battery performances. The protein-directed assembly of nanoparticles in solution is a complex process involving the complicated interactions among proteins, particles, and solvent molecules. In this paper, we investigate the effects of coating proteins and specific solvent environments on the assembled porous structures. Comprehensive dissipative particle dynamics (DPD) simulations have been implemented to explore the molecular interactions and uncover the fundamental mechanisms in a gelatin-directed self-assembly of carbon black particles under different solvent conditions. Our simulations show that compact triple-strand "rod-like" structures are formed in water while loose curved "sheet-like" structures are formed in an acetic acid/water mixture. The structural difference is mainly due to the redistribution of the charges on the gelatin side chains under specific acid-solvent conditions. The strong and flexible "sheet-like" structures lead to a homogenous porous structure with high porosity and with large functionalized surfaces. Our simulations results can reasonably explain the experimental observations; this work demonstrates the great potential of DPD as a powerful tool in guiding future experimental design and optimization.

Introduction

The advanced energy storage devices (ESDs) are in great demand due to the rapid development of portable electronics, electric vehicles, and storage systems for solar and wind energies.[1−3] High energy-density rechargeable battery systems, such as lithium-ion batteries, lithium–sulfur batteries, and lithium–air batteries are of great interest in the development of next-generation ESDs.[4−9] However, there are still persistent critical issues to be addressed, especially for the ultrahigh energy battery systems.[3,10,11] The notable volume change of high-capacity electrodes such as silicon and sulfur, calls for unique electrode architectures or fabrication processes in order to stabilize the electrodes.[12] More severely, the electrochemical intermediates of sulfur (polysulfides) that are able to dissolve in liquid electrolytes bring about a shuttle effect, which significantly deteriorates the battery performance.[2,13−17] To resolve these issues, tremendous efforts have been devoted to the fabrication of nanostructured porous electrodes enabling fast ionic/electronic transport and buffering of volume expansion.[18−22] Meanwhile, functional battery interlayers with rational porous structures have been demonstrated to have good abilities to facilitate ionic transport and suppress the diffusion/migration of polysulfides.[23−28] Therefore, the generation of good and functional porous structures plays a vital role in achieving high energy and power densities for ESDs.

A variety of strategies have been implemented in the fabrication of porous structures for different battery applications, but it remains challenging to rationally design and fabricate a porous structure with the desired properties and functionalities.[21,29−32] Recently, Fu et al.[19,23] reported successful fabrications of a three-dimensional porous nanostructured electrode and an interlayer through a protein-directed self-assembly process. As reported, unique porous structures with protein functionalized surfaces can facilitate the ionic transport, increase the binding between electrode and electrolyte, and effectively trap the dissolved polysulfides. The overall electrochemical performance of the battery has been greatly improved by employing these protein-based porous structures. The work has demonstrated that the porous structures generated from the protein-directed self-assembled process pose great potential in the development of high-performance ESDs. In addition, it was also found in experiments that the final porous structures were strongly dependent on the coating materials and the specific solvent environments during fabrication.

Although porous structures from protein-coated nanoparticles have been fabricated and characterized in experiments,[19,23] the picture of a protein-directed self-assembly process under different conditions in a solution state is largely unclear and difficult to measure in experiments. The assembly process is complex and involves complicated interactions among solvent molecules, proteins, and particles. Understanding the molecular scale interactions and the fundamental mechanisms during a self-assembly process in the solution state is significantly important for further design and optimization of porous structures for battery applications. Simulations with certain molecular details play key roles in uncovering the molecular mechanisms involved in such complicated systems. Among different simulation algorithms, the dissipative particle dynamics (DPD) method is an effective mesoscale coarse-grained simulation technique with reasonable accuracy and efficiency. The DPD method has been proven to be a useful approach to investigate the phase behavior of soft matter systems.[33−36] Recently, the DPD method has been successfully employed to study the self-assembly process of micelles, fibers, nanorods, vesicles.[37−41] and polymer systems.[42−46]

In this work, we systematically investigate the self-assembly process of protein-coated carbon black (CB) particles in the solution state through comprehensive DPD simulations. We will explore the effects from different coating materials [polyvinylidene fluoride (PVDF) polymer or gelatin protein], surface coating density, and solvent environments [N-methyl-2-pyrrolidone (NMP), deionized (DI) water or acetic acid (AA)/DI mixture], on the self-assembled structures. The simulation results and conclusions may provide guidelines for further improvement of the fabrication process and eventually lead to better performance of ESDs.

Materials and Methods

Experimental Methods

For completeness, we have redone our previous experiments[23] including fabrication of the CB nanocomposites with different coating materials and characterization of the morphology and porosity of the final self-assembled porous structures. The detailed descriptions of the experiments are provided in the Supporting Information.

DPD Algorithm

The simulations presented in this paper are based on the DPD method.[34] The DPD is a mesoscopic coarse-grained simulation method for soft materials. In DPD simulations, a small group of atoms are coarse-grained into a single bead and the dynamics of each bead is governed by Newton’s equation of motion. The total force acting on bead i is calculated from three types of pairwise forces:Here, C is the conservative force, D is the dissipative force, and R is the random force from bead j within a cut-off distance rC. Each term in eq can be calculated aswhere a is the maximum repulsion between beads i and j, the value can be calculated from the liquid compressibility and solubility (see eq below); is the unit vector and r is the distance between beads i and j; = – is the relative velocity; is the weight function; σ is noise strength and γ is the friction coefficient, and they follow a relation: σ2 = 2γkBT; kBT is the unit of energy, and the standard values of 3.0 and 4.5 are used for σ and γ, respectively; ζ denotes a random number with zero mean and unit variance. All the variables are nondimensionalized in the DPD simulations by the mass of solvent beads m0, the cut-off distance rC, and energy unit kBT.

Coarse-Grained Models

In this paper, we simulate the self-assembly of PVDF- or gelatin-coated CB particles in different solutions. Figure shows the schematic representations of different types of coarse-grained beads for all the components that have been studied in our simulations. Specifically, two water molecules are grouped into one DPD bead and denoted as type (W) (Figure a). As a result, the effective volume of a DPD bead is ∼60 Å3. Accordingly, the interaction cut-off distance rC ≈ 5.65 Å. CB particles are constructed as a sphere with a diameter of 3.0 by setting the DPD bead (CB) on a FCC lattice with a lattice constant of α = 0.34 nm. The CB particles are modeled as rigid body throughout the simulations (Figure b). For solution molecules, each AA molecule is grouped into one bead (A) (Figure c), and each NMP molecule is treated as two beads connected by a spring force with a force constant of 100. The NMP molecule is represented as the DPD bead type (NMP) (Figure d).

Figure 1

Schematic representation of molecular structures and the coarse-grained DPD models for (a) water, (b) CB particles, (c) AA molecules, (d) NMP molecules, (e) PVDF monomer, and (f) chemical structure of gelatin. In the following figures, unless otherwise stated, the color codes are the same as this figure.

The PVDF polymer (−[C2H2F2]−) is modeled as bead-spring chains that consist of 9 consecutive beads (P). One monomer of the PVDF polymer has an effective volume of ∼60 Å3.[47] Therefore, each bead represents one monomer −[C2H2F2]– (Figure e). The gelatin structure is more complicated as shown in Figure f. The chemical structure of gelatin is arranged in the following sequence: -Ala-Gly-Pro-Arg-Gly-Glu-4Hyp-Gly-Pro-.[48] In our model following the mapping procedures in refs,[49,50] glycine with an effective residue volume around 60 Å3 will be treated as a single bead[51] and denoted as the skeletal bead (S). Accordingly, the −NH–CH2–CO– structure connecting adjacent residues will be represented as skeletal beads (S) as well. In addition, proline and hydroxyproline are modeled as skeletal beads (S) bonded to a side bead (T). Arginine is a positively charged side chain and is mapped as three connecting beads, where the positively charged bead (H) is attached to a noncharged side bead (T) and then connected with the skeletal bead (S). The negatively charged glutamic acid is modeled as a skeletal bead (S) connected with a negatively charged bead (H). Therefore, a typical structure of gelatin is represented as 9 connected skeleton beads (S) plus some side-chain beads (T) or (H).

Two consecutive beads (i and j) in PVDF and gelatin are connected by a spring force B = Kb(r – r0), where Kb and r0 are the spring constant and equilibrium bond length, respectively. Here, Kb = 100, r0 = 0.6rC is used for gelatin[49] and Kb = 20, r0 = 0.85rC is used for PVDF.[42] The harmonic angle style potential Uθ = Kθ(θ – θ0)2 is applied on two consecutive bonds for both gelatin and polymer chains, where Kθ, θ, and θ0 are the bending constant, inclination angle, and equilibrium angle. For three consecutive skeleton beads (S) of gelatin or three consecutive PVDF beads (P), the bending constant and equilibrium angle are set as Kθ = 6 and θ0 = 180° to keep the chain structures. For the angle between the S–S bond and S–T bond on gelatin, the parameters chosen as Kθ = 4.5 and θ0 = 90° were taken for the skeleton beads and side-chain beads of gelatin. Gelatin chains and polymer chains were anchored to the surface of a CB nanoparticle by a spring force with Kb = 200, r0 = 0.2. Counterions were added to neutralize the system in the simulation, and represented as bead (e).

Force-Field Parameters

In this work, the conservative interaction parameter between the same type of particles is set to a = 25kBT/rC, and for two beads of different types (e.g., A-type beads and B-type beads), the repulsive parameters can be calculated as[34]where aAB is the repulsive parameters between A-type beads and B-type beads and χAB is the Flory–Huggins interaction parameter. The value of χAB will depend on the solubility parameter of each component, the temperature, and the average volume of the beads:[52]where Vbead is the average volume of a DPD bead; δA and δB are the solubility parameters of bead type A and B, respectively. The repulsive parameters implemented in this work are either calculated from eqs and 6 or adopted from published papers.[49,53] A summary of all repulsion parameters is presented in Table below.

Table 1

Interaction Parameters in Terms of kBT/rC between Different Beads Used in the Simulations, All Symbols Are the Same as Those Used in Figure

	A	CB	H	S	T	W	e	P	NMP
A	25	25	40	33	25	61	61
CB		25	40	33	25	62	62	25	25
H			25	33	33	25	25
S				25	33	30	30
T					25	38	38
W						25	25
E							25
P								25	25
NMP									25

The electrostatic interactions between charged beads are modeled using the “smeared charge” approach with the Slater-type charge density distribution which has an exponential decay. This model can remove the divergence of the Coulombic interactions at r = 0.[54,55]Here, λ is the decay length of the charge. The pairwise electrostatic force E between charged particles i and j at short range with a truncation of 3.0 in real space can be computed as

At long range, the electrostatic interaction of smeared charges reduces to the Coulombic potential and the standard Ewald summation with the integer cut-off 16 in reciprocal space. The value of λ is chosen such that the reciprocal of the decay length β = 1/λ = 0.929. The charge distribution on the protein side chains is determined from the pH values of the solution and will be discussed later.

For mixed solvents such as the AA/DI mixture in this study, the conventional DPD method has significant drawbacks in predicting the physical properties.[56] It is due to the fact that the repulsive parameter a is calculated by solubility parameter δ, where the contribution of hydrogen bonds is implicitly included based on the Hansen approach. Recently, Kacar and de With[57] have solved this problem by adding a Morse potential term VMorse(r) = D0[e–2ϵ(–2 e–ϵ(] to the conservative DPD potential to explicitly model the hydrogen bond interactions between mixture solvent beads (A and W in this work). Here, D0 and ϵ represent the depth and width of the potential, and r0 is the hydrogen bond equilibrium distance. The value of ϵ is set as 2. The cut-off value for the hydrogen bond is taken as the same as that of the nonbonded cut-off of DPD potential. We have performed all-atom molecular dynamics simulations of a AA/DI mixture with a weight ratio of 8:2, the hydrogen bond equilibrium distance is estimated from the radial distribution functions (RDFs) between the center of mass of the hydrogen bond molecules as r0 = 0.45. The value of D0 was estimated from the potential differences of the mixture and the pure components[56] as D0 ≈ 12kBT.

Simulation Set Up

Three self-assembly systems: PVDF-coated CB particles (PVDF/CB) in NMP, gelatin-coated CB particles (gelatin/CB) in DI water, and gelatin-coated CB particles (gelatin/CB) in an AA/water (AA/DI) mixture were prepared through PACKMOL[58] according to the experimental conditions. PVDF/CB in the NMP solvent was simulated with the standard DPD method. Gelatin/CB in DI water and gelatin/CB in the AA/DI mixture solvent were simulated using the modified DPD algorithm with electrostatic force and Morse potential included (Section ). In our simulations, the velocity-Verlet integration algorithm is used to update the positions of each bead, and the integration time step Δt is set to 0.02τ, here is the time unit in our simulations. All the variables in DPD simulations are nondimensionalized by the cut-off distance rC, bead mass m0, and energy kBT. The size of the simulation box is 48rC × 48rC × 48rC, which contains a total of 331 776 beads with the number density of ρ = 3/rC3. All simulations were performed in the NVE ensembles, and the periodic boundary conditions are adopted in all three directions. The DPD units can be converted into SI units based on the mapping scheme: , m0 ≈ 60 g/mol, and τ ≈ 2.8 ps. All simulations in this work were carried out by using the LAMMPS[59] for 500 000 timesteps. Three independent realizations in each case were performed for statistical consistency.

Results and Discussion

Our previous experiments[23] demonstrated that the final porous structures are highly dependent on both the coating materials and solvent environments. (Please see the Supporting Information for related descriptions and discussions.) The molecular interactions and self-assembly behaviors of coated particles in different solvents play critical roles in the determination of the final porous structures after evaporation. Therefore, we perform a series of DPD simulations to systematically explore the fundamental mechanisms involved in the self-assembly of PVDF- or gelatin-coated CB particles. Mimicking the experimental conditions, three cases with different coatings and different solvents have been implemented: (1) PVDF/CB in NMP; (2) gelatin/CB in DI water; and (3) gelatin/CB in the AA/DI mixture. In each case we run three independent realizations for statistical consistency. All realizations yield similar results, therefore only one of the realizations is shown in the following sections for each case.

DPD Studies on PVDF-Coated CB Particles in NMP

We first setup our simulations to study the self-assembly of PVDF-coated CB particles in NMP (the simulation details can be found in the Materials and Methods section). The number of the coating PVDF chains is varied from 1 to 3 to account for the weak interaction between PVDF and CB. Figure shows the equilibrium snapshots of the systems at the end of the simulations. As shown, the PVDF-coated particles are well dispersed in the simulation box regardless of the coating number. In addition, we did not observe any assembled structures throughout the simulations. NMP is well known as a good organic solvent to disperse CB in the experiment, which is in good agreement with the observations from the simulations.

Figure 2

DPD simulations of the PVDF-coated CB particle assembly in NMP. The equilibrium distribution of CB particles with (a) 1 PVDF chain coating, (b) 2 PVDF chains coating, and (c) 3 PVDF chains coating. NMP beads are not shown for clarity.

DPD Studies on Gelatin-Coated CB Particles in DI Water

The experiments (see Figure S1 and associated discussions in Supporting Information) indicated that in DI water the interactions between gelatin and CB are weak and a large amount of gelatin chains cannot be attached on the CB particle surface. Therefore, a small number (1–3) of gelatin chains are attached on the particles in our DPD simulations. As shown in Figure , a typical structure of gelatin contains two charged functional residues (arginine and glutamic acid) on side chains, both of which are denoted as (H) beads in our coarse grained model. The distribution of charges on the residues is sensitive to the pH value of the solvent. In this work the values of the charges are determined by the pKa values of amino acids and the pH values according to ref (60).where the positive sign accounts for basic amino acids (such as arginine on gelatin) and negative sign accounts for acidic amino acids (such as glutamic acid on gelatin).[61] According to this equation, arginine has a positive charge, while glutamic acid has a negative charge. As indicated in Figure a. The overall charge on each chain is neutral. Figure b–d show the equilibrium structures as the number of gelatin coating changes from one to three. As illustrated, clear self-assembled “rod-like” structures are formed in all three cases. The number of gelatin coatings also affects the structure details. As shown in Figure b, small solitary clusters and several “rod-like” segments are formed with only 1 gelatin chain coated onto the CB surface. The length of the “rod-like” segment is short and some of the segments are connected through a small aggregate. When the gelatin coating number is increased to two, as illustrated in Figure c, the length of “rod-like” segments is increased, and all segments are connected by a few ball-like clusters. As the gelatin coating is further increased, as shown in Figure d for three gelatin coating, longer rods with less branches are formed.

Figure 3

DPD simulations of the gelatin-coated CB particle assembly in DI water. (a) Charge distribution on gelatin side chains in DI water. The equilibrium distribution of CB particles with (b) 1 gelatin coating, (c) 2 gelatin coating, and (d) 3 gelatin coating. Water beads are not shown for clarity.

DPD Studies on Gelatin-Coated CB Particles in the AA/DI Mixture

The situation changes dramatically in the AA/DI mixture at pH = 2. First of all, according to eq glutamic acid has zero charge, but arginine remains positively charged. Consequently, the net charge of the gelatin is no longer neutral, but positive charged in acid solution. It is well known that AA is an effective denaturant helping to open up the protein structures and expose protein functional residues.[19,23,62] As a result, we expect stronger gelatin–CB interactions and more gelatin chains on each particle surface. Therefore, in the AA/DI mixture we study the assembly of CB particles with three, four, and six gelatin coatings. Figure shows the self-assembled structures at equilibrium states for each case. As shown, different from the gelatin/CB in the DI water system, “sheet-like” structures are formed in AA/DI solution. As the gelatin coating is increased, the sheet surface becomes larger, thinner, and more curved.

Figure 4

DPD simulations of the gelatin-coated CB particle assembly in the AA/DI mixture. (a) Charge distribution on gelatin side chains in the AA/DI mixture. The equilibrium distribution of CB particles with (b) 3 gelatin coating, (c) 4 gelatin coating, and (d) 6 gelatin coating. Water and AA beads are not shown for clarity.

Simulations demonstrate that the assembly structures of the CB particles strongly depend on the coating materials and solvent environments. To further characterize the different structures, we calculate the RDFs of CB particles for all the cases. The RDF is defined as , where d(r) is the ensemble averaged number of particles around a referenced particle within the volume of a shell from r → r + Δr; ρ = N/V is the number of particles per unit volume; V is the volume of the system and N is the total number of particles in the system. As shown in Figure a, the RDF for PVDF/CB in NMP is flat, indicating that the particles are well dispersed in solution and there is no clear assembled structures. However, for all gelatin/CB systems as clearly shown in Figure , the distributions peak at rpeak ≈ 3.8 indicates the accumulation of particles. For gelatin/CB in DI water, the distribution peak gpeak ≈ 90 is observed when there is only one gelation coating. When the coating number is increased, the peak value is slightly decreased to ∼75. For gelatin/CB in the AA/DI mixture, the rpeak values shift slightly to the right, and the gpeak values are much smaller compared with the DI water cases. The larger value of rpeak indicates that the particles are more loosely packed, and the assembled structures outspread more at a smaller value of gpeak. Figure b,c illustrate the structural details for the two gelatin/CB systems. As is shown, in DI water, the “rod-like” strings are always triple-strand and the CB particles are closely stacked together. In the AA/DI mixture, the outstretched “sheet-like” surfaces are mostly formed with single layer loosely packed CB particles. These structural features are consistent with the corresponding RDF distributions.

Figure 5

Details of the self-assembled structures. (a) RDF of CB particles for each simulated case. (b) Compact triple-strand “rod-like” structures formed for gelatin/CB in the case of DI water. (c) Loose curved “sheet-like” structures formed for gelatin/CB in the case of the AA/DI mixture. All solvent and gelatin beads are not shown.

Solvent and Charge Effects on Final Self-Assembled Structures

Both the solution and the redistribution of the charges on gelatin side chains may contribute to the structural differences for gelatin/CB systems. To differentiate the effects, we have studied the assembly of pure CB particles in the DI water and AA/DI mixture. The final self-assembled structures of bare CB particles in both solutions are illustrated in Figure a. As shown in the inset, large irregular agglomerations were formed in both solvent conditions. Moreover, both RDF profiles showed similar peaks at the similar locations (marked by green dashed squares and denoted as the 1st peak and 2nd peak). The 1st peak in the AA/DI mixture is slightly shifted to the right because the AA molecules have stronger attraction with the CB particle (see Table ). Both RDF profiles and morphologies in Figure a indicate that the self-assembled structures from these two solvent environments are statistically similar, suggesting that the structural difference is not because of the interactions from different solvent molecules. To clarify the effect from the charges, we arbitrarily set the charges on the side chain of gelatin to zero and re-ran the case for six gelatin-coated CB particles in the AA/DI mixture. As shown in Figure b, all the particles accumulate into a cylindrical structure. Based on Figure , the dominant factor leading to the “sheet-like” structures is the charge redistribution on the side chains under different solvent environments.

Figure 6

Solvent and charge effects on self-assembled structures in solution. (a) RDF of pure CB particles’ self-assembled structures in water (I) and in the AA/DI mixture (II). Inset shows the structures at the equilibrium state. (b) Self-assembled structure when charges on gelatin are set to zero in AA solution. All solvent beads are not shown.

Our DPD simulations show that distinctive structures are formed with different coatings in different solutions. The structures at the solution state will directly impact the final porous structures at the solid state after evaporation. Our simulation results show that the PVDF-coated CB particles are well dispersed in NMP as illustrated in Figure a, therefore during evaporation the accumulation of CB particles is random and uncontrollable (see Figure d), in the end large cracks may appear at different locations (see Figure S1a in Supporting Information). For gelatin-coated CB particles in DI water, interconnected and closely packed “rod-like” structures are formed (see Figure b). As illustrated in Figure e, during the drying process, the rods can easily break into sticks and randomly stack together. More importantly, due to the weak interactions between gelatin and CB particles, excessive amounts of free gelatins are in the solution and become solid blocking the porous structures after evaporation. As a result, the porosity is extremely low in this case (∼12% see Figure S1b in Supporting Information). For gelatin-coated CB particles in the AA/DI mixture, the particles assemble into curved “sheet-like” structures. These structures are stable and flexible because of the strong gelatin–gelatin interactions. It is most likely that those curved sheets will form homogeneous porous structures (porosity of ∼45% see Figure S1c in the Supporting Information) with large channels as demonstrated in Figure c after evaporation. Additionally, as indicated in Figure c,d, the large sheet surfaces are covered with numerous charges (shown in red color), producing repelling forces during evaporation, which will further stabilize the self-assembled structure. As shown in Figure f, a well-maintained structure will be obtained even at a dry state. The homogeneous porous structures with functionalized surfaces work especially well as an interlayer between the cathode and the separator in Li–S batteries. As demonstrated by the simulations and experiments in ref (23), the resulting interlayer provides more pathways for Li+ transport, and the functionalized surfaces effectively trap the polysulfides suppressing the “shuttle” effects.

Figure 7

Schematic illustration of the structures in the solution state (a–c) and in dry state (d–f) after evaporation.

Conclusions

Our experiments show that the final porous structures of the CB particles strongly depend on the coating materials and solvent environment. Coating the CB particles with gelatin proteins and dispersing them in an AA/DI mixture result in a homogeneous porous structure with high porosity. The accumulation of protein-coated particles in solution is a complex and protein-directed self-assembly process, which involves complicated interactions among solvent molecules, proteins, and CB particles. Here, we have implemented the comprehensive DPD simulations to investigate the self-assembly processes of CB particles with different coatings (PVDF or gelatin) at different solutions (NMP, DI water, or AA/DI mixture). Our DPD simulations show that PVDF-coated CB particles are well dispersed in the organic NMP regardless of the coating density. Gelatin-coated CB particles result in triple-strand and compact “rod-like” structures in DI water, while gelatin-coated CB particles in the AA/DI mixture lead to thin and curved “sheet-like” structures. The redistribution of the charges on gelatin side chains under different solvent environments plays key roles in the determination of the assembled structures. The structures in solution directly impact the final porous structures after removal of the solvent. The simulation results are able to reasonably explain the experimental observations. Both simulation and experiments demonstrate that homogeneous high porosity structures with functionalized surfaces can be obtained through a protein-directed self-assembly process under appropriate conditions, and DPD simulations represents a powerful tool in guiding future experimental design and optimization. Finally, in DPD coarse graining several atoms into one bead and grouping intermolecular forces into several short-range forces make the simulations efficient and especially suitable for study of mesoscale self-assembling processes. However, due to the coarse graining the detailed molecular interactions dictating the process cannot be explicitly explored, and this information can be further investigated in future work using molecular simulations.