Hysteresis Elimination for an Anisotropic Liquid-Crystal Model via Molecule Design and Replica-Exchange Optimization.

The phenomenon of hysteresis in simulations, in which a system's current state is correlated to previous states and inhibits the transition to a more stable phase, may often lead to misleading results in physical chemistry. In this study, in addition to the replica exchange method (REM), a novel approach was taken by combining an evolution strategy based on the evolutionary principles of nature to predict phase transitions for the Hess-Su liquid-crystal model. In this model, an anisotropy term is added to the simple 6-12 Lennard-Jones model to intuitively reproduce the behavior of liquid crystals. We first applied the pressure-temperature REM to the Hess-Su model and optimized the replica spacing for the energy distribution to gain the maximum advantage from the REM. We then used the same approach as for the Hamiltonian REM, seeking to optimize the replica spacing in the same way. Based on both results, we attempted to predict this coarse-grained liquid-crystal model's exact phase transition point. In the Hamiltonian REM, replicas were prepared with different molecular aspect ratios corresponding to the values of the anisotropy terms in the potential function. The Hess-Su liquid-crystal model, which undergoes a direct transition from the nematic to the solid phase without going through a smectic phase, is a challenging research target for understanding phase transitions. Despite the tremendous computational difficulty in overcoming the strong hysteresis present in this system, our method could predict the phase transition point clearly and significantly reduce the extent of hysteresis. Our approach is beneficial when simulating more complex systems and, above all, shows great potential for more accurate and efficient phase transition predictions in the field of molecular simulation in the future.

Novel approach for designing order parameters of clathrate hydrate structures by graph neural network.

Clathrate hydrates continue to be the focus of active research efforts due to their use in energy resources, transportation, and storage-related applications. Therefore, it is crucial to define their essential characteristics from a molecular standpoint. Understanding molecular structure in particular is crucial because it aids in understanding the mechanisms that lead to the formation or dissociation of clathrate hydrates. In the past, a wide variety of order parameters have been employed to classify and evaluate hydrate structures. An alternative approach to inventing bespoke order parameters is to apply machine learning techniques to automatically generate effective order parameters. In earlier work, we suggested a method for automatically designing novel parameters for ice and liquid water structures with Graph Neural Networks (GNNs). In this work, we use a GNN to implement our method, which can independently produce feature representations of the molecular structures. By using the TeaNet-type model in our method, it is possible to directly learn the molecular geometry and topology. This enables us to build novel parameters without prior knowledge of suitable order parameters for the structure type, discover structural differences, and classify molecular structures with high accuracy. We use this approach to classify the structures of clathrate hydrate structures: sI, sII, and sH. This innovative approach provides an appealing and highly accurate replacement for the traditional order parameters. Furthermore, our method makes clear the process of automatically designing a universal parameter for liquid water, ice, and clathrate hydrate to analyze their structures and phases.

Prediction of Water Diffusion in Wide Varieties of Polymers with All-Atom Molecular Dynamics Simulations and Deep Generative Models.

Permeation through polymer membranes is an important technology in the chemical industry, and in its design, the self-diffusion coefficient is one of the physical quantities that determine permeability. Since the self-diffusion coefficient sensitively reflects intra- and intermolecular interactions, analysis using an all-atom model is required. However, all-atom simulations are computationally expensive and require long simulation times for the diffusion of small molecules dissolved in polymers. MD-GAN, a machine learning model, is effective in accelerating simulations and reducing computational costs. The target systems for MD-GAN prediction were limited to polyethylene melts in previous studies; therefore, this study extended MD-GAN to systems containing copolymers with branches and successfully predicted water diffusion in various polymers. The correlation coefficient between the predicted self-diffusion coefficient and that of the long-time simulation was 1.00. Additionally, we found that incorporating statistical domain knowledge into MD-GAN improved accuracy, reducing the mean-square displacement prediction outliers from 14.6% to 5.3%. Lastly, the distribution of latent variables with embedded dynamics information within the model was found to be strongly related to accuracy. We believe that these findings can be useful for the practical applications of MD-GAN.

