A density functional study of the photocatalytic degradation of polycaprolactone by the decatungstate anion in acetonitrile solution.

A recent experimental study has reported that decatungstate [W10O32]4- can degrade various polyesters in the presence of light and molecular oxygen [Li et al., Nanoscale, 2023, 15, 15038]. We apply density functional theory to the photocatalyst-polycaprolactone model complex in acetonitrile solution and elucidate the degradation mechanisms and catalytic cycle. We consider hydrogen atom transfer (HAT) and single electron transfer (SET) mechanisms. The potential energy profiles show that the former proceeds exergonically in a single step but that the latter involves a subsequent proton transfer and finally yields HAT products as well. Oxygenated polymer species can regain the transferred hydrogen and regenerate the reduced photocatalyst. We propose a photocatalytic cycle that realizes both the photocatalyst regeneration and the polymer degradation.

A new dynamic Monte Carlo method satisfying n-particle diffusion equation with position-dependent diffusion coefficient, free energy, and intermolecular interactions.

A dynamic Monte Carlo (MC) method recently proposed by us [Nagai et al., J. Chem. Phys. 156, 154506 (2022)] to describe single-particle diffusion of a molecule in a heterogeneous space with position-dependent diffusion coefficient and free energy is generalized here to n-particle dynamics, where n molecules diffuse in heterogeneous media interacting via their intermolecular potential. Starting from the master equation, we give an algebraic proof that the dynamic MC transition probabilities proposed here produce particle trajectories that satisfy the n-particle diffusion equation with position-dependent diffusion coefficient D0i(ri), free energy F1i(ri), and intermolecular interactions Vij(ri, rj). The MC calculations based on this method are compared to molecular dynamics (MD) calculations for two-dimensional heterogeneous Lennard-Jones test systems, showing excellent agreement of the long-distance global diffusion coefficient between the two cases. Thus, the particle trajectories produced by the present MC transition probabilities satisfy the n-particle diffusion equation, and the diffusion equation well describes the long-distance trajectories produced by the MD calculations. The method is also an extension of the conventional equilibrium Metropolis MC calculation for homogeneous systems with a constant diffusion coefficient to the dynamics in heterogeneous systems with a position-dependent diffusion coefficient and potential. In the present method, interactions and dynamics of the real systems are coarse-grained such that the calculation cost is drastically reduced. This provides an approach for the investigation of particle dynamics in very complex and large systems, where the diffusing length is of sub-micrometer order and the diffusion time is of the order of milliseconds or more.

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.

Prediction of self-diffusion coefficients of chemically diverse pure liquids by all-atom molecular dynamics simulations.

Molecular self-diffusion coefficients underlie various kinetic properties of the liquids involved in chemistry, physics, and pharmaceutics. In this study, 547 self-diffusion coefficients are calculated based on all-atom molecular dynamics (MD) simulations of 152 diverse pure liquids at various temperatures employing the OPLS4 force field. The calculated coefficients are compared with experimental data (424 extracted from the literature and 123 newly measured by pulsed-field gradient nuclear magnetic resonance). The calculations well agree with the experimental values. The determination coefficient and root mean square error between the observed and calculated logarithmic self-diffusion coefficients of the 547 entries are 0.931 and 0.213, respectively, demonstrating that the MD calculation can be an excellent industrial tool for predicting, for example, molecular transportation in liquids such as the diffusion of active ingredients in biological and pharmaceutical liquids. The self-diffusion coefficients collected in this study are compiled into a database for broad researches including artificial intelligence calculations.

Global diffusion of hydrogen molecules in the heterogeneous structure of polymer electrolytes for fuel cells: Dynamic Monte Carlo combined with molecular dynamics calculations.

Using our recently developed dynamic Monte Carlo (MC) method [Nagai et al., J. Chem. Phys. 156, 154506 (2022)], we investigated the global diffusion of hydrogen molecules over structural heterogeneities of polymer electrolyte membranes in fuel cells. The three-dimensional position-dependent free energies and the diffusion constants of the hydrogen molecules, required by the present dynamic MC calculations, were taken from our previous study [Nagai et al., J. Chem. Phys. 156, 044507 (2022)] and newly evaluated in this work, respectively. The calculations enabled evaluating the hydrogen dynamics over long-time scales, including global diffusion constants. Based on the calculated global diffusion constants and free energies, the permeability of hydrogen molecules was estimated via the solubility-diffusion model. The estimated values were in good agreement with the reported experimental data, thus validating the present methodology. The analysis of the Monte Carlo trajectories indicated that the main permeation paths are located in the polymer and interfacial phases, although the water phase may make a non-negligible contribution to mass transport.

