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.

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.

Exact long-range Coulombic energy calculation for net charged systems neutralized by uniformly distributed background charge using fast multipole method and its application to efficient free energy calculation.

In molecular dynamics (MD) calculations of the free energies of ions and ionic molecules, we often encounter net charged molecular systems where the electrical neutrality condition is broken. This charge causes a problem in the evaluation of long-range Coulombic interactions under periodic boundary conditions. A standard remedy for this problem is to consider a hypothetical homogeneous background charge density to neutralize the total system. Here, we present a new expression for the evaluation of Coulombic interactions for such systems including background charge using the fast multipole method (FMM). Furthermore, an efficient scheme is developed to evaluate solute-solvent interaction energies using the FMM, reducing the computational burden for the far-field part. We calculate the hydration free energies of Mg2+, Na+, and Cl- ions dissolved in a neutral solvent using the new expression. The calculated free energies show good agreement with the results obtained using the well-established particle mesh Ewald method. This demonstrates the validity of the proposed expression. This work should make a contribution to highly parallelized MD calculations for large-scale charged systems (particularly, those with over million particles).

Molecular dynamics study of solubilization of cyclohexane, benzene, and phenol into mixed micelles composed of sodium dodecyl sulfate and octaethylene glycol monododecyl ether.

Molecular dynamics calculations of a mixed micelle composed of sodium dodecyl sulfate (SDS) and octaethylene glycol monododecyl ether (C12 E8 ) were performed for six compositions (SDS/C12 E8 = 100/0, 80/20, 60/40, 40/60, 20/80, and 0/100) to investigate the composition dependence of the mixed micelle structure and solubilization of cyclohexane, benzene, and phenol molecules by the micelle. The radial density distribution of the hydrophilic polyoxyethylene (POE) group of C12 E8 as a function of distance from the micelle center is very sharp for micelles with high SDS content because the POE group captures a Na+ ion in solution and wraps around it to form a compact crown-ether-like complex. The hydrophobic dodecyl groups of SDS and C12 E8 were separately distributed in the mixed micelle core. ΔG(r) evaluated for each solute showed that despite the structural changes of the micelle the binding strength of the solute molecules to the micelle did not change significantly. © 2019 Wiley Periodicals, Inc.

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.

Spherical harmonics analysis of surface density fluctuations of spherical ionic SDS and nonionic C12E8 micelles: A molecular dynamics study.

The surface structure and its fluctuation of spherical micelles were investigated using a series of density correlation functions newly defined by spherical harmonics and Legendre polynomials based on the molecular dynamics calculations. To investigate the influence of head-group charges on the micelle surface structure, ionic sodium dodecyl sulfate and nonionic octaethyleneglycol monododecylether (C12E8) micelles were investigated as model systems. Large-scale density fluctuations were observed for both micelles in the calculated surface static structure factor. The area compressibility of the micelle surface evaluated by the surface static structure factor was tens-of-times larger than a typical value of a lipid membrane surface. The structural relaxation time, which was evaluated from the surface intermediate scattering function, indicates that the relaxation mechanism of the long-range surface structure can be well described by the hydrostatic approximation. The density fluctuation on the two-dimensional micelle surface has similar characteristics to that of three-dimensional fluids near the critical point.

Molecular dynamics study of the potential of mean force of SDS aggregates.

In our previous study, all-atomistic molecular dynamics (MD) calculations have been carried out for the aggregation of ionic sodium dodecyl sulfate in water [S. Kawada et al., Chem. Phys. Lett. 646, 36 (2016)]. Aggregates of 20-30 dodecyl sulfate ions were formed within a short MD run for 10 ns. However, further aggregation did not occur despite a long MD calculation for more than 100 ns. This suggests that strong electrostatic repulsive interactions between the aggregates prevent the fusion of the aggregates. In the present study, mean force and potential of mean force acting between two aggregates with aggregation number N = 30 were evaluated as a function of their separation by MD calculations. The repulsive force becomes strong with decreasing distance between the two aggregates before they merge into one. An origin of the repulsive force is an electric double layer formed by the sulfate group and counter sodium ions. Strength of the repulsive force is in good agreement with the theoretical value given by the Derjaguin-Landau-Verwey-Overbeek (DLVO) theory. Once the aggregates establish contact, the force between them turns to be a large attractive force that can be explained by the interfacial tension. In order to form a single micelle from the two aggregates, it is necessary for them to climb over a free energy barrier of 23 kJ/mol. Once, the barrier is overcome, the micelle is stabilized by ∼200 kJ/mol. The time constant of aggregation evaluated from the calculated free energy barrier was about 28 µs at the concentration in our previous study.