Graph neural networks classify molecular geometry and design novel order parameters of crystal and liquid.

Molecular dynamics simulation produces three-dimensional data on molecular structures. The classification of molecular structure is an important task. Conventionally, various order parameters are used to classify different structures of liquid and crystal. Recently, machine learning (ML) methods have been proposed based on order parameters to find optimal choices or use them as input features of neural networks. Conventional ML methods still require manual operation, such as calculating the conventional order parameters and manipulating data to impose rotational/translational invariance. Conversely, deep learning models that satisfy invariance are useful because they can automatically learn and classify three-dimensional structural features. However, in addition to the difficulty of making the learned features explainable, deep learning models require information on large structures for highly accurate classification, making it difficult to use the obtained parameters for structural analysis. In this work, we apply two types of graph neural network models, the graph convolutional network (GCN) and the tensor embedded atom network (TeaNet), to classify the structures of Lennard-Jones (LJ) systems and water systems. Both models satisfy invariance, while GCN uses only length information between nodes. TeaNet uses length and orientation information between nodes and edges, allowing it to recognize molecular geometry efficiently. TeaNet achieved a highly accurate classification with an extremely small molecular structure, i.e., when the number of input molecules is 17 for the LJ system and 9 for the water system, the accuracy is 98.9% and 99.8%, respectively. This is an advantage of our method over conventional order parameters and ML methods such as GCN, which require a large molecular structure or the information of wider area neighbors. Furthermore, we verified that TeaNet could build novel order parameters without manual operation. Because TeaNet can recognize extremely small local structures with high accuracy, all structures can be mapped to a low-dimensional parameter space that can explain structural features. TeaNet offers an alternative to conventional order parameters because of its novelty.

Construction of full-atomistic polymer amorphous structures using reverse-mapping from Kremer-Grest models.

We propose a method to build full-atomistic (FA) amorphous polymer structures using reverse-mapping from coarse-grained (CG) models. In this method, three models with different resolutions are utilized, namely the CG1, CG2, and FA models. It is assumed that the CG1 model is more abstract than the CG2 model. The CG1 is utilized to equilibrate the system, and then sequential reverse-mapping procedures from the CG1 to the CG2 models and from the CG2 to the FA models are conducted. A mapping relation between the CG1 and the FA models is necessary to generate a polymer structure with a given density and radius of chains. Actually, we have used the Kremer-Grest (KG) model as the CG1 and the monomer-level CG model as the CG2 model. Utilizing the mapping relation, we have developed a scheme that constructs an FA polymer model from the KG model. In the scheme, the KG model, the monomer level CG model, and the FA model are successively constructed. The scheme is applied to polyethylene (PE), cis 1,4-polybutadiene (PB), and poly(methyl methacrylate) (PMMA). As a validation, the structures of PE and PB constructed by the scheme were carefully checked through comparison with those obtained using long-time FA molecular dynamics (MD) simulations. We found that both short- and long-range chain structures constructed by the scheme reproduced those obtained by the FA MD simulations. Then, as an interesting application, the scheme is applied to generate an entangled PMMA structure. The results showed that the scheme provides an efficient and easy way to construct amorphous structures of FA polymers.

Impact of free energy of polymers on polymorphism of polymer-grafted nanoparticles.