Dynamic Monte Carlo calculation generating particle trajectories that satisfy the diffusion equation for heterogeneous systems with a position-dependent diffusion coefficient and free energy.

A series of new Monte Carlo (MC) transition probabilities was investigated that could produce molecular trajectories statistically satisfying the diffusion equation with a position-dependent diffusion coefficient and potential energy. The MC trajectories were compared with the numerical solution of the diffusion equation by calculating the time evolution of the probability distribution and the mean first passage time, which exhibited excellent agreement. The method is powerful when investigating, for example, the long-distance and long-time global transportation of a molecule in heterogeneous systems by coarse-graining them into one-particle diffusive molecular motion with a position-dependent diffusion coefficient and free energy. The method can also be applied to many-particle dynamics.

Three-dimensional free-energy landscape of hydrogen and oxygen molecules in polymer electrolyte membranes: Insight into diffusion paths.

Polymer electrolyte membranes, for example, the Nafion™ membranes, used in the fuel cells are responsible for separating reactive gas molecules as well as for the efficient exchange of protons. Although control of the permeation of the gases is important to enhance the fuel cell performance, the mechanism by which hydrogen and oxygen molecules permeate through the membranes remains unclear. To clarify the mechanism, we investigated the three-dimensional free-energy landscape of hydrogen and oxygen molecules in Nafion membranes with various water contents focusing on relevant diffusion paths. Low-free-energy paths are found mainly in the polymer phase and the interfacial region between the polymer and water phases. Thus, the path of the transportation may be attributed to the polymer phase and interfacial phases. However, the free-energy value in the aqueous phase is only slightly higher (∼1-2 kBT) than that in the other two phases, which indicates that a secondary contribution from the aqueous phase is expected. The free-energy landscape in the polymer and interfacial phases was found rugged, while it is comparatively flat in the water phase. We also found that an increase in water content brings about a smoother free-energy landscape in the polymer and interfacial phases. The decreased ruggedness may facilitate the gas diffusivity. These observations help understand the molecular mechanism of the gas diffusion in the membranes.

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.

All-atom molecular dynamics study of hepatitis B virus containing pregenome RNA in solution.

Immature hepatitis B virus (HBV) captures nucleotides in its capsid for reverse transcription. The nucleotides and nucleotide analog drugs, which are triphosphorylated and negatively charged in the cell, approach the capsid via diffusion and are absorbed into it. In this study, we performed a long-time molecular dynamics calculation of the entire HBV capsid containing pregenome RNA to investigate the interactions between the capsid and negatively charged substances. Electric field analysis demonstrated that negatively charged substances can approach the HBV capsid by thermal motion, avoiding spikes. The substances then migrate all over the floor of the HBV capsid. Finally, they find pores through which they can pass through the HBV capsid shell. Free energy profiles were calculated along these pores for small ions to understand their permeability through the pores. Anions (Cl-) showed higher free energy barriers than cations (Na+ and K+) through all pores, and the permeation rate of Cl- was eight times slower than that of K+ or Na+. Furthermore, the ions were more stable in the capsid than in the bulk water. Thus, the HBV capsid exerts ion selectivity for uptake and provides an environment for ions, such as nucleotides and nucleotide analog drugs, to be stabilized within the capsid.

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.

Wettability of a Poly(vinylidene fluoride) Surface by a Pure Good Solvent and a Good Solvent/Nonsolvent Mixture: All-Atom Molecular Dynamics Study.

