Constant Pressure Molecular Dynamics Simulations Research Articles

Prediction of quaternary hydrides based on densest ternary sphere packings

We exhaustively search quaternary metal hydrides based on the (13-2-1) and (13-3-1) structures that are two of the putative densest ternary sphere packings in cubic systems [R. Koshoji et al., Phys. Rev. E 104, 024101 (2021)]. The 73 304 candidate hydrides are generated by substituting the small spheres with hydrogen atoms, and the medium, large, and fourth spheres with metallic atoms. Specifically, the substitution of the small spheres with hydrogen atoms gives the unconventional hydrogen sublattices. We screen unstable hydrides in the candidates through geometrical optimizations, constant pressure molecular dynamics simulations, and phonon calculations under hydrostatic pressure of 10 GPa, and identify 23 hydrides with static and dynamic stability, including ${\mathrm{H}}_{12}{\mathrm{ScY}}_{2}\mathrm{La}$ and ${\mathrm{H}}_{12}\mathrm{Ti}{\mathrm{Ni}}_{3}\mathrm{Ba}$. The superconducting transition temperatures ${T}_{c}$, calculated by density functional theory for superconductors, are found to be 5.7 and 6.7 K for the selected two hydrides ${\mathrm{H}}_{12}{\mathrm{ScY}}_{2}\mathrm{Ca}$ and ${\mathrm{H}}_{12}{\mathrm{ScY}}_{2}\mathrm{Sr}$, respectively. We expect that the 23 candidates of hydrides screened by the exhaustive search for 73 304 hydrides provide a guideline to narrow down the search space for trials in the experimental synthesis of quaternary metal hydrides.

Effect of ion size disparity on the thermal hysteresis of ionic liquids

By means of computer simulations, the effect of varying anion to cation size ratio on the liquid–solid transition and on the thermal hysteresis behavior of a model ionic liquid is investigated. Constant pressure molecular dynamics simulations are carried out by considering different values of anion to cation size ratios at the same number density and 1 atm pressure. In this study, we find that changing ion size ratio affects the solidification or melting temperatures, latent heat, of the model ionic liquid. Within the considered range of anion to cation size ratios, solidification/melting temperature as well as area of thermal hysteresis loop decreases with increasing size ratio and thus illustrates the possibility to fine tune thermal hysteresis loop, which is relevant for phase change materials based thermal energy storage applications.

Assessing the mechanical properties of molecular materials from atomic simulation

<abstract> <p>We review approaches to deriving mechanical properties from atomic simulations with a special emphasis on temperature-dependent characterization of polymer materials. The complex molecular network of such materials implies only partial, rather local ordering stemming from the entanglement of molecular moieties or covalent bonding of network nodes, whereas the polymer strands between the nodes may undergo nm-scale reorganization during thermal fluctuations. This not only leads to a strong temperature-dependence of the elastic moduli, but also gives rise to visco-elastic behavior that complicates characterization from molecular dynamics simulations. Indeed, tensile-testing approaches need rigorous evaluation of strain-rate dependences, provoking significant computational demands. Likewise, the use of fluctuations observed from unbiased constant-temperature, constant-pressure molecular dynamics simulation is not straight-forward. However, we suggest pre-processing from Fourier-filtering prior to taking Boltzmann-statistics to discriminate elastic-type vibrations of the simulation models for suitable application of linear-response theory.</p> </abstract>

Structural and Dynamical Properties of Alkaline Earth Metal Halides in Supercritical Water: Effect of Ion Size and Concentration.

Constant temperature-constant pressure molecular dynamics simulations have been performed for aqueous alkaline earth metal chloride [M2+-Cl- (M = Mg, Ca, Sr, and Ba)] solutions over a wide range of concentrations (0.27-5.55 m) in supercritical (SC) and ambient conditions to investigate their structural and dynamical properties. A strong influence of the salt concentration is observed on the ion-ion pair correlation functions in both ambient and SC conditions. In SC conditions, significant clustering is observed in the 0.27 m solution, whereas the reverse situation is observed at room temperature and this is also supported by the residence times of the clusters. The concentration and ion size (cation size) seem to have opposite effects on the average number of hydrogen bonds. The simulation results show that the self-diffusion coefficients of water, cations, and the chloride ion increase with increasing temperature, whereas they decrease with increasing salt concentration. The cluster size distribution shows a strong density dependence in both ambient and SC conditions. In SC conditions, cluster sizes display a near-Gaussian distribution, whereas the distribution decays monotonically in ambient conditions.