Colloidal crystals have gathered wide attention as a model material for optical applications because of their feasibility in controlling the propagation of light by their crystal structure and lattice spacing as well as the simplicity of their fabrication. However, due to the simple interaction between colloids, the colloidal crystal structures that can be formed are limited. It is also difficult to adjust the lattice spacing. Furthermore, colloidal crystals are fragile compared to other crystals. In this study, we focused on polymer-grafted nanoparticles (PGNP) as a possible solution to these unresolved issues. We expected that PGNPs, composed of two distinct layers (the hard core of a nanoparticle and the soft corona of grafted polymers on the surface), will demonstrate similar behaviors as star polymers and hard spheres. We also predicted that PGNPs may exhibit polymorphism because the interaction between PGNPs strongly depends upon their grafting density and the length of the grafted polymer chains. Moreover, we expected that crystals made from PGNPs will be structurally tough due to the entanglement of grafted polymers. From exploration of crystal polymorphs of PGNPs by molecular dynamics simulations, we found face-centered cubic (FCC)/hexagonal close-packed (HCP) and body-centered cubic (BCC) crystals, depending on the length of the grafted polymer chains. When the chains were short, PGNPs behaved like hard spheres and crystals were arranged in FCC/HCP structure, much like the phase transition observed in an Alder transition. When the chains were long enough, the increase in the free energy of grafted polymers was no longer negligible and crystals were arranged in BCC structure, which has a lower density than FCC/HCP. When the chains were not too short or long, FCC/HCP structures were first observed when the volume fraction of system was small, but a phase transition occurred when the system was further compressed and the crystals arranged themselves in a BCC structure. These results most likely have laid strong foundations for future simulations and experimental studies of PGNP crystals.

MD-GAN with multi-particle input: the machine learning of long-time molecular behavior from short-time MD data.

Molecular dynamics simulation is a method of investigating the behavior of molecules, which is useful for analyzing a variety of structural and dynamic properties and mechanisms of phenomena. However, the huge computational cost of large-scale and long-time simulations is an enduring problem that must be addressed. MD-GAN is a machine learning-based method that can evolve part of the system at any time step, accelerating the generation of molecular dynamics data [Endo et al., Proceedings of the AAAI Conference on Artificial Intelligence, 2018, 32]. For the accurate prediction of MD-GAN, sufficient information on the dynamics of a part of the system should be included with the training data. Therefore, the selection of the part of the system is important for efficient learning. In a previous study, only one particle (or vector) of each molecule was extracted as part of the system. The effectiveness of adding information from other particles to the learning process is investigated in this study. When the dynamics of three particles of each molecule were used in the polyethylene experiment, the diffusion was successfully predicted using the training data with a time length of approximately 40%, compared to the single-particle input. Surprisingly, the unobserved transition of diffusion in the training data was also predicted using this method. The reduced cost for the generation of training MD data achieved in this study is useful for accelerating MD-GAN.

An Efficient Random Number Generation Method for Molecular Simulation.

We propose a new random number generation method, which is the fastest and the simplest of its kind, for use with molecular simulation. We also discuss the possibility of using this method with various other numerical calculations. To demonstrate the significant increases in calculation speeds that can be gained by using our method, we present a comparison with prior methods for dissipative particle dynamics (DPD) simulations. The DPD method uses random numbers to reproduce thermal fluctuations of molecules. As such, an efficient method to generate random numbers in parallel computing environments has been widely sought after. Several random number generation methods have been developed that use encryption. In this study, we establish for the first time that random numbers with desirable properties exist in the particle coordinates used in DPD calculations. We propose a method for generating random numbers without encryption that utilizes this source of randomness. This is an innovative method with minimal computational cost, since it is not dependent on a complicated random number generation algorithm or an encryption process. Furthermore, our method may lead to faster random number generation for many other physical and chemical simulations.

Optimal Replica-Exchange Molecular Simulations in Combination with Evolution Strategies.

We have incorporated Evolution Strategies into the Replica-Exchange Monte Carlo simulation method to predict the phase behavior of several example fluids. The replica-exchange method allows one system to exchange temperatures with its neighbors to search for the most stable structure relatively efficiently in a single simulation. However, if the temperature intervals of the replicas are not positioned carefully, there is an issue that local exchange does not occur. Our results for a simple Lennard-Jones fluid and the liquid-crystal Yukawa model demonstrate the utility of the approach when compared to conventional methods. When Evolution Strategies were applied to the Replica-Exchange Monte Carlo simulation, the problem of a significant localized decrease in exchange probability near the phase transition was avoided. By obtaining the optimal temperature intervals, the system efficiently traverses a broader parameter space with a small number of replicas. This is equivalent to accelerating molecular simulations with limited computational resources and can be useful when attempting to predict the phase behavior of complex systems.

Correlation between ordering and shear thinning in confined OMCTS liquids.

