Homogenization of mechanical properties for 2-D PBX-9501 model
Molecular dynamics simulations yield mechanical properties of the Estane binder and HMX crystals but cannot be used to extract mechanical properties of the PBX-9501 on the scale of micrometers. The MPM is a promising tool for obtaining properties of PBX-9501 from the properties of its pure components, e.g. performing property homogenization.
A number of stress relaxation, creep, and dynamic-mechanical experiments have been performed on the pure materials and composites representative of PBX-9501. The composites have been generated by Monte-Carlo packing of 2-D cylinders having the distribution representative of the range of the distribution of HMX particles in PBX-9501. The resulting composite is shown in Fig. .The 160000 material points have been used to represent polymer and HMX in MPM. Typical MPM runs were 107 timesteps. The stress relaxation experiments are the most suitable for the property homogenization task. Results for two stress relaxation experiments are shown in Fig. . The shear strain was created by applying shear on one side of the composite as show in Figure 1, while keeping the one side of the composite fixed. Stress on the material points have been found to exhibit large oscillations. Our preliminary computational experiments indicated that the shear modulus for the composite not only increases in magnitude and shows significantly longer relaxation times, but the shape of the relaxation becomes significantly different. For example, if pure polymer has a single exponential decay of G(t)=G0exp(-t/t), the shape of the composite has a stretched exponential decay G(t)=G0compexp(-(t/tcomp)b), with b「0.25, therefore, a series of exponents will be required to describe G(t) in the form used in the 3-D ASCI MPM code.

![]()
Force field development of Plasticizers in binder materials
Plasticisers composed of bis(2,2-dinitropropyl)acetal (BDNPA) and bis(2,2- dinitropropyl) formal (BDNPF) have found widespread application in energetic formulations. For example, US Army applications such as the M900 tank program, LOVA gun propellant, PBXN-106 used to fill the US Navy 5.54 gun contains BDNPA/F and the HMX-based explosives modeled by the C-SAFE center use BDNPA/F as an energetic plasticizer. We developed a classical force field for the BDNPA.
Because quantum chemistry calculations on the MP2/aug-cc-pvDz//B3LYP/aug-cc-pvDz level, required for predicting conformational energetics with chemical accuracy, are computationally expensive taking 15-20 days for each point to complete, we followed the previously tested methodology of developing valence force field on shorter model compounds representative of the parts of BDNPA and transferring the obtained valence to the BDNPA compound. Relative conformational energetics of model compounds shown in Figure X has been studies at the MP2/aug-cc-pvDz//B3LYP/aug-cc-pvDz level. Partial charges for the model compounds were obtain by the restrictive fit to the electrostatic potential around each model compound for the lowest energy geometries, equilibrium bond length and bending distances were fitted to the lowest energy conformers. Torsional potential parameters were fit to obtain the best representation of the most important torsional isomerization paths as shown in Figure Y.
The ability of the force field to predict conformational energetics of the selected conformers of BDNPA calculated at the MP2/aug-cc-pvDz//B3LYP/aug-cc-pvDz is in progress.
![]()
![]()

![]()


