An exa-scale high-performance molecular dynamics simulation program: MODYLAS.

A new version of the highly parallelized general-purpose molecular dynamics (MD) simulation program MODYLAS with high performance on the Fugaku computer was developed. A benchmark test using Fugaku indicated highly efficient communication, single instruction, multiple data (SIMD) processing, and on-cache arithmetic operations. The system's performance deteriorated only slightly, even under high parallelization. In particular, a newly developed minimum transferred data method, requiring a significantly lower amount of data transfer compared to conventional communications, showed significantly high performance. The coordinates and forces of 101 810 176 atoms and the multipole coefficients of the subcells could be distributed to the 32 768 nodes (1 572 864 cores) in 2.3 ms during one MD step calculation. The SIMD effective instruction rates for floating-point arithmetic operations in direct force and fast multipole method (FMM) calculations measured on Fugaku were 78.7% and 31.5%, respectively. The development of a data reuse algorithm enhanced the on-cache processing; the cache miss rate for direct force and FMM calculations was only 2.74% and 1.43%, respectively, on the L1 cache and 0.08% and 0.60%, respectively, on the L2 cache. The modified MODYLAS could complete one MD single time-step calculation within 8.5 ms for the aforementioned large system. Additionally, the program contains numerous functions for material research that enable free energy calculations, along with the generation of various ensembles and molecular constraints.

Algorithm to minimize MPI communications in the parallelized fast multipole method combined with molecular dynamics calculations.

In the era of exascale supercomputers, large-scale, and long-time molecular dynamics (MD) calculations are expected to make breakthroughs in various fields of science and technology. Here, we propose a new algorithm to improve the parallelization performance of message passing interface (MPI)-communication in the MPI-parallelized fast multipole method (FMM) combined with MD calculations under three-dimensional periodic boundary conditions. Our approach enables a drastic reduction in the amount of communication data, including the atomic coordinates and multipole coefficients, both of which are required to calculate the electrostatic interaction by using the FMM. In communications of multipole coefficients, the reduction rate of communication data in the new algorithm relative to the amount of data in the conventional one increases as both the number of FMM levels and the number of MPI processes increase. The aforementioned rate increase could exceed 50% as the number of MPI processes becomes larger for very large systems. The proposed algorithm, named the minimum-transferred data (MTD) method, should enable large-scale and long-time MD calculations to be calculated efficiently, under the condition of massive MPI-parallelization on exascale supercomputers.

Molecular dynamics study of lipid bilayers modeling outer and inner leaflets of plasma membranes of mouse hepatocytes. I. Differences in physicochemical properties between the two leaflets.

Outer and inner leaflets of plasma cell membranes have different lipid compositions, and the membrane properties of each leaflet can differ from each other significantly due to these composition differences. However, because of the experimental difficulty in measuring the membrane properties for each leaflet separately, the differences are not well understood at a molecular level. In this study, we constructed two lipid bilayer systems, modeling outer and inner leaflets of plasma membranes of mouse hepatocytes based on experimental composition data. The ion concentration in the interlamellar water phase was also set to match the concentration in extra- and intracellular fluids. The differences in physical properties between the outer and inner leaflets of mouse hepatocyte cell membrane models were investigated by performing 1.2 µs-long all-atomistic molecular dynamics calculations under physiological temperature and pressure conditions (310.15 K and 1 atm). The calculated electron density profiles along the bilayer normal for each model bilayer system captured well the asymmetric feature of the experimental electron density profile across actual cell plasma membranes, indicating that our procedure of modeling the outer and inner leaflets of the cell plasma membranes was satisfactory. We found that compared to the outer leaflet model, the inner leaflet model had a very bulky and soft structure in the lateral direction. To confirm the differences, membrane fluidity was measured from the lateral diffusivity and relaxation times. The fluidity was significantly higher in the inner leaflet model than in the outer leaflet model. We also discuss two topics that are of wide interest in biology, i.e., the interdigitation of acyl tails of lipid molecules between two monolayers and the lateral concentration fluctuation of lipid molecules in the bilayers.