This study investigated the wettability of poly(vinylidene fluoride) (PVDF) surfaces by a good pure solvent and a good solvent/nonsolvent mixture based on all-atom molecular dynamics (MD) simulations. In particular, droplets of pure N-methyl-2-pyrrolidone (NMP) and of mixed NMP/water molecules were brought into contact with both crystalline and amorphous PVDF surfaces. The contact angles of the macroscopic droplets on the crystalline surface were higher and those on the amorphous surface were lower than the experimental values. As the PVDF sheet surface is a mixture of crystalline and amorphous phases, the experimental contact angles being between those on crystalline and amorphous surfaces is reasonable. On the crystalline surface, the decrease in the contact angle with increasing NMP concentration in the droplets can be explained by the increase in the NMP density near the solid-liquid interface. On the amorphous surface, however, the contact angle is strongly affected by the swelling of PVDF by the mixed droplets at high NMP concentrations. The solvation free energy of PVDF in NMP is greater than that in water, suggesting that this may be a driving force of the swelling of the amorphous PVDF. Furthermore, when the Cassie equation for mixed crystalline and amorphous surfaces was assumed, the calculated contact angle corresponded well with the experimental value.

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.

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 of Polyrotaxane in Solution Investigated by Quasi-Elastic Neutron Scattering and Molecular Dynamics Simulation: Sliding Motion of Rings on Polymer.

In this study, we investigated the molecular dynamics of polyrotaxane (PR), composed of α-cyclodextrins (CDs) and a poly(ethylene glycol) (PEG) axial chain, in solution by means of quasi-elastic neutron scattering (QENS) measurements and full-atomistic molecular dynamics (MD) simulations. From QENS experiments, we estimated the diffusion coefficients of CD and PEG monomers in PR, which are in quantitative agreement with those obtained by MD simulations. By analyzing the simulation results, we succeeded, for the first time, in observing and quantifying the sliding motion of CD along a PEG chain. The diffusion coefficient for the sliding motion is almost 6 times lower than that of the translational diffusion of CD in PR at room temperature. The retardation of the sliding motion is caused by the energy barrier on PEG produced by molecular interactions between CD and PEG. We propose a simple equation to describe the diffusion coefficient of the sliding dynamics in PR by combining the Einstein-Stokes diffusion model and a one-dimensional jump diffusion model. This work provides a general strategy for the molecular designs to control the sliding motion in PR.

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.

Development of dissociative force field for all-atomistic molecular dynamics calculation of fracture of polymers.

A dissociative force field for all-atomistic molecular dynamics calculations has been developed to investigate impact fracture of polymers accompanying dissociation of chemical bonds of polymer main chain. Energy of dimer molecules was evaluated as a function of both bond-length b and bond-angle θ by CASPT2 calculations, whose quality is enough to describe dissociation of chemical bonds. Because we found that the bond dissociation energy D decreases with increasing bond-angle, we employed the Morse-type function VBond (b, θ) = {D - VAngle (θ)}[1 - exp{-α(b - b0 ) - ß(b - b0 )2 }] where a quartic function VAngle (θ) = k1 (θ - θ0 ) + k2 (θ - θ0 )2 + k3 (θ - θ0 )3 + k4 (θ - θ0 )4 . This function reproduced well the CASPT2 potential energy surface in a wide range of b and θ. The parameters have been obtained for four popular glassy polymers, polyethylene, poly(methyl methacrylate), poly(styrene), and polycarbonate. © 2019 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.

Heteroaryldihydropyrimidines Alter Capsid Assembly By Adjusting the Binding Affinity and Pattern of the Hepatitis B Virus Core Protein.

Hepatitis B virus (HBV) infections are a major global health concern, for which heteroaryldihydropyrimidines (HAPs) have been developed. HAPs accelerate and/or result in aberrant capsid assembly; however, their effect on the assembly mechanism is unknown. This study aimed to compare the effects of three representative HAPs on core protein dimer assembly through molecular dynamics simulations and free energy calculations. Molecular docking and equilibrium simulations showed that different HAPs bind at the same binding site and are involved in different interactions. The observed conformational changes in HAPs deter the calculation of binding affinity. Herein, the reduced free energy perturbation/Hamiltonian replica exchange molecular dynamics method was used to enhance sampling during binding affinity calculations, indicating consistency between the binding free energies of HAPs and pEC50. Furthermore, binding pattern analysis revealed that the tetramer could sample flat structures after binding HAPs. The present results suggest a mechanism wherein HAPs accelerate capsid assembly by increasing the binding affinity of dimers, leading to aberrant assembly by altering the binding orientation of dimers.