Isobaric Simulations Research Articles

Extraction of Thermal Properties of Organic Ablative Materials Using Molecular Dynamics Simulations

In this study, molecular dynamics method is used for extracting thermal properties of materials relevant to the ablative type of thermal protection system. The methodology and results are discussed for extracting thermal conductivity and specific heat capacity of ablative materials for their virgin and char states in a manner that is defined in the experimental data sets. To select a suitable nonreactive force field, Green–Kubo simulations are performed over phenol, toluene, and benzene materials using Amber, Dreiding, and all-atom optimized potentials for liquid simulations (OPLSAA) force fields. Next, crosslinked and non-crosslinked polymers are created, and pyrolysis simulations are performed over them using ReaxFF to obtain virgin and char states of the polymers. Thermal conductivity is extracted using the Green–Kubo method, whereas specific heat capacity is extracted by using the enthalpy gradient obtained by performing isobaric simulations over virgin and char materials. Later, by introducing the concept of an artificial composite material, the extracted thermal properties are used in macroscale one-dimensional heat conduction simulations with pyrolysis using an in-house finite element method (FEM)-based thermal response solver.

Exploring Librational Pathways with on-the-Fly Machine-Learning Force Fields: Methylammonium Molecules in MAPbX3 (X = I, Br, Cl) Perovskites.

Two seemingly similar crystal structures of the low-temperature (∼100 K) MAPbX3 (X = I, Br, Cl) perovskites, but with different relative methylammonium (MA) ordering, have appeared as representatives of this orthorhombic phase. Distinguishing them by X-ray diffraction experiments is difficult, and conventional first-principles-based molecular dynamics approaches are often too computationally intensive to be feasible. Therefore, to determine the thermodynamically stable structure, we use a recently introduced on-the-fly machine-learning force field method, which reduces the computation time from years to days. The molecules exhibit a large degree of anharmonic motion depending on temperature: that is, rattling, twisting, and tumbling. We observe the crystal’s “librational pathways” while slowly heating it in isothermal–isobaric simulations. Marked differences in the thermal evolution of structural parameters allow us to determine the real structure of the system via a comparison with experimentally determined crystal structures.

Non-equilibrium hybrid insertion/extraction through the 4th dimension in grand-canonical simulation

The process of inserting/deleting a particle during grand-canonical Monte-Carlo (MC) simulations is investigated using a novel, original technique: the trial event is made of a short nonequilibrium molecular dynamics (MD) trajectory during which a coordinate w along a 4th dimension is added to the particle in the course of insertion/deletion and is forced to decrease from large values down to zero (for insertion) or increased from 0 up to large values (for extraction) at imposed vw velocity. The probability of acceptation of the whole MC move is controlled by the chemical potential and the external work applied during the trajectory. Contrary to the standard procedures which create/delete suddenly a particle, the proposed technique gives time to the fluid environment to relax during the gradual insertion/extraction before the acceptation decision. The reward for this expensive trial move is a gain of many orders of magnitude in the success rate. The power and wide domain of interest of this hybrid "H4D" algorithm which marries stochastic MC and nonequilibrium deterministic MD flavors are briefly illustrated with hard sphere, water, and electrolyte systems. The same approach can be easily adapted in order to measure the chemical potential of a solute particle immersed in a fluid during canonical or isobaric simulations. It then becomes an efficient application of the Jarzynski theorem for the determination of solvation free energy.

Elastic Properties of Swelling Clay Particles at Finite Temperature upon Hydration

The knowledge of mechanical properties of swelling clays on the scale of the clay layer is of crucial importance to build predictive hydromechanical macroscopic constitutive laws in a bottom-up approach. In this work, we computed the elastic properties of particles of a swelling clay (i.e., montmorillonite) by molecular dynamics simulations. Because of the softness of the material in the direction orthogonal to the clay layer with increasing water content, the computation of the whole stiffness tensor was not trivial: we needed to implement the elastic bath method,1 which allows us to attenuate thermal strains and compute the stiffness via strain fluctuations in isothermal–isobaric simulations. We investigated the effect of water content, temperature, and the interlayer cation on the elastic properties of swelling clays. In particular, we showed that the out-of-plane stiffness coefficients computed at 300 K differed significantly from the same coefficients computed at 0 K. The out-of-plane coefficients were very sensitive to both temperature and water content. In contrast, the dependence of the in-plane coefficients to temperature was slight. Moreover, the decrease of the in-plane coefficients with water content could be entirely explained by the change of geometry of the system due to the swelling of the interlayer space. The interlayer cation impacted the elastic properties of montmorillonite only in the driest states.

Global phase behaviour of polyphilic tapered dendrons