We conducted extensive quantum chemistry calculations on estane model compounds. The quantum chemistry data were used to parameterize a fully-atomistic force field for classical MD simulations of estane. Currently, MD simulations of melts of relatively long estane chains as well as mixtures of shorter estane segments are under way. The focus of this study is to obtain a fundamental understanding of local and mesoscale structure, conformational properties, self-assembly of rigid estane segments as well as influence of these properties on dynamical and viscoelastic behavior of estane melts at various thermodynamic conditions. In addition to atomistic simulations we have initiated coarse-grained simulations of estane. In our coarse-grained representation each bead in a polymer chain represents several atoms allowing us to significantly expand the length and time scales accessible in MD simulation as well as investigate a broader range of parameter space. Using umbrella-sampling techniques that allow us to calculate the potential of mean force between various segments of estane chain in the melt it will be possible to map to a coarse-grained model that reproduces structural, thermodynamical and conformational features of estane melt predicted from MD simulations using the fully atomistic potential.
In order to carry out atomistically detail modeling of PBX-9501 we initiated high level quantum chemistry calculations for the model compounds of the plasticizer molecules (CH3-C(NO2)2-CH2-O-CH(R)-O-CH2-C(NO2)2-CH3 where R= CH3 or H). In addition to understanding the bulk properties of polymeric binder (estane+plasticizer) and crystalline organic explosive (HMX) it is very important to understand the influence of the interface between these components on the structural, conformational and dynamical properties of the constituents as well as composite properties overall. We also initiated MD simulations of an HMX/estane interface using atomistically detail potential functions. In Fig. 3 a snapshot of an a-HMX slab surrounded by an estane melt is illustrated. While in the future the interfaces between most important planes of b and d-HMX and estane will be investigated, the current study of a-HMX/estane systems allows us to investigate two chemically very different HMX surfaces: an O and N exposed surface and a CH2 exposed surface, therefore providing insight into the importance of crystal orientation and the nature of the surface-exposed crystal plane on the polymer behaviour at the interface.
In order to gain fundamental understanding of the influence of filler particles and their interaction with the polymer matrix on viscoelastic and dynamic properties of the polymer chains as well as the composite overall, we have conducted extensive simulations of particle filled polymer composites using MD simulations with a coarse-grained model. In this study it was found that the viscosity of the composite scales proportionally to the specific surface area of the filler particles. It was also found that repulsive interactions (relatively to polymer-polymer interaction) between filler and the polymer results in a speed up of polymer dynamics at the interface. Stronger attractive interaction between particle and polymer results in significant slowing down of polymer dynamics as well as significant heterogeneity of polymer dynamics on all length scales. The latter complicates application of conventional theoretical models (e.g. Rouse analysis) in describing polymer matrix dynamical and viscoelastic properties. While investigated particle filled polymer composites do not represent PBX-9501 explicitly the fundamental knowledge of the influence of filler-polymer interactions on composite and polymer matrix properties will be very valuable for future simulations of PBX-9501.
In addition to simulations of particle filled polymers, simulations of coarse-grained systems with self assembly has been performed to understand the underlying physics of self-assembly and correlations between structural and dynamical properties in self-assembling polymers. In particularly, the attention was focus on simulations of polymer chains in which 10% of the beads have a strong attraction with each other ("sticky" beads) and therefore tend to aggregate in micelles or clusters. Such aggregation results in formation of a continuous polymer network where polymer chains are "cross-linked" by self-assembled clusters which is similar to characteristics of estane melts where MDI rigid units tend to aggregate. We have focused on uncovering the correlations between single-molecule events (hoping of sticky beads from cluster to cluster, dynamics of sticky beads inside the cluster, etc.) with macroscopic viscoelastic properties of the system (complex shear modulus, viscosity, etc.). Also some effort was dedicated to understanding the influence of different coarse grained models (lattice Monte-Carlo, implicit solvent MD simulation and explicit solvent MD simulation) on the predicted properties of self-assembled systems.
Development of Parallel Tempering Methodology for Simulating Self-Assembling Polymers (Estane)

![]()
In order to overcome slow equilibration in polymeric systems, especially in systems with self-assembly like estane where MDI segments have a strong tendency to form nanoscale clusters, Smith's group has implemented and continue to develop one of the latest simulation methodologies - parallel tempering (PT). Parallel tempering (PT) is an algorithm that has been recently developed to dramatically speed up the equilibration and sampling of configurational phase space in MD and MC simulations. The algorithm requires parallel simulation of N identical systems at series of different thermodynamic conditions, e.g. temperatures (T1, T2, T3,...TN) with occasional switching the coordinates between the systems. Systems at higher temperatures are able to explore phase space much faster than systems at lower temperatures. An appropriate Monte Carlo protocol for switching the coordinates between systems based upon statistical mechanics allows each system to be simulated at all temperatures. After some simulation time each system visit higher temperatures and return to lower temperatures with a significantly different configuration from the original one. The PT algorithm is schematically illustrated in Fig. 4.

![]()
We has successfully demonstrated that PT simulations can significantly speed up equilibration and sampling of conformational properties in polymer melts as well as in self-assembling polymer networks. In Fig. 5 this is illustrated by comparing polymer end-to-end autocorrelation function (ACF) obtained from PT and conventional MD simulations. While tempering on temperature was found to be potentially effective, over the last year the drawbacks of using global tempering parameter ( changes in which influence all interactions in the system simultaneously) such as temperature or density were realized. In particularly, even small changes in the tempering parameters (e.g. DT~2o) significantly shifts the distribution of potential energy of large systems and therefore requires one to have a large number of replicas (on the order of 100) to provide a chain of sufficiently overlapped distributions between the lowest (temperature of interest) and the highest (temperature at which relaxation of the system is relatively fast) temperatures. Having a hundred replicas, results in a long waiting time in order for systems to travel from the lowest to the highest temperature - a necessary condition for relaxation of each system. While such PT simulations would still be computationally efficient in comparison with conventional MD wall-clock time limitations make PT runs essentially impractical.
Figure
6
In order to overcome this bottleneck Smith's group in collaboration with W. Paul (Mainz, Germany) is working on enhancing PT methodologies by developing tempering parameters that create localized perturbations and target relaxation limiting process. In particular, they implemented PT methodology with tempering on an external field that promotes some of the polymer chain atoms to be phantom (not interacting with other intermolecular atoms) or equivalently be displaced to a 4th dimension as illustrated in Fig. 6. This procedure creates localized regions of reduced density that allow conformational and structural (local) relaxation that otherwise would be arrested by tight intermolecular packing. Because this is a localized perturbation the overall change in potential energy of the system is not significant and therefore relatively few replicas (up to 20) are required to cover a range of the tempering parameter (strength of external field) from a value where no atoms are promoted to the 4th dimension (system of interest) to values where a significant fraction of atoms (up to 20%) are in the 4th dimension and dynamics of the system is enhanced by several orders of magnitude. Smith's group has implemented and validated this technique and, currently, is optimizing the performance and conducting production runs for a polybutadiene melt in the vicinity of the glass transition.