Elucidating water dynamics in MgCl2 hydrates from molecular dynamics simulation

The water mobility in single crystals of MgCl2 hydrates has been investigated with molecular dynamics. Standard force fields have been benchmarked for molecular dynamics simulations of MgCl2 hydrates. To provide a reliable molecular mechanics model, force fields are selected on their ability to reproduce the structure of MgCl2·6H2O at 300 K. The selected force fields are then tested on their ability to also reproduce the structures of the different hydrates (n = 12, 8, 6, 4, 2, 1) and available thermodynamic data. For the currently best force-field available, constant-temperature, constant-pressure molecular dynamics simulations are preformed to elucidate the mechanisms of hydrate water mobility in perfect single crystals of the tetra- and hexahydrate. Long range water diffusion was not observed; each water molecule remained in the coordination sphere of its original Mg2+ host. However, collective ring-like motions of four water molecules at once within the coordination shell of a Magnesium ion were observed.

Equilibrium structures of anisometric, quadrupolar particles confined to a monolayer.

We investigate the structural properties of a two-dimensional system of ellipsoidal particles carrying a linear quadrupole moment in their center. These particles represent a simple model for a variety of uncharged, non-polar conjugated organic molecules. Using optimization tools based on ideas of evolutionary algorithms, we first examine the ground state structures as we vary the aspect ratio of the particles and the pressure. Interestingly, we find, besides the intuitively expected T-like configurations, a variety of complex structures, characterized with up to three different particle orientations. In an effort to explore the impact of thermal fluctuations, we perform constant-pressure molecular dynamics simulations within a range of rather low temperatures. We observe that ground state structures formed by particles with a large aspect ratio are in particular suited to withstand fluctuations up to rather high temperatures. Our comprehensive investigations allow for a deeper understanding of molecular or colloidal monolayer arrangements under the influence of a typical electrostatic interaction on a coarse-grained level.

On the composition dependence of thermodynamic, dynamic and dielectric properties of water-methanol model mixtures. Molecular dynamics simulation results

We have investigated thermodynamic and dynamic properties as well as the dielectric constant of water-metha\-nol model mixtures in the entire range of composition by using constant pressure molecular dynamics simulations at ambient conditions. The SPC/E and TIP4P/Ew water models are used in combination with the OPLS united atom modelling for methanol. Changes of the average number of hydrogen bonds between particles of different species and of the fractions of differently bonded molecules are put in correspondence with the behavior of excess mixing volume and enthalpy, of self-diffusion coefficients and rotational relaxation times. From the detailed analyses of the results obtained in this work, we conclude that an improvement of the description of an ample set of properties of water-methanol mixtures can possibly be reached, if a more sophisticated, carefully parameterized, e.g., all atom, model for methanol is used. Moreover, exploration of parametrization of the methanol force field, with simultaneous application of different combination rules for methanol-water cross interactions, is required.

Composition dependence of thermodynamic, dynamic and dielectric properties of water–methanol model mixtures. Molecular dynamics simulation results with the OPLS-AA model for methanol

We have studied thermodynamic and dynamic properties, as well as the dielectric constant, of water–methanol model mixtures in the entire range of composition, by using constant pressure molecular dynamics simulations at room temperature and ambient pressure. The SPC-E and TIP4P-Ew water models are explored in combination with the OPLS-AA (all atom) model of methanol. Changes of the average number of hydrogen bonds between particles of different species and of the fractions of differently bonded molecules are put in correspondence with the evolution of excess mixing properties, of self-diffusion coefficients, reorientational relaxation times and power spectra. We were looking for the existence of correlations between properties of a different kind, in order to find reasons of satisfactory or dissatisfactory performance of the models applied. It seems that any improvement of the description of water–methanol mixtures requires a more sophisticated, carefully parametrized model for methanol and possibly, certain deviations in the combination rules.