Extension of the fast multipole method for the rectangular cells with an anisotropic partition tree structure.

The fast multipole method (FMM) is an order N method for the numerically rigorous calculation of the electrostatic interactions among point charges in a system of interest. The FMM is utilized for massively parallelized software for molecular dynamics (MD) calculations. However, an inconvenient limitation is imposed on the implementation of the FMM: In three-dimensional case, a cubic MD unit cell is hierarchically divided by the octree partitioning under isotropic periodic boundary conditions along three axes. Here, we extended the FMM algorithm adaptive to a rectangular MD unit cell with different periodicity along the axes by applying an anisotropic hierarchical partitioning. The algorithm was implemented into the parallelized general-purpose MD calculation software designed for a system with uniform distribution of point charges in the unit cell. The partition tree can be a mixture of binary and ternary branches, the branches being chosen arbitrarily with respect to the coordinate axes at any levels. Errors in the calculated electrostatic interactions are discussed in detail for a selected partition tree structure. The extension enables us to execute MD calculations under more general conditions for the shape of the unit cell, partition tree, and boundary conditions, keeping the accuracy of the calculated electrostatic interactions as high as that with the conventional FMM. An extension of the present FMM algorithm to other prime number branches, such as 5 and 7, is straightforward.

Fast multipole method for three-dimensional systems with periodic boundary condition in two directions.

We derived a new expression for the electrostatic interaction of three-dimensional charge-neutral systems with two-dimensional periodic boundary conditions (slab geometry) using a fast multipole method (FMM). Contributions from all the image cells are expressed as a sum of real and reciprocal space terms, and a self-interaction term. The reciprocal space contribution consists of two parts: zero and nonzero terms of the absolute value of the reciprocal lattice vector. To test the new expressions, electrostatic interactions were calculated for a randomly placed charge distribution in a cubic box and liquid water produced by molecular dynamics calculation. The accuracy could be controlled by the degree of expansion of the FMM. In the present expression, the computational complexity of the electrostatic interaction of N-particle systems is order N, which is superior to that of the conventional two-dimensional periodic Ewald method for a slab geometry and the particle mesh Ewald method with a large empty space at an interface of the unit cell. © 2020 Wiley Periodicals, Inc.

Molecular Behavior of Linear Alkylbenzene Sulfonate in Hydrated Crystal, Tilted Gel, and Liquid Crystal Phases Studied by Molecular Dynamics Simulation.