We present a coarse-grained molecular dynamics study into the fluid phase behaviour of a model dendron like polyphilic molecule. Our model consists of 11 spherical beads, held in a semi flexible planar tapered arrangement, with a single apex bead and four interior beads defined to have self attraction, whilst the remaining periphery beads and all cross interactions are made softly repulsive, in essence bestowing an incompatibility between the three bead types. A parametric study was performed covering around 1000 state points in a domain of temperature and pressure spanning from regions where the fully isotropic fluid is observed, to the onset of crystallisation. Three micro-segregated thermotropic liquid crystalline phases are detected; a hexagonal columnar (HC) mesogenic phase, a supramolecular sphere body-centred cubic (SSbcc) mesogenic phase and a supramolecular sphere fluid (SSf) mesogenic phase. Each phase is characterized by order parameters based on dendron clustering distributions and spatial orientations, facilitating the phase diagram mapping. The mesogenic nature of the HC, SSbcc and SSf structured phases was confirmed visually, presenting dendron diffusion along and around the columns in the former, and intra and inter supramolecular sphere dendron diffusion in the latter two. HC and SSbcc phases were obtained by self-assembly by quenching from isotropic states, proving the stability and reproducibility of the liquid crystalline phases. Tracking of appropriate metrics gives insight into their formation mechanisms. Temperature driven first order phase transitions are seen from the HC phase to the SSf phase at higher pressures, and to the SSbcc and then SSf phase at lower pressures. For this latter transition, a Gibbs Duhem integration is used to trace the coexistence curve from a point determined by thermodynamic integration. The other transitions, including a continuous transition from the SSf to fully isotropic phase, are estimated from isothermal and isobaric simulation runs. It is of note that both the HC and SSbcc mesogenic phases resemble those reported experimentally however, the intermediate liquid-phase is an undiscovered prediction of our model. The simulations presented give physical insight on the relation between the apparent distinct phases observed. Our results shed light on the underlying relationships between these phases of self-organized supramolecular soft matter and showcase how the different experimental instances are local manifestations of an underlying generic phase diagram, which we present here.

Nosé–Hoover nonequilibrium dynamics and statistical mechanics

At equilibrium Nosé's 1984 revolutionary thermostat idea linked Newton's mechanics with Gibbs' statistical mechanics. His work expanded the scope of isothermal and isobaric simulations. Nosé–Hoover dynamics has subsequently facilitated the simulation and detailed understanding of nonequilibrium problems. The fractal phase-space distributions, and their close link to the Lyapunov spectrum, provide a novel explanation of irreversibility and a rich field for exploration.

Simulating Fluid–Solid Equilibrium with the Gibbs Ensemble

The Gibbs ensemble is employed to simulate fluid–solid equilibrium for a shifted-force Lennard-Jones system. This is achieved by generating an accurate canonical Helmholtz free-energy model of the (defect-free) solid phase. This free-energy model is easily generated, with accuracy limited only by finite-size effects, by a single isothermal–isobaric simulation at a pressure not too far from coexistence for which the chemical potential is known. We choose to illustrate this method at the known triple-point because the chemical potential is easily calculated from the coexisting gas. Alternatively, our methods can be used to locate fluid–solid coexistence and the triple-point of pure systems if the chemical potential of the solid phase can be efficiently calculated at a pressure not too far from the actual coexistence pressure. Efficient calculation of the chemical potential of solids would also enable the Gibbs ensemble simulation of bulk solid–solid equilibrium and the grand-canonical ensemble simulation of bulk solids.

Hydrogen fluoride phase behavior and molecular structure: Ab initio derived potential models

Several variations of ab initio based molecular models for hydrogen fluoride (HF) are examined by Monte Carlo molecular simulation to determine their bulk-phase properties. The models are taken from the literature, and represent fits of functional forms to the potential energy surface of the HF dimer as given by ab initio computational chemistry calculations. For one of these models, we examine three variations for bulk-phase modeling. In particular, we consider first the effect of including versus neglecting an Ewald sum for the long-range dipole–dipole interactions; second, we examine a modification of the form for the short range repulsive region of the potential; and third, we add three-body contributions to the energy via an available 12-dimensional potential for the trimer, again representing a fit to ab initio energy calculations. The simulations examine the density (via isothermal–isobaric simulation) and radial distribution function (via canonical–ensemble simulations) each at two state points where corresponding experimental data are available. We also examine vapor–liquid coexistence properties, considering the saturation densities, heat of vaporization, and vapor pressure from 225 K to states approaching (but not closely) each model’s critical point. Inclusion of the three-body energy is the only variation that has any beneficial effect on the radial distribution function as compared to experiment, and this variation also gives good results for the vapor pressure, and significantly raises the critical point toward the experimental value. However this model also grossly overestimates the liquid-phase coexistence density. In almost all regards none of the models or variations can be considered to give a satisfactory representation of the bulk-phase behavior. Improvements to the models require more careful attention to the balance between repulsive and attractive pair interactions at short range.