(Invited) Bias Potential Controlled First-Principles Calculations in Batteries and Energy Storage Devices

For batteries and energy storage devices, elucidating a microscopic detail of an electrochemical reaction that takes place at the electrode and electrolyte interface is important to improve the performance of the electrochemical devices. Although a first-principles calculation is a realistically feasible way to study the electrochemical reaction, it is still difficult to conduct a direct simulation of those devices in working conditions. There are two major difficulties in practical simulations: First, it is difficult to describe the electric field near the electrode surface which is generated in the electric double layer (EDL); second, it is difficult to control the applied bias potential during the simulation. The former is known to be incompatible with the periodic boundary condition used in the first-principles simulation. The latter requires a technique to control the Fermi energy of the electrode on the fly. In 2006, we proposed a method to incorporate the electric field in the EDL with the periodic boundary condition by solving the Poisson equation using Green’s function method, which is called effective screening medium (ESM) method [1]. Using this technique, we successfully simulated an electrochemical reaction which is driven by the strong electric field at platinum-water interface [2]. To address the second issue, we have developed a simulation scheme for performing a first-principles calculation at a constant electrode potential in which we can control the Fermi energy of the electrode, i.e. bias potential, by connecting the system to a potentiostat [3]. We call our constant electrode potential technique as the “constant-μ” method, where μ denotes the Fermi energy (chemical potential for electrons in a metallic system). The most intuitive way to control the Fermi energy of the system is to fix the Fermi energy to the target Fermi energy (μext) during the self-consistent field (SCF) calculation. But this strategy will fail within a few SCF steps. To realize a stable and robust calculation, we employed the extended system technique which is usually used in the constant temperature and constant pressure molecular dynamics simulation. We assume that the charge in the simulation cell is a dynamical variable; we call this the fictitious charge particle (FCP) [3]. This assumption means that the system is conceptually connected to a charge reservoir and is governed by the following grand potential instead of the total energy: Ω = E tot – μext n FCP, where E tot, μext and n FCP indicate the total energy, the target Fermi energy, and the charge of the FCP. By taking the derivatives of the grand potential Ω with respect to atomic position (τ) and n FCP, we can calculate the generalized forces acting on, respectively, atom, F atom=–dE tot/dτ, and the FCP, F FCP=–(μ–μext), where μ is the instantaneous Fermi energy. Since the deviation of the Fermi energy from the target potential is acting as the restoring force for the FCP, F FCP drives n FCP so that the Fermi energy of the system goes to μext. Once we can define the forces, we can combine the constant-μ method with both dynamical and static first-principles calculations. If one needs to simulate dynamical aspect of the electrochemical system, we can conduct a molecular dynamics (MD) simulation with the constant-μ method. By introducing a constrained dynamics such as the blue-moon ensemble method [4], we can calculate the free energy difference and activation free energy of an electrochemical reaction. In the static calculation, minimization of the grand potential results in a generalized optimization problem with atomic and FCP systems. Besides the geometry optimization process for atoms, we need to optimize n FCP along with F FCP so that the Fermi energy of the system corresponds to the target Fermi energy. If one needs to calculate the grand potential difference and activation energy of a reaction, we can combine the nudged elastic band (NEB) method [5] with the constant-μ method. In the presentation, I will give an overview of the method and show some applications for dynamic and static calculations, such as the solvation/desolvation process of the Li-ion and an atomic diffusion on electrode surface under a bias potential.

Fluctuations in Systems of Hexatic Smectic B Liquid Crystals

Systems of hexatic smectic B liquid crystal are investigated by constant-pressure and constant-temperature molecular dynamics simulations. By using a symplectic integrator designed for soft-matter under a condition where the fluctuation of the system is large, spontaneous transfer among the metastable states is observed. By analyzing the transfer in both directions, distributions of work values corresponding to the forward and reverse processes are obtained.

Origin of Hopping Dynamics in Hexatic Smectic B and Smectic A Liquid Crystals