The lamellar phase produced by surfactants with water exhibits several subphases, such as hydrated crystal (Lc), gel (Lß), tilted gel (Lß'), and liquid crystal (Lα) phases, depending on temperature, pressure, and hydration. The dynamics of the surfactant molecules in these phases are still unclear. In the present study, we investigate the translational and conformational dynamics of sodium linear alkylbenzene sulfonate (LAS) molecules in the Lc, Lß', and Lα phases. In the Lα phase, the lateral diffusion of LAS is as fast as that found for phospholipid bilayers in the Lα phase. The diffusion coefficient was undetectably small for the Lc and Lß' phases. The conformation of LAS in the Lα phase relaxes very rapidly, whereas those in the Lc and Lß' phases relax very slowly. The time scale of the relaxations greatly depends on the segment of the LAS molecule for the latter two phases. The relaxation time for the SO3- head group and benzene ring of LAS was much longer than that for alkyl chains. Conformational pattern analyses of LAS alkyl chains revealed that the high fraction of the gauche conformation for the odd-numbered C-C bonds aligns the chain parallel to the bilayer normal and is the main origin of the different relaxation times for different segments in the chain. In the Lc, Lß', and Lα phases, the orientations of the SO3- group and the benzene ring are locked by the salt bridge among SO3- groups and sodium ions. As a result, the orientational order found for the C-C bonds in the LAS alkyl chains is kept even in the Lα phase.

Structure of Hydrated Crystal (Lc), Tilted Gel (Lß'), and Liquid Crystal (Lα) Phases of Linear Alkylbenzene Sulfonate (LAS) Studied by X-ray Diffraction and Molecular Dynamics Simulation.

Precipitation of crystals in aqueous surfactant solutions hampers their use as liquid detergents. Understanding the molecular structure that causes the crystal formation is of key importance in designing a new surfactant. In this study, we first identified the hydrated crystal (Lc) structure of linear alkylbenzene sulfonate (LAS) at ambient temperature by combining X-ray diffraction (XRD) with all-atom molecular dynamics (MD) simulation. It is known that the Lc domain involves only a stoichiometrically small number of water molecules, and the water number of existing interlamellar layers may change when the system shows the thermally induced phase transition. Therefore, we explored possible crystal structures with five different hydration levels of 0, 0.5, 1, 2, and 4 water per LAS using MD simulations with the d-spacing measured by XRD. Among them, only the diffraction pattern calculated for the monohydrate LAS well coincided with the experimental diffraction pattern at 300 K. A structural change of the monohydrate from the Lc phase to the tilted gel (Lß') phase was also observed by heating from 300 to 360 K in the MD simulations, where a cross-sectional structure of the alkyl groups of LAS molecules changed from a face-centered rectangular lattice to a hexagonal one. Further heating to 400 K resulted in a disordered liquid crystal (Lα) phase. Thus, configurational changes of the detergent during thermally induced phase transitions were well explored and characterized by MD simulations with the aid of XRD data. This approach can be broadly applied where we explore the phase behavior of hydrated crystals of any surfactant species.

Pressure tensor for electrostatic interaction calculated by fast multipole method with periodic boundary condition.

A microscopic expression of the pressure tensor using the fast multipole method (FMM) with periodic boundary conditions has been derived. The pressure tensor calculated using this expression has been compared with that obtained using the Ewald method with high accuracy. The precision of the pressure tensor can be controlled as a function of expansion order p of FMM. Using the calculated pressure tensor, the constant pressure molecular dynamics calculation with fully fluctuating cell can be performed for anisotropic systems such as crystals, metals, liquid crystals, glasses, polymers, lipid bilayers, and interfacial regions between two phases. © 2018 Wiley Periodicals, Inc.

Evaluation of atomic pressure in the multiple time-step integration algorithm.

In molecular dynamics (MD) calculations, reduction in calculation time per MD loop is essential. A multiple time-step (MTS) integration algorithm, the RESPA (Tuckerman and Berne, J. Chem. Phys. 1992, 97, 1990-2001), enables reductions in calculation time by decreasing the frequency of time-consuming long-range interaction calculations. However, the RESPA MTS algorithm involves uncertainties in evaluating the atomic interaction-based pressure (i.e., atomic pressure) of systems with and without holonomic constraints. It is not clear which intermediate forces and constraint forces in the MTS integration procedure should be used to calculate the atomic pressure. In this article, we propose a series of equations to evaluate the atomic pressure in the RESPA MTS integration procedure on the basis of its equivalence to the Velocity-Verlet integration procedure with a single time step (STS). The equations guarantee time-reversibility even for the system with holonomic constrants. Furthermore, we generalize the equations to both (i) arbitrary number of inner time steps and (ii) arbitrary number of force components (RESPA levels). The atomic pressure calculated by our equations with the MTS integration shows excellent agreement with the reference value with the STS, whereas pressures calculated using the conventional ad hoc equations deviated from it. Our equations can be extended straightforwardly to the MTS integration algorithm for the isothermal NVT and isothermal-isobaric NPT ensembles. © 2017 Wiley Periodicals, Inc.

Detailed Structural Analysis of a Self-Assembled Vesicular Amphiphilic NCN-Pincer Palladium Complex by Using Wide-Angle X-Ray Scattering and Molecular Dynamics Calculations.

Wide-angle X-ray scattering experiments and all-atomistic molecular dynamics calculations were performed to elucidate the detailed structure of bilayer vesicles constructed by self-assembly of an amphiphilic palladium NCN-pincer complex. We found an excellent agreement between the experimental and calculated X-ray spectra, and between the membrane thickness determined from a TEM image and that calculated from an electron-density profile, which indicated that the calculated structure was highly reliable. The analysis of the simulated bilayer structure showed that in general the membrane was softer than other phospholipid bilayer membranes. In this bilayer assemblage, the degree of alignment of complex molecules in the bilayer membrane was quite low. An analysis of the electron-density profile shows that the bilayer assemblage contains a space through which organic molecules can exit. Furthermore, the catalytically active center is near this space and is easily accessible by organic molecules, which permits the bilayer membrane to act as a nanoreactor. The free energy of permeation of water through the bilayer membrane of the amphiphilic complex was 12 kJ mol-1 , which is much lower than that for phospholipid bilayer membranes in general. Organic molecules are expected to pass though the bilayer membrane. The self-assembled vesicles were shown to be catalytically active in a Miyaura-Michael reaction in water.

Molecular dynamics study of lipid bilayers modeling the plasma membranes of mouse hepatocytes and hepatomas.

Molecular dynamics (MD) calculations of lipid bilayers modeling the plasma membranes of normal mouse hepatocytes and hepatomas in water have been performed under physiological isothermal-isobaric conditions (310.15 K and 1 atm). The changes in the membrane properties induced by hepatic canceration were investigated and were compared with previous MD calculations included in our previous study of the changes in membrane properties induced by murine thymic canceration. The calculated model membranes for normal hepatocytes and hepatomas comprised 23 and 24 kinds of lipids, respectively. These included phosphatidylcholine, phosphatidylethanolamine, phosphatidylserine, phosphatidylinositol, sphingomyelin, lysophospholipids, and cholesterol. We referred to previously published experimental values for the mole fraction of the lipids adopted in the present calculations. The calculated structural and dynamic properties of the membranes such as lateral structure, order parameters, lateral self-diffusion constants, and rotational correlation times all showed that hepatic canceration causes plasma membranes to become more ordered laterally and less fluid. Interestingly, this finding contrasts with the less ordered structure and increased fluidity of plasma membranes induced by thymic canceration observed in our previous MD study.

Dissociation of methane hydrate in aqueous NaCl solutions.

Molecular dynamics simulations of the dissociation of methane hydrate in aqueous NaCl solutions are performed. It is shown that the dissociation of the hydrate is accelerated by the formation of methane bubbles both in NaCl solutions and in pure water. We find two significant effects on the kinetics of the hydrate dissociation by NaCl. One is slowing down in an early stage before bubble formation, and another is swift bubble formation that enhances the dissociation. These effects arise from the low solubility of methane in NaCl solution, which gives rise to a nonuniform spatial distribution of solvated methane in the aqueous phase. We also demonstrate that bubbles form near the hydrate interface in dense NaCl solutions and that the hydrate dissociation proceeds inhomogeneously due to the bubbles.

Effect of bubble formation on the dissociation of methane hydrate in water: a molecular dynamics study.

We investigate the dissociation of methane hydrate in liquid water using molecular dynamics simulations. As dissociation of the hydrate proceeds, methane molecules are released into the aqueous phase and eventually they form bubbles. It is shown that this bubble formation, which causes change in the methane concentration in the aqueous phase, significantly affects the dissociation kinetics of methane hydrate. A large system size employed in this study makes it possible to analyze the effects of the change in the methane concentration and the formation of bubbles on the dissociation kinetics in detail. It is found that the dissociation rate decreases with time until the bubble formation and then it turns to increase. It is also demonstrated that methane hydrate can exist as a metastable superheated solid if there exists no bubble.

Molecular dynamics study of lipid bilayers modeling the plasma membranes of normal murine thymocytes and leukemic GRSL cells.

Molecular dynamics (MD) calculations for the plasma membranes of normal murine thymocytes and thymus-derived leukemic GRSL cells in water have been performed under physiological isothermal-isobaric conditions (310.15K and 1 atm) to investigate changes in membrane properties induced by canceration. The model membranes used in our calculations for normal and leukemic thymocytes comprised 23 and 25 kinds of lipids, respectively, including phosphatidylcholine, phosphatidylethanolamine, phosphatidylserine, phosphatidylinositol, sphingomyelin, lysophospholipids, and cholesterol. The mole fractions of the lipids adopted here were based on previously published experimental values. Our calculations clearly showed that the membrane area was increased in leukemic cells, and that the isothermal area compressibility of the leukemic plasma membranes was double that of normal cells. The calculated membranes of leukemic cells were thus considerably bulkier and softer in the lateral direction compared with those of normal cells. The tilt angle of the cholesterol and the conformation of the phospholipid fatty acid tails both showed a lower level of order in leukemic cell membranes compared with normal cell membranes. The lateral radial distribution function of the lipids also showed a more disordered structure in leukemic cell membranes than in normal cell membranes. These observations all show that, for the present thymocytes, the lateral structure of the membrane is considerably disordered by canceration. Furthermore, the calculated lateral self-diffusion coefficient of the lipid molecules in leukemic cell membranes was almost double that in normal cell membranes. The calculated rotational and wobbling autocorrelation functions also indicated that the molecular motion of the lipids was enhanced in leukemic cell membranes. Thus, here we have demonstrated that the membranes of thymocyte leukemic cells are more disordered and more fluid than normal cell membranes.

MODYLAS: A Highly Parallelized General-Purpose Molecular Dynamics Simulation Program for Large-Scale Systems with Long-Range Forces Calculated by Fast Multipole Method (FMM) and Highly Scalable Fine-Grained New Parallel Processing Algorithms.

Our new molecular dynamics (MD) simulation program, MODYLAS, is a general-purpose program appropriate for very large physical, chemical, and biological systems. It is equipped with most standard MD techniques. Long-range forces are evaluated rigorously by the fast multipole method (FMM) without using the fast Fourier transform (FFT). Several new methods have also been developed for extremely fine-grained parallelism of the MD calculation. The virtually buffering-free methods for communications and arithmetic operations, the minimal communication latency algorithm, and the parallel bucket-relay communication algorithm for the upper-level multipole moments in the FMM realize excellent scalability. The methods for blockwise arithmetic operations avoid data reload, attaining very small cache miss rates. Benchmark tests for MODYLAS using 65 536 nodes of the K-computer showed that the overall calculation time per MD step including communications is as short as about 5 ms for a 10 million-atom system; that is, 35 ns of simulation time can be computed per day. The program enables investigations of large-scale real systems such as viruses, liposomes, assemblies of proteins and micelles, and polymers.

A molecular dynamics study of the lateral free energy profile of a pair of cholesterol molecules as a function of their distance in phospholipid bilayers.

Free energy profile of a pair of cholesterol molecules in a leaflet of 1-palmitoyl-2-oleoyl-phosphatidylcholine (POPC) bilayers in the liquid-crystalline phase has been calculated as a function of their lateral distance using a combination of NPT-constant atomistic molecular dynamics calculations (P = 1 atm and T = 310.15 K) and the thermodynamic integration method. The calculated free energy clearly shows that the two cholesterol molecules form a dimer separated by a distance of 1.0-1.5 nm in POPC bilayers. Well depth of the free energy profile is about 3.5 kJ/mol, which is comparable to the thermal energy k(B)T at 310.15 K. This indicates that the aggregation of cholesterol molecules in the bilayers depends on the temperature as well as the concentration of the system. The free energy function obtained here may be used as a reference when coarse grained potential model is investigated for this two-component system. Local structure of POPC molecules around two cholesterol molecules has also been investigated.

All-atom molecular dynamics study of a spherical micelle composed of N-acetylated poly(ethylene glycol)-poly(gamma-benzyl L-glutamate) block copolymers: a potential carrier of drug delivery systems for cancer.

An all-atom molecular dynamics simulation of a spherical micelle composed of amphiphilic N-acetylated poly(ethylene glycol)-poly(gamma-benzyl L-glutamate) (PEG-PBLG-Ac) block copolymers was performed in aqueous solution at 298.15 K and 1 atm. Such copolymers have received considerable attention as carriers in drug delivery systems. In this study, we used copolymers consisting of 11 EG units and 9 BLG units as models. Starting from the copolymers arranged spherically, the calculation predicted an equilibrium state consisting of a slightly elliptical micelle structure with a hydrophobic PBLG inner core and a hydrophilic PEG outer shell. The micelle structure was dynamically stable during the simulation, with the PEG blocks showing a compact helical conformation and the PBLG blocks an alpha-helix form. Multiple hydrogen bonds with solvent water molecules stabilized the helical conformation of the PEG blocks, leading to their hydration as shown by longer residence times of water molecules near the PEG ether oxygen atoms compared with that of bulk water. Some water molecules have also been found distributed within the hydrophobic core; they showed continuous exchange with bulk water during the simulation. Those molecules existed mostly as a cluster in spaces between the copolymers, forming hydrogen bonds among themselves as well as with the hydrophobic core through hydrophilic groups such as esters and amides. The water molecules forming hydrogen bonds with the micelle may play an important role in the stabilization of the micelle structure.

Hydrogen-bonded clusters on the vapor/ethanol-aqueous-solution interface.

The structure of the vapor/ethanol-aqueous-solution interface has been carefully investigated focusing on an intermolecular hydrogen bond (HB) and molecular clusters bound by HBs. This paper is a continuation of our previous molecular dynamics (MD) study (Langmuir 2005, 21, 10885), and all analysis was performed based on five independent adsorption-equilibrated configurations of a slab of ethanol solution at 298.15 K, where the ethanol mole fraction of the solution, chi(e), is 0.0052, 0.012, 0.024, 0.057, and 0.12, respectively. The geometrical definition of HB enabled the detection of the HB between ethanol-ethanol, ethanol-water, and water-water molecules. The variations of the density of HB and the coordination number of HB across the vapor/solution interface were analyzed. Analysis on the density of HB reveals that a monolayer of adsorbed ethanol can be classified into two parts where ethanol molecules prefer to form HBs with each other and ethanol molecules prefer to form HBs with water molecules. Despite chi(e), the coordination number of ethanol-ethanol HB monotonically increases toward the vapor region, while those of ethanol-water and water-water HBs monotonically decrease. In addition, the variation of the mean size of both ethanol one-component clusters and ethanol/water binary clusters across the interface were analyzed. The mean size of an ethanol one-component cluster and that of an ethanol/water binary cluster are expressed as a maximum at the interface. These behaviors are linked with the size distributions of both one-component and binary clusters. A relatively large system in this calculation also enables detailed discussion about the molar dependency of the bulk structural properties of an ethanol solution.

Two-dimensional supercritical behavior of an ethanol monolayer: a molecular dynamics study.

The two-dimensional (2D) supercritical behavior of an ethanol monolayer formed at the vapor/liquid interface of an ethanol solution has been investigated by a molecular dynamics (MD) calculations with a combination of the OPLS-UA and SPC/E potential models. A 100 A thick slab of ethanol solution was placed at the volume center of the rectangular unit cell by 10 A thick nonabsorbate water surfaces. With such an initial configuration, five independent 15 ns NVT constant MD calculations were carried out under 298.15 K, in which the initial ethanol mole fraction of the bulk solution layer was set to 0.010, 0.022, 0.045, 0.10, and 0.20, respectively. The 2D radial distribution function (rdf) of an adsorbed ethanol molecule showed that the ethanol monolayer could be regarded as a 2D fluid where the adsorbed ethanol molecule had an effective 2D diameter of 4.65 A. On the basis of the 2D rdf result, 2D cluster analysis was carried out from the perspective of the percolation theory. It is confirmed that the critical area occupation probability density, the critical exponents, and the fractal dimension of both nonpercolating and percolating clusters satisfied their nature of universality. Therefore, we concluded that an ethanol monolayer formed at the vapor/liquid interface of ethanol solution behaves as a 2D supercritical fluid at 298.15 K.