A configurational-bias approach for the simulation of inner sections of linear and cyclic molecules

A novel Monte Carlo technique, the rebridging configurational-bias (RCB) method, has been developed to rearrange inner sections of chain molecules having strong intramolecular interactions along the backbone. The ability of sampling inner chain segments is important for the simulation of systems with low concentration of chain ends, such as polymers and molecules with cyclic structures. In the RCB method, inner segments are removed and then regrown site-by-site in a configurational-bias fashion. However, a short preliminary simulation is carried out to calculate weighting functions based on histograms of the separation distance between pairs of sites along the chain; these functions are used to bias the trial positions of growing inner sites so as to promote efficient chain closure. A look-ahead search scheme is employed for the sampling of the last two sites to increase the overall acceptance rate of the reconstruction process. The validity and performance of the new algorithm were tested by studying linear alkane systems of various chain lengths. Fast conformational equilibration was observed, from the level of local bond orientations to the level of the chain end-to-end vector orientations. Cyclic alkanes containing from 8 to 100 carbon atoms were also studied by using the RCB method. Our results for the conformational properties of cyclooctane are generally consistent with previous molecular dynamics (MD) results and with experimental data. The vapor–liquid coexistence curve of cyclooctane was also mapped out by using the RCB method in conjunction with a histogram reweighting technique for the analysis of isothermal–isobaric simulations.

Extended continuum configurational bias Monte Carlo methods for simulation of flexible molecules

The continuum configurational bias (CCB) Monte Carlo method has been extended to perform elementary moves that involve the rearrangement of inner segments of flexible chains. When regrowing inner sites, the continuity with the rest of the chain is ensured by disregarding those configurations that would imply an unrealistic elongation of the bonds once the chain is reconstructed. The formalism presented here also allows the simulation of branched chains and crosslinked-network structures. The Monte Carlo elementary moves proposed in this work are used in conjunction with an alternative method of preferential sampling in which the segments to be rearranged are chosen from a preselected region of space. The performance and capabilities of the new moves are compared to those of standard CCB and crank-shaft algorithms for simulation of melts and solutions of hard-sphere chains at high densities. Our results indicate that the methods presented here provide a fast relaxation of the bond orientation and the end-to-end orientation autocorrelation functions. Our isobaric simulations for homopolymer chains of up to 51 sites and for concentrated solutions of chain molecules in the monomer are consistent with previously reported data obtained by approximate molecular dynamics methods and by conventional Monte Carlo methods. However, small disagreements with existing data are identified at high densities. These PV results are also compared to the predictions of two recent equations of state. This comparison shows the presence of some small but systematic deviations.

Orientational glassy phases of C60 and neopentane. A Monte Carlo study

Based on an isothermal, isobaric simulation the structure and properties of the plastic crystalline phases of C60 and neopentane have been examined. Instantaneous cooling of the plastic crystalline phases of both C60 and neopentane leads to orientational glassy phases. These are accompanied by significant slowing down of reorientational motion. Constant pressure quench experiments on C60 yield a glass transition temperature of around 80 K.

A dynamic Monte Carlo method suitable for molecular simulations

A dynamic Monte Carlo (MC) method for simulations has been presented which is similar to the well-known Brownian dynamics method, but is suitable for many systems and can readily be adopted to isobaric simulations. This dynamic MC has been applied to several systems and shows similar accuracy and improved efficiency compared with routine methods. It samples configuration space approximately according to the equilibrium probability density and is applicable to NPT simulations of rather complex systems, also with rigid constraints.

Nebula boundary conditions in Monte Carlo simulations

A new alternative to common periodic boundary conditions in Monte Carlo simulations is proposed. A very large system is simulated with elementary one-particle moves performed preferentially in the central region of a simulation cell and local quantities are measured in the same way. Using the Lennard-Jones liquid as a test system it is found that the new method gives more accurate values of thermodynamic quantities than a comparable standard isochoric simulation in periodic boundary conditions and that it is not worse than an isobaric simulation.

Microscopic descriptions of helium at high pressure

Three methods: the correlated cell model, isobaric simulations and self consistent phonons are applied to the calculation of the high density phase diagram of 4He and compared. This leads to a good description of the high density solid insulators that can help us to understand recent high pressure measurements.