Origin of hopping dynamics in liquid crystals is investigated by constant-pressure and constant-temperature molecular dynamics simulations. Hopping dynamics is observed in both hexatic smectic B (HexB) and smectic A (SmA) liquid crystals, however in different directions. At the crystal-Hex B transition, a second order transtition occurs where the layer thickness increases continuously. At the HexB-SmA transition, melting inside the layers is observed. The hopping occurs in the direction normal to the direction with loose density, i.e., xy-direction in HexB and z-direction in SmA.

Temperature dependence of the microscopic structure and density anomaly of the SPC/E and TIP4P-Ew water models. Molecular dynamics simulation results

We have investigated temperature trends of the microscopic structure of the SPC/E and TIP4P-Ew water models in terms of the pair distribution functions, coordination numbers, the average number of hydrogen bonds, the distribution of bonding states of a single molecule as well as the angular distribution of molecules by using the constant pressure molecular dynamics simulations. The evolution of the structure is put in correspondence with the dependence of water density on high temperatures down to the region of temperatures where the system becomes supercooled. It is shown that the fraction of molecules with three and four bonds determine the maximum density for both models. Moreover, the temperature dependence of the dielectric constant is obtained and analyzed.

Network Analysis of Free Energy Landscape of Metastable States of Hexatic Smectic B Liquid Crystal

The topology of the free-energy landscape of a model system, which gives rise to multiple metastable states of hexatic smectic B (HexB) liquid crystals, is investigated using network theory. Directed and weighted networks of HexB states are constructed from a series of dynamical data calculated by constant-pressure and constant-temperature molecular dynamics simulations. The k-shell decomposition is extended to directed networks, and the networks of HexB metastable states are analyzed. Singular values of the weighted adjacency matrix, with elements consisting of the weight of the directed edge, are used to distinguish important vertices for evolution.

Melting of icosahedral nickel clusters under hydrostatic pressure.

The thermal stabilities and melting behavior of icosahedral nickel clusters under hydrostatic pressure have been studied by constant-pressure molecular dynamics simulation. The potential energy and Lindemann index are calculated. The overall melting temperature exhibits a strong dependence on pressure. The Lindemann index of solid structure before melting varies slowly and is almost independent of pressure. However, after the clusters melt completely, the Lindemann index at the overall melting point strongly depends on pressure. The overall melting temperature is found to be increasing nonlinearly with increasing pressure, while the volume change during melting decreases linearly with increasing pressure. Under a high pressure and temperature environment, similar angular distributions were found between liquid and solid structures, indicating the existence of a converging local structure.

Confinement-Induced Stiffening of Thin Elastomer Films: Linear and Nonlinear Mechanics vs Local Dynamics

Constant-pressure molecular-dynamics simulations have been carried out of a bead–rod model polymer confined between two attractive crystalline substrates. Three different substrate–substrate separations (i.e., film thicknesses) were used and two different polymer–substrate interaction strengths. The density profiles show a monomer layering close to the polymer–substrate interface. A higher density was found in this region compared to the middle bulklike layers of the films. The dependence of the film-averaged density on temperature and thickness was measured for all polymer films. Decreasing the film thickness leads to an increase of this density and of the film-averaged glass-transition temperature. Layer-resolved analysis of the segmental dynamics of the thickest films shows a gradient of the mobility upon approach of the polymer–substrate interface, while the middle-layer dynamics exhibits bulklike behavior. With decreasing film thickness, these gradients become overlapping. All polymer films were deformed uniaxially normal to the substrates beyond their linear viscoelastic regime; their elastic moduli but also their secant moduli at larger strain amplitudes were extracted. In the linear regime, the stiffness was found to increase with decreasing film thickness; this correlates well with the layer-resolved segmental dynamical behavior. A strong drop of the stiffness was observed upon increasing the deformation amplitude; this drop was more pronounced for stiffer films. It is shown that the drop in stiffness can be qualitatively explained by the drop of the relaxation times as well as by the increased heterogeneity of the dynamics in different film layers upon deformation. The thickness dependencies of the structural, dynamical, and mechanical properties become more pronounced with stronger adsorption to the substrates.