A molecular dynamics study of local pressures and interfacial tensions of SDS micelles and dodecane droplets in water.

To obtain the radial (normal) and lateral (transverse) components of the local pressure tensor, PN(R) and PT(R), respectively, and the interfacial tension of micelles, molecular dynamics (MD) calculations were performed for spherical sodium dodecyl sulfate (SDS) micelles. The local pressure tensor was calculated as a function of radial distance R using the Irving-Kirkwood formula. Similar MD calculations were also carried out for an n-dodecane droplet in water to compare the differences in the local pressure and interfacial tension values with those of the micelles. The calculated interfacial tensions were 20 ± 5 and 44 ± 10 mN/m for the SDS micelles and dodecane droplets, respectively. The excess free energies due to the interfacial tension were 340 and 1331 kJ/mol for the SDS micelle and dodecane droplet, respectively. The micelles are stabilized by 991 kJ/mol by covering their hydrophobic cores with hydrophilic groups. The dodecane droplet has a large interfacial tension caused by the zero or positive values of PN(R) - PT(R) at all values of R. In contrast, the small interfacial tension in the SDS micelles comes from the negative PN(R) - PT(R) values over a wide range of R. The pressure difference between the inside and outside of the oil droplet and its interfacial tension well satisfies the Laplace equation. However, the hydrophobic core of the SDS micelle is quite different from the liquid alkane, and the SDS micelles do not follow Laplace's picture. Decomposing the interfacial tension into contributions from various interactions, it is found that those between charged and polar groups dominate the interfacial tension of the SDS micelles. The positive electrostatic potential (1.3 V) on the micelle surface and the negative potential (-0.15 V) on the oil droplet contribute to the interfacial tensions by 19 and 0.5 mN/m, respectively. Thus, the interfacial tension of the SDS micelles is produced by electrostatic interactions, in contrast to the dodecane droplet.

A molecular dynamics study of the breathing and deforming modes of the spherical ionic SDS and nonionic C12E8 micelles.

In order to investigate shape of the micelles and its thermal fluctuations, molecular dynamics calculations have been performed for spherical ionic sodium dodecyl sulfate (SDS) and nonionic octaethyleneglycol monododecyl ether (C12E8) micelles. New statistical functions suitable for extracting the fluctuations of the shape of the spherical micelles were defined using spherical harmonics and Legendre polynomials. The breathing and deforming modes of the SDS and C12E8 micelles were analyzed in detail based on these new functions. The elastic nature of the micelle core was also discussed. The present analysis gives a new molecular picture that the micelle shape is a superposition of the various kinds of breathing and deforming modes, and each mode has a specific relaxation time of the shape fluctuation.

Drug binding and mobility relating to the thermal fluctuation in fluid lipid membranes.

Drug binding and mobility in fluid lipid bilayer membranes are quantified in situ by using the multinuclear solution NMR combined with the pulsed-field-gradient technique. One-dimensional and pulsed-field-gradient (19)F and (1)H NMR signals of an anticancer drug, 5-fluorouracil (5FU) are analyzed at 283-313 K in the presence of large unilamellar vesicles (LUVs) of egg phosphatidylcholine (EPC) as model cell membranes. The simultaneous observation of the membrane-bound and free 5FU signals enables to quantify in what amount of 5FU is bound to the membrane and how fast 5FU is moving within the membrane in relation to the thermal fluctuation of the soft, fluid environment. It is shown that the mobility of membrane-bound 5FU is slowed down by almost two orders of magnitude and similar to the lipid movement in the membrane, the movement closely related to the intramembrane fluidity. The mobility of 5FU and EPC is, however, not similar at 313 K; the 5FU movement is enhanced in the membrane as a result of the loose binding of 5FU in the lipid matrices. The membrane-bound fraction of 5FU is approximately 0.1 and almost unaltered over the temperature range examined. It is also independent of the 5FU concentration from 2 to 30 mM with respect to the 40-50 mM LUV. The free energy of the 5FU binding is estimated at -4 to -2 kJ/mol, the magnitude always close to the thermal fluctuation, 2.4-2.6 kJ/mol.