Despite decades of extensive research, the behavior of confined liquids, particularly in the mixed/boundary lubrication regime, remains unelucidated. This can be attributed to several factors, including the difficulty to make direct experimental observations of the behavior of lubricant molecules under nonequilibrium conditions, the high computational cost of molecular simulations to reach steady state, and the low signal-to-noise ratio at extremely low shear rates corresponding to actual operating conditions. In this regard, we studied the correlation between the structure formation and shear viscosity of octamethylcyclotetrasiloxane confined between two mica surfaces in a mixed/boundary lubrication regime. Three different surface separations-corresponding to two-, three-, and five-layered structures-were considered to analyze the effect of confinement. The orientational distributions with one specific peak for n = 2 and two distributions, including a parallel orientation with the surface normal for n > 2, were observed at rest. The confined liquids exhibited a distinct shear-thinning behavior independent of surface separations for a relatively low shear rate, γ̇≲108s-1. However, the shear viscosities at γ̇≲108s-1 depended on the number of layered structures. Newtonian behavior was observed with further increase in the shear rate. Furthermore, we found a strong correlation between the degree of molecular orientation and the shear viscosity of the confined liquids. The magnitude of the shear viscosity of the confined liquids can primarily be determined by the degree of molecular orientation, and shear thinning originates from the vanishing of specific orientational distributions with increasing shear rate.

Water molecules in CNT-Si3N4 membrane: Properties and the separation effect for water-alcohol solution.

Water confined in carbon nanotubes (CNTs) has been intensively studied because of its unique properties and potential for various applications and is often embedded in silicon nitride (Si3N4) membranes. However, the understanding of the influence of Si3N4 on the properties of water in CNTs lacks clarity. In this study, we performed molecular dynamics simulations to investigate the effect of the Si3N4 membrane on water molecules inside CNTs. The internal electric field generated in the CNTs by the point charges of the Si3N4 membrane changes the structure and dynamical properties of water in the nanotubes, causing it to attain a disordered structure. The Si3N4 membrane decreases the diffusivity of water in the CNTs; this is because the Coulomb potential energy (i.e., electrostatic interaction) of water decreases owing to the presence of Si3N4, whereas the Lennard-Jones potential energy (i.e., van der Waals interaction) does not change significantly. Furthermore, electrostatic interactions make the water structure more stable in the CNTs. As a result, the Si3N4 membrane enhances the separation effect of the water-methanol mixture with CNTs in the presence of an external electric field. Furthermore, the threshold of the external electric field strength to induce water-methanol separation with CNTs is reduced owing to the presence of a silicon nitride membrane.

Phase behaviors of deeply supercooled bilayer water unseen in bulk water.

Akin to bulk water, water confined to an isolated nanoslit can show a wealth of new 2D phases of ice and amorphous ice, as well as unusual phase behavior. Indeed, 2D water phases, such as bilayer hexagonal ice and monolayer square ice, have been detected in the laboratory, confirming earlier computational predictions. Herein, we report theoretical evidence of a hitherto unreported state, namely, bilayer very low density amorphous ice (BL-VLDA), as well as evidence of a strong first-order transition between BL-VLDA and the BL amorphous ice (BL-A), and a weak first-order transition between BL-VLDA and the BL very low density liquid (BL-VLDL) water. The diffusivity of BL-VLDA is typically in the range of 10-9 cm2/s to 10-10 cm2/s. Similar to bulk (3D) water, 2D water can exhibit two forms of liquid in the deeply supercooled state. However, unlike supercooled bulk water, for which the two forms of liquid can coexist and merge into one at a critical point, the 2D BL-VLDL and BL high-density liquid (BL-HDL) phases are separated by the highly stable solid phase of BL-A whose melting line exhibits the isochore end point (IEP) near 220 K in the temperature-pressure diagram. Above the IEP temperature, BL-VLDL and BL-HDL are indistinguishable. At negative pressures, the metastable BL-VLDL exhibits a spatially and temporally heterogeneous structure induced by dynamic changes in the nanodomains, a feature much less pronounced in the BL-HDL.