Phase Diagram for a Lennard-Jones System Obtained through Constant-Pressure Molecular Dynamics Simulations

The phase diagram for a Lennard-Jones system was estimated using conventional NPH molecular dynamics (MD) simulations. Standard periodic boundaries were assumed for unit cells containing 1000 molecules. An elongated unit cell with both solid and vacuum sections was found to be suitable for the NPH MD simulations when calculating both the melting and vapor pressure curves under pressures lower than the critical one. Under high pressures, a unit cell with both solid and liquid sections was used as the initial configuration to obtain the melting temperature. The results of these simulations were compared with the phase transition point given by the reported equations of state.

Low concentrated hydroxyectoine solutions in presence of DPPC lipid bilayers: A computer simulation study

The influence of hydroxyectoine on the properties of the aqueous solution in presence of DPPC lipid bilayers is studied via semi-isotropic constant pressure (NPT) Molecular Dynamics simulations. We investigate the solvent–co-solute behavior in terms of Kirkwood–Buff integrals as well as hydrogen bond life times for an increasing hydroxyectoine concentration up to 0.15mol/L. The observed preferential exclusion mechanism identifies hydroxyectoine as a kosmotropic osmolyte. Our findings with regard to the DPPC lipid bilayer indicate an increase of the surface pressure as well as the solvent accessible surface area in presence of higher hydroxyectoine concentrations. The results are in agreement to the outcome of recent experiments. With this study, we are able to validate the visibility of co-solute–solute–solvent effects for low and physiologically relevant osmolyte concentrations.

Simple model for chain packing and crystallization of soft colloidal polymers

We study a simple bead-spring polymer model exhibiting competing crystallization and glass transitions. Constant-pressure molecular dynamics simulations are employed to study phase behavior and morphological order. For adequately slow quench rates, chain systems exhibit a first-order phase transition (crystallization) below a critical temperature T=T(cryst). We observe the formation of close-packed crystallites of FCC and/or HCP order, separated by domain walls, twin defects, and amorphous regions. Such crystal structures closely resemble the corresponding ordered morphologies of athermal polymer packings: fully flexible chains retain random-walk-like configurations in the crystalline state and do not form lamellae, while semiflexible chains do form lamellae. The model presented here is well suited to the modeling of granular and colloidal polymers, in particular for elucidating the factors that dictate the formation of specific ordered morphologies.

Local order evolution of liquid Cu during glass transition under different pressures: A molecular dynamics study

Based on the second-moment approximation of tight-binding scheme, constant-pressure molecular dynamics simulations are performed for liquid Cu during the glass transition under different pressures. By means of pair analysis technique and bond orientational order analysis we find that the dominant bond pairs are those related to fcc and hcp crystalline order not those representing icosahedral short-range order (ISRO) when the systems enter into glass transition region. Although these two kinds of bond pairs compete with each other, the system tends towards a mixture of crystalline bond pairs during glass formation. The effect on various bond pairs brought about by higher pressure is much less for liquids than for glasses. The experimental observation of a shoulder on the second peak of the structure factor for supercooled liquids might not merely attribute to ISRO, since supercooled liquid Cu exhibits such a shoulder, but does not display an enhanced icosahedral symmetry.

Pressure-induced structural transition of ZnO nanocrystals studied with molecular dynamics

We have studied the pressure-induced structural transition of ZnO nanocrystals using constant pressure molecular dynamics simulations for finite system. We have observed the transition from the fourfold coordination wurtzite to the sixfold coordination rocksalt structure, and the process of transition is strongly dependent on the morphology of the nanocrystals. It is found that the perfect faceted ZnO nanocrystals undergo wurtzite to rocksalt transition with a perfect fivefold h-MgO structure as the intermediate status. But for the faceted ones without perfect surface structure, as the number of the atoms removed from the (001) and (001¯) surface edge increases, the local morphology will become more similar to spherical. The nanocrystal will receive equal stress from every direction and it will be more difficult to compress the structure along only c axis as the perfect faceted ZnO nanocrystal. In this situation, only partial structure experiences intermediate fivefold coordination structure or even no intermediate fivefold coordination structure exists dependent on the surface disorder level.