Lateral diffusion of lipids separated from rotational and translational diffusion of a fluid large unilamellar vesicle.

A new method to separate lateral diffusion of lipids in spherical large unilamellar vesicles from the rotational and the translational diffusion of the vesicle as a whole is proposed. The lateral diffusion coefficient DL is obtained as a time-dependent part of the observed diffusion coefficient in vesicles of 800-nm diameters, by systematically changing the diffusion time interval of the high-field-gradient NMR measurement. Although the lipid is in a confined space, the DL of 1,2-dipalmitoyl-sn-glycero-3-phosphocholine is (1.5±0.6)×10(-11) m(2) s(-1) in the fluid state at 45°C, more than one order of magnitude faster than the rotational and the translational diffusion coefficients of the vesicle by the hydrodynamic continuum model. The method provides a potential for quantifying the lateral diffusion of lipids and proteins in fluid bilayer vesicles as model cell membranes in a natural manner.

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.

Binding of hydrophobic fluorinated bisphenol A to large unilamellar vesicles of egg phosphatidylcholine.

The kinetics of membrane binding and dissociation of fluorinated bisphenol A (FBPA, (CF(3))(2)C(C(6)H(4)OH)(2)) is quantified by 1D (19)F NMR spectra in situ. Although the bound and free components are in fast exchange, the rate constants and bound fraction is nonetheless determined from an analysis of the spectra. The analysis relies on the expression of 1D NMR signal intensity by a set of Bloch equations with exchange terms. The time span of the binding and dissociation of hydrophobic FBPA to large unilamellar vesicles of egg phosphatidylcholine (EPC) is 10(-3) to 10(-2) s. The rates of FBPA binding and dissociation are kept constant per EPC molecule even at different concentrations of the vesicle. The free energy of FBPA transfer is -20 ± 2 kJ/mol at 303 K. The process is entropy-driven. The efficiency of FBPA transfer is enhanced by a factor of 7 × 10(4), as compared with the hydrophilic 5-fluorouracil.

Kinetics of binding and diffusivity of leucine-enkephalin in large unilamellar vesicle by pulsed-field-gradient 1H NMR in situ.

The kinetics of binding, the diffusivity, and the binding amount of a neuropeptide, leucine-enkephalin (L-Enk) to lipid bilayer membranes are quantified by pulsed-field-gradient (PFG) 1H NMR in situ. The peptide signal is analyzed by the solution of the Bloch equation with exchange terms in the presence of large unilamellar vesicles (LUVs) as confined, but fluid model cell membranes. Even in the case that the membrane-bound and the free states of L-Enk cannot be distinguished in the one-dimensional NMR spectrum, the PFG technique unveils the bound component of L-Enk after the preferential decay of the free component at the high field gradient. In 100-nm diameter LUVs consisting of egg phosphatidylcholine, the rate constants of the peptide binding and dissociation are 0.040 and 0.40 s-1 at 303 K. This means that the lifetime of the peptide binding is of the order from second to ten-second. The diffusivity of the bound L-Enk is 5×10-12m2/s, almost 60 times as restricted as the movement of free L-Enk at 303K. One-tenth of 5mM L-Enk is bound to 40mM LUV. The binding free energy is calculated to be -2.9 kJ/mol, the magnitude close to the thermal fluctuation, 2.5 kJ/mol. The result demonstrates the potential of PFG 1H NMR to quantify molecular dynamics of the peptide binding to membranes.

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.

A molecular dynamics study of free energy of micelle formation for sodium dodecyl sulfate in water and its size distribution.

Free energy of micelle formation has been evaluated for spherical sodium dodecyl sulfate (SDS) in water by a thermodynamic integration method combined with a series of large-scale molecular dynamics calculations following the chemical species model. In particular, free energy change delta mu(n+1)0 with respect to the addition of one surfactant molecule to the spherical micelle of size n was obtained as a function of n. The free energy profile showed a minimum followed by a maximum, which corresponds to a peak in the size distribution. The calculated peak size n = 57 near its critical micelle concentration is in good agreement with the experimental averaged aggregation number n = 55-75 of the SDS micelle. The distribution showed a rather sharp peak, indicating that the size is almost a monodisperse one. The size is likely to be insensitive to the total concentration of the surfactant.