Phase Transitions and Hysteresis for a Simple Model Liquid Crystal by Replica-Exchange Monte Carlo Simulations.

In this work, the advantages of applying the temperature and pressure replica-exchange method to investigate the phase transitions and the hysteresis for liquid-crystal fluids were demonstrated. In applying this method to the commonly used Hess-Su liquid-crystal model, heat capacity peaks and points of phase co-existence were observed. The absence of a smectic phase at higher densities and a narrow range of the nematic phase were reported. The identity of the crystalline phase of this system was found to a hexagonal close-packed solid. Since the nematic-solid phase transition is strongly first order, care must be taken when using this model not to inadvertently simulate meta-stable nematic states at higher densities. In further analysis, the Weighted Histogram Analysis Method was applied to verify the precise locations of the phase transition points.

Mechanism for H2 diffusion in sII hydrates by molecular dynamics simulations.

Among the many different types of molecules that form clathrate hydrates, H2 is unique as it can easily diffuse into and out of clathrate cages, a process that involves the physical-chemical interactions between guest (H2) and host (water) molecules, and is unlike any other molecular system. The dynamic and nano-scale process of H2 diffusion into binary structure II hydrates, where the large cages are occupied by larger molecules, was studied using molecular dynamics simulation. As the H2 molecules diffused from one cage to another, two types of diffusion processes were observed: (i) when moving between a pair of large cages, the H2 molecules pass through the central part of the hexagonal rings; (ii) however, when the H2 molecules move from a large cage to a small one, it requires one of the pentagonal rings to partially break, as this allows the H2 molecule to pass through the widened space. While the diffusion of H2 molecules between large cages was found to occur more frequently, the presence of SF6 molecules in the large cages was found to inhibit diffusion. Therefore, in order to attain higher H2 storage capacities in binary hydrates, it is suggested that there is an optimal number of large cages that should be occupied by SF6 molecules.

Evidence of low-density and high-density liquid phases and isochore end point for water confined to carbon nanotube.

Possible transition between two phases of supercooled liquid water, namely the low- and high-density liquid water, has been only predicted to occur below 230 K from molecular dynamics (MD) simulation. However, such a phase transition cannot be detected in the laboratory because of the so-called "no-man's land" under deeply supercooled condition, where only crystalline ices have been observed. Here, we show MD simulation evidence that, inside an isolated carbon nanotube (CNT) with a diameter of 1.25 nm, both low- and high-density liquid water states can be detected near ambient temperature and above ambient pressure. In the temperature-pressure phase diagram, the low- and high-density liquid water phases are separated by the hexagonal ice nanotube (hINT) phase, and the melting line terminates at the isochore end point near 292 K because of the retracting melting line from 292 to 278 K. Beyond the isochore end point (292 K), low- and high-density liquid becomes indistinguishable. When the pressure is increased from 10 to 600 MPa along the 280-K isotherm, we observe that water inside the 1.25-nm-diameter CNT can undergo low-density liquid to hINT to high-density liquid reentrant first-order transitions.

Separation of water-alcohol mixtures using carbon nanotubes under an electric field.

Carbon nanotubes (CNTs) are a promising candidate for separation membranes because of their ability to transport substances at very high flow rates. However, there is a tradeoff between achieving a high selectivity using small pore sizes and the reduction of water flux. Here, using molecular dynamics simulations, we report that CNTs can effectively separate water-methanol mixtures under an electric field. Without an electric field and under piston pressure, both water and methanol flow through a CNT, resulting in no separation effect. In contrast, under an electric field and high piston pressure, CNTs allow selective water permeation while rejecting the permeation of methanol molecules. This separation effect is caused by the ordered structures of water molecules in the CNT. A high filtering effect is observed under the conditions of high methanol concentration in the solution or even with large-diameter CNTs up to 3.39 nm. As long as the ordered structure of water in the CNTs can be maintained, the strong filtering effect can be maintained.

Ordering in clusters of uniaxial anisotropic particles during homogeneous nucleation and growth.

The nucleation process of anisotropic particles often differs from that of their spherically symmetric counterparts. Despite a large body of work on the structure of droplets of anisotropic particles, their formation process remains poorly understood. In this study, homogeneous nucleation of uniaxial anisotropic particles was studied. Through structural analysis of cluster development and the formation free energy during the nucleation stage, it was revealed that the nucleation of uniaxial particles begins from highly ordered states. There is, however, a marked decrease in orientational order within the cluster before critical nucleus size is attained. Further investigation on variations in the molecular interactions demonstrates how droplet elongation and the direction of the nematic ordering director relative to the axis of elongation can both be controlled according to the nature of the molecular anisotropy.

Cage occupancies, lattice constants, and guest chemical potentials for structure II hydrogen clathrate hydrate from Gibbs ensemble Monte Carlo simulations.

In this paper, equilibrium properties of structure II hydrates of hydrogen were determined from Monte Carlo simulations in the isothermal-isobaric Gibbs ensemble. Water and hydrogen molecules are described by the TIP4P/Ice and Silvera-Goldman models, respectively. The use of the Gibbs ensemble has many key advantages for the simulation of hydrates. By the separation of hydrogen vapor and hydrate phases into their own domains, coupled with transfer moves of hydrogen molecules between domains, cage occupancies were determined. Furthermore, the choice of this ensemble also allows equilibrium lattice constants and guest molecule chemical potentials to be straightforwardly estimated. Results for hydrogen mass fractions indicate reasonable agreement with prior simulation data and theoretical models, while detailed analysis of cage occupancy distributions and neighboring cage pair occupancy combinations gives valuable insight into the behavior of this hydrate at the inter-cage scale. These results will aid in the construction of theoretical models, for which knowledge of the occupancy of neighboring cages is of great importance. In support of previous experimental and theoretical works, we also find evidence of double occupancy of a few small cages inside of the hydrate stability zone, albeit at very high pressures; approximately 0.1% of small cages are doubly occupied at 300 MPa, for temperatures of 225 K and 250 K.

Molecular Dynamics Simulation of Ice Crystal Growth Inhibition by Hexadecyl-trimethyl-ammonium Bromide.

Recent experiments have found hexadecyl-trimethyl-ammonium bromide (CTAB) to have superior ice nucleation inhibition properties [ J. Phys. Chem. B 121, 6580]. The mechanism of how the inhibition takes place remains unclear. Therefore, molecular dynamics was used to simulate ice crystallization of a water/CTAB/ice system. The ice crystallization rate for a pure water system was compared for the basal [0001], first prism [101̅0], and secondary prism plane [112̅0], where the basal plane grew the slowest followed by the first prism plane. When CTAB was added to the ice-liquid water system, crystallization was clearly impeded. Even when ice starts growing away from the CTAB molecule, the hydrophilic head would at some point protrude and get caught in the water/ice interface. Once the head of the CTAB was encapsulated in the advancing interface, the hydrophobic body would wriggle around and disrupt the formation of hydrogen bond networks that are essential for ice growth. When the interface clears the length of the body of the CTAB molecule, ice crystallization resumes at its normal pace. In summary, the inhibition of ice growth is a combination of the hydrophilic head acting as an anchor and the dynamic motion of the hydrophobic tail hindering stable hydrogen bonding for ice growth.

Structural determinants in the bulk heterojunction.

Photovoltaics is one of the key areas in renewable energy research with remarkable progress made every year. Here we consider the case of a photoactive material and study its structural composition and the resulting consequences for the fundamental processes driving solar energy conversion. A multiscale approach is used to characterize essential molecular properties of the light-absorbing layer. A selection of bulk-representative pairs of donor/acceptor molecules is extracted from the molecular dynamics simulation of the bulk heterojunction and analyzed at increasing levels of detail. Significantly increased ground state energies together with an array of additional structural characteristics are identified that all point towards an auxiliary role of the material's structural organization in mediating charge-transfer and -separation. Mechanistic studies of the type presented here can provide important insights into fundamental principles governing solar energy conversion in next-generation photovoltaic devices.