Full Molecular Dynamics Simulations Research Articles

Multiscale simulation of nanodrop over surfaces with varying hydrophilicity

Macroscale solvers primarily depend on empirical constitutive relations, which curbs their applicability and reliability. On the other hand, though the microscale solver is much more accurate, it cannot be used for engineering scale domains because of the high computational cost. As an alternative, in the present study, a new adiabatic two-phase multiscale solver has been developed by coupling Finite volume and molecular dynamics solvers with the help of the domain decomposition method. Due to the finiteness of the molecular subdomain, molecules near the boundary miss forces from the missing molecules. An effort has been made to develop a new boundary force model which mimics the effect of those missing molecules. The spreading of a drop has been simulated over surfaces with varying hydrophilicity using the solver and compared with the result of full molecular dynamics simulation. Further, the computational time has been compared with the full molecular dynamics simulation for different drop sizes to show the usefulness of the developed model. Drop translation over biphilic surfaces has also been simulated with the new solver.

Intelligent dissipative particle dynamics: Bridging mesoscopic models from microscopic simulations via deep neural networks

We develop a bottom-up coarse graining framework to construct mesoscopic force fields directly from microscopic dynamics by a machine learning technique. This scheme, called the intelligent dissipative particle dynamics (IDPD), holds various consistencies, including configuration, phase, static and dynamic consistencies that are incorporated into a deep neural network. The network is trained with full atomic data in a special way for preserving the natural symmetries of the system, and generates both the conservative and non-conservative force fields in the mesoscopic system. Subsequently, the force fields are used in the dissipative particle dynamics (DPD) that is considered as the effective mesoscopic dynamics resulting from the Mori–Zwanzig projection of the underlying atomic dynamics. To show the validity and applicability of IDPD, we coarse grain star polymer and methane fluids respectively, from full molecular dynamics (MD) simulations. Quantitative comparisons between the MD and IDPD indicate that the IDPD is able to reproduce both the static and dynamic properties of underlying MD system, with the aids of a pressure correction and a diffusion rescaling, even when coarse graining multiple molecules into a particle. It follows that the IDPD is competent for preserving the dynamic properties and coarse graining non-bonded molecules, which are common challenges for most of coarse graining methods.

First-Principles Modeling of Protein/Surface Interactions. Polyglycine Secondary Structure Adsorption on the TiO2 (101) Anatase Surface Adopting a Full Periodic Approach.

Computational modeling of protein/surface systems is challenging since the conformational variations of the protein and its interactions with the surface need to be considered at once. Adoption of first-principles methods to this purpose is overwhelming and computationally extremely expensive so that, in many cases, dramatically simplified systems (e.g., small peptides or amino acids) are used at the expenses of modeling nonrealistic systems. In this work, we propose a cost-effective strategy for the modeling of peptide/surface interactions at a full quantum mechanical level, taking the adsorption of polyglycine on the TiO2 (101) anatase surface as a test case. Our approach is based on applying the periodic boundary conditions for both the surface model and the polyglycine peptide, giving rise to full periodic polyglycine/TiO2 surface systems. By proceeding this way, the considered complexes are modeled with a drastically reduced number of atoms compared with the finite-analogous systems, modeling the polypeptide structures at the same time in a realistic way. Within our modeling approach, full periodic density functional theory calculations (including implicit solvation effects) and ab initio molecular dynamics (AIMD) simulations at the PBE-D2* theory level have been carried out to investigate the adsorption and relative stability of the different polyglycine structures (i.e., extended primary, β-sheet, and α-helix) on the TiO2 surface. It has been found that, upon adsorption, secondary structures become partially denatured because the peptide C═O groups form Ti-O═C dative bonds. AIMD simulations have been fundamental to identify these phenomena because thermal and entropic effects are of paramount importance. Irrespective of the simulated environments (gas phase and implicit solvent), adsorption of the α-helix is more favorable than that of the β-sheet because in the former, more Ti-O═C bonds are formed and the adsorbed secondary structure results less distorted with respect to the isolated state. Under the implicit water solvent, additionally, adsorbed β-sheet structures weaken with respect to their isolated states as the H-bonds between the strands are longer due to solvation effects. Accordingly, the results indicate that the preferred conformation upon adsorption is the α-helix over the β-sheet.

Toxicity inhibition strategy of microplastics to aquatic organisms through molecular docking, molecular dynamics simulation and molecular modification

In the present study, the combined toxic effect of microplastics and their additives (five) on aquatic organisms (zebrafish) was studied using full factorial design method, molecular docking, and molecular dynamics (MD) simulation technology. The aquatic toxicity control programmer was designed to improve the optimal combination of plasticizer and microplastics based on the design of environment-friendly phthalic acid ester (PAE) derivatives. First, a total of 64 groups of microplastic-additives were designed using the full factorial design method. Next, the microplastic-additives and aquatic receptor protein were docked together, and the binding energy of these complexes was calculated using the MD simulation method. The results revealed that the aquatic toxicity effects of different microplastic-additive combinations were variable; therefore, the optimal combination of microplastics exhibiting the lowest aquatic toxicity effect could be screened out. Base on the analyzing the bonding effect and surrounded amino acid residues between the microplastic additives and receptor protein, the main driving forces for the binding of the microplastic-additive and the protein were hydrophobic force, hydrogen bonding force and electrostatic force. The main effects and the second-order interaction of the microplastic-additives combination were analyzed using the fixed-effect model. The main additives that affect the aquatic toxicity of the microplastics can be known. In addition, based on the MD simulation of the molecular replacement of PAE derivatives, the optimal level of component combination of low aquatic toxicity effect of microplastics was constructed.

A multiscale volume of fluid method with self-consistent boundary conditions derived from molecular dynamics

Molecular dynamics (MD) and volume of fluid (VOF) are powerful methods for the simulation of dynamic wetting at the nanoscale and macroscale, respectively, but the massive computational cost of MD and the sensitivity and uncertainty of boundary conditions in VOF limit their applications to other scales. In this work, we propose a multiscale simulation strategy by enhancing VOF simulations using self-consistent boundary conditions derived from MD. Specifically, the boundary conditions include a particular slip model based on the molecular kinetic theory for the three-phase contact line to account for the interfacial molecular physics, the classical Navier slip model for the remaining part of the liquid–solid interface, and a new source term supplemented to the momentum equation in VOF to replace the convectional dynamic contact angle model. Each slip model has been calibrated by the MD simulations. The simulation results demonstrate that with these new boundary conditions, the enhanced VOF simulations can provide consistent predictions with full MD simulations for the dynamic wetting of nanodroplets on both smooth and pillared surfaces, and its performance is better than those with other VOF models, especially for the pinning–depinning phenomenon. This multiscale simulation strategy is also proved to be capable of simulating dynamic wetting above the nanoscale, where the pure MD simulations are inaccessible due to the computational cost.

Full molecular dynamics simulations of molecular liquids for single-beam spectrally controlled two-dimensional Raman spectroscopy.

Single-beam spectrally controlled (SBSC) two-dimensional (2D) Raman spectroscopy is a unique 2D vibrational measurement technique utilizing trains of short pulses that are generated from a single broadband pulse by pulse shaping. This approach overcomes the difficulty of 2D Raman spectroscopy in dealing with small-signal extraction and avoids complicated low-order cascading effects, thus providing a new possibility for measuring the intramolecular and intermolecular modes of molecular liquids using fifth-order 2D Raman spectroscopy. Recently, for quantitatively investigating the mode-mode coupling mechanism, Hurwitz et al. [Opt. Express 28, 3803 (2020)] have developed a new pulse design for this measurement to separate the contributions of the fifth- and third-order polarizations, which are often overlapped in the original single-beam measurements. Here, we describe a method for simulating these original measurements and the new 2D Raman measurements on the basis of a second-order response function approach. We carry out full molecular dynamics simulations for carbon tetrachloride and liquid water using an equilibrium-nonequilibrium hybrid algorithm, with the aim of explaining the key features of the SBSC 2D Raman spectroscopic method from a theoretical point of view. The predicted signal profiles and intensities provide valuable information that can be applied to 2D spectroscopy experiments, allowing them to be carried out more efficiently.

Methanol carbonylation over acid mordenite: Insights from ab initio molecular dynamics and machine learning thermodynamic perturbation theory

In this work we present a detailed ab initio study of the carbonylation reaction of methoxy groups in the zeolite mordenite, as it is the rate determining step in a series of elementary reactions leading to ethanol. For the first time we employ full molecular dynamics simulations to evaluate free energies of activation for the reactions in side pockets and main channels. Results show that the reaction in the side pocket is preferred and, when dispersion interactions are taken into account, this preference becomes even stronger. This conclusion is confirmed using multiple levels of density functional theory approximations with (PBE-D2, PBE-MBD, and vdW-DF2-B86R) or without (PBE, HSE06) dispersion corrections. These calculations, that in principle would require several demanding molecular dynamics simulations, were made possible at a minimal computational cost by using a newly developed approach that combines thermodynamic perturbation theory with machine learning.

Lattice model for active flows in microchannels.

We introduce a one-dimensional lattice model to study active particles in narrow channel connecting finite reservoirs. The model describes interacting run-and-tumble swimmers exerting pushing forces on neighboring particles, allowing the formation of long active clusters inside the channel. Our model is able to reproduce the emerging oscillatory dynamics observed in full molecular dynamics simulations of self-propelled bacteria [Paoluzzi et al., Phys. Rev. Lett. 115, 188303 (2015)PRLTAO0031-900710.1103/PhysRevLett.115.188303] and allows us to extend in a simple way the analysis to a wide range of system parameters (box length, number of swimmers), taking into account different physical conditions (presence or absence of tumbling, different forms of the entrance probability into the channel). We find that the oscillatory behavior is suppressed for short channels length L<L^{*} and for high tumbling rates λ>λ^{*}, with threshold values L^{*} and λ^{*} which in general depend on physical parameters. Moreover, we find that oscillations persist by using different entrance probabilities, which, however, affect the oscillation properties and the filling dynamics of reservoirs.

Concurrent Multiscale Simulations of Rough Lubricated Contact of Aluminum Single Crystal

In this paper, a concurrent multiscale simulation strategy coupling atomistic and continuum models was proposed to investigate the three-dimensional contact responses of aluminum single crystal under both dry and lubricated conditions. The Hertz contact is performed by using both the multiscale and full molecular dynamics (MD) simulations for validation. From the contact area, kinetic energy and stress continuity aspects, the multiscale model shows good accuracy. It can also save at least five times the computational time compared with the full MD simulations for the same domain size. Furthermore, the results of lubricated contact show that the lubricant molecules could effectively cover the contact surfaces; thereby separating the aluminum surfaces and bearing the support loads. Moreover, the surface topography could be protected by the thin film formed by the lubricant molecules. It has been found that the contact area decreases obviously with increasing the magnitude of load under both dry and lubricated contacts. Besides, a decrease in contact area is also seen when the number of lubricant molecules increases. The present study has confirmed that the dimension of lubricated contacts could be greatly expanded during the simulation using the proposed multiscale method without sacrificing too much computational time and accuracy.

Evaluating Self-Assembly Propensity of Tetra-Peptide using MD and Machine Learning

Evaluation of the self-assembly propensity such as aggregation and amyloid formation of a peptide is important to elucidate protein property. It is also useful in protein engineering for avoiding the aggregation of designed proteins. And, self-assembly structures with high stability are used for nanomaterials. As measuring aggregation propensity of many peptides using experimental methods are difficult, molecular dynamics (MD) simulation has been used for predicting self-assembling peptides. It has been reported that di-peptides (202=400 types) and tri-peptides (203=8000 types) aggregation propensity can be evaluated using MD simulation. However, such methods are difficult to be applied to tetra-peptides due to a large sequence space (204=160,000 types). In this study, we combined the MD simulation and machine learning to accelerate the evaluation of aggregation propensity. We conducted full atomistic MD simulation for 319 peptides and used them for training a machine-learning model to search all possible tetra-peptides. The two scores were used as the measures for self-assembly propensity: AP score related to the accessible surface area of peptides and Beta score related to the possibility of β-sheets structure. We validated that machine learning can accurately predict the result of MD simulation with Pearson's correlation coefficients of 0.88 (AP score) and 0.86 (Beta score). Top 100 predictions suggest that the histidine at C-terminus promotes the formation of a β-sheets structure. We experimentally confirmed that a selected peptide with N-terminal proline predicted to have high AP score and Beta score was aggregated and then changed to fibrous crystals. The crystal structure determined using X-ray diffraction data showed β-sheets structure. These results demonstrate that our framework combining machine learning and MD simulation is a promising strategy for evaluating a property of peptides and proteins.

Numerical study on the immiscible two-phase flow in a nano-channel using a molecular-continuum hybrid method

The present study investigated the two-phase immiscible Couette flow in which a lower smooth or rough plate is stationary and an upper plate is moving, using the hybrid method which combines the simulation using the lattice Boltzmann method in the bulk region with the molecular dynamics simulation near the lower smooth or rough surface. The computational time required using the hybrid method was compared with that using the full molecular dynamics simulation, showing big savings in the computational time required by using the hybrid method. The boundary conditions formed on the lower surface depend on the surface energy, viscosity ratio and the geometry of the lower plate, resulting in the slip, no-slip and locking boundary conditions. The results corresponding to the case of the rough surface were compared with those for the smooth surface in order to investigate the effect of the roughness on the flow behavior on the lower plate.

Anisotropic Thermal Boundary Resistance across 2D Black Phosphorus: Experiment and Atomistic Modeling of Interfacial Energy Transport.

Interfacial thermal boundary resistance (TBR) plays a critical role in near-junction thermal management of modern electronics. In particular, TBR can dominate heat dissipation and has become increasingly important due to the continuous emergence of novel nanomaterials with promising electronic and thermal applications. A highly anisotropic TBR across a prototype 2D material, i.e., black phosphorus, is reported through a crystal-orientation-dependent interfacial transport study. The measurements show that the metal-semiconductor TBR of the cross-plane interfaces is 241% and 327% as high as that of the armchair and zigzag direction-oriented interfaces, respectively. Atomistic ab initio calculations are conducted to analyze the anisotropic and temperature-dependent TBR using density functional theory (DFT)-derived full phonon dispersion relation and molecular dynamics simulation. The measurement and modeling work reveals that such a highly anisotropic TBR can be attributed to the intrinsic band structure and phonon spectral transmission. Furthermore, it is shown that phonon hopping between different branches is important to modulate the interfacial transport process but with directional preferences. A critical fundamental understanding of interfacial thermal transport and TBR-structure relationships is provided, which may open up new opportunities in developing advanced thermal management technology through the rational control over nanostructures and interfaces.

Accelerating multiscale modelling of fluids with on-the-fly Gaussian process regression

We present a scheme for accelerating hybrid continuum-atomistic models in multiscale fluidic systems by using Gaussian process regression as a surrogate model for computationally expensive molecular dynamics simulations. Using Gaussian process regression, we are able to accurately predict atomic-scale information purely by consideration of the macroscopic continuum-model inputs and outputs and judge on the fly whether the uncertainty of our prediction is at an acceptable level, else a new molecular simulation is performed to continually augment the database, which is never required to be complete. This provides a substantial improvement over the current generation of hybrid methods, which often require many similar atomistic simulations to be performed, discarding information after it is used once. We apply our hybrid scheme to nano-confined unsteady flow through a high-aspect-ratio converging–diverging channel, and make comparisons between the new scheme and full molecular dynamics simulations for a range of uncertainty thresholds and initial databases. For low thresholds, our hybrid solution is highly accurate—around that of thermal noise. As the uncertainty threshold is raised, the accuracy of our scheme decreases and the computational speed-up increases (relative to a full molecular simulation), enabling the compromise between accuracy and efficiency to be tuned. The speed-up of our hybrid solution ranges from an order of magnitude, with no initial database, to cases where an extensive initial database ensures no new MD simulations are required.

Entropy production and Vlasov equation for self-gravitating systems

The evolution of a self-gravitating system to a non-equilibrium steady state occurs through a process of violent relaxation. In the thermodynamic limit the dynamics of a many body system should be governed by the Vlasov equation. Recently, however, a question was raised regarding the validity of Vlasov equation during the process or violent relaxation. In this paper we will explore the entropy production during the relaxation process using N-body molecular dynamics simulations. We will show that the entropy production time grows as , with and in the limit , entropy will remain constant, consistent with the Vlasov equation. Furthermore, we will show that the mean field dynamics constructed on the basis of the Vlasov equation is in excellent agreement with the full molecular dynamics simulations, justifying the applicability of Vlasov equation during the violent relaxation phase of evolution.

A Local Region Molecular Dynamics Simulation Method for Nanoscale Sliding Contacts

Computational efficiency and accuracy always conflict with each other in molecular dynamics (MD) simulations. How to enhance the computational efficiency and keep accuracy at the same time is concerned by each corresponding researcher. However, most of the current studies focus on MD algorithms, and if the scale of MD model could be reduced, the algorithms would be more meaningful. A local region molecular dynamics (LRMD) simulation method which can meet these two factors concurrently in nanoscale sliding contacts is developed in this paper. Full MD simulation is used to simulate indentation process before sliding. A criterion called contribution of displacement is presented, which is used to determine the effective local region in the MD model after indentation. By using the local region, nanoscale sliding contact between a rigid cylindrical tip and an elastic substrate is investigated. Two two-dimensional MD models are presented, and the friction forces from LRMD simulations agree well with that from full MD simulations, which testifies the effectiveness of the LRMD simulation method for two-dimensional cases. A three-dimensional MD model for sliding contacts is developed then to show the validity of the LRMD simulation method further. Finally, a discussion is carried out by the principles of tribology. In the discussion, two two-dimensional full MD models are used to simulate the nanoscale sliding contact problems. The results indicate that original smaller model will induce higher equivalent scratching depth, and then results in higher friction forces, which will help to explain the mechanism how the LRMD simulation method works. This method can be used to reduce the scale of MD model in large scale simulations, and it will enhance the computational efficiency without losing accuracy during the simulation of nanoscale sliding contacts.

A Study on Slip Characteristics Using Hybrid Particle-Continuum Method

An effective boundary potential has been proposed to solve nonperiodic boundary condition (NPBC) of hybrid method. The optimized hybrid method is applied to investigate the influences of the channel height and solid–liquid interaction parameters on slip characteristics of Couette flows in micro/nanochannels. By changing the channel height, we find that the relative slip lengths show the obvious negative correlation with the channel height and fewer density oscillations are generated near the solid wall in the larger channel height. Moreover, we continue to investigate the solid–liquid interaction parameters, including the solid–liquid energy scales ratio (C1) and solid–liquid length scales ratio (C2). The results show that the solid–liquid surface changes from hydrophobic to hydrophilic with the increase of C1, the arrangement of liquid particles adjacent to the solid particles is more disorganized over the hydrophobic solid–liquid surface compared with the hydrophilic surface, and the probability of the liquid particles that appear near the solid particles becomes smaller. Meanwhile, the relative slip lengths are minimum when the liquid and solid particles have the same diameter. Furthermore, the relative slip lengths follow a linear relationship with the shear rate when the solid–liquid interaction parameters change. The plenty computational time has been saved by the present hybrid method compared with the full molecular dynamics simulation (FMD) in this paper.

A scale-up nanoporous membrane centrifuge for reverse osmosis desalination without fouling

A scale-up nanoporous membrane centrifuge is designed and modeled. It can be used for nanoscale scale separation including reverse osmosis desalination. There are micron-size pores on the wall of the centrifuge and nanoscale pores on local graphene membrane patches that cover the micron-size pores. In this work, we derived the critical angular velocity required to counter-balance osmosis force, so that the reverse-osmosis (RO) desalination process can proceed. To validate this result, we conducted a large scale (four million atoms) full atom molecular dynamics (MD) simulation to examine the critical angular velocity required for reverse osmosis at nanoscale. It is shown that the analytical results derived based on fluid mechanics and the simulation results observed in MD simulation are consistent and well matched. The main advantage of such nanomaterial based centrifuge is its intrinsic anti-fouling ability to clear [Formula: see text] and [Formula: see text] ions accumulated at the vicinity of the pores due to the Coriolis effect. Analyses have been conducted to study the relation between osmotic pressure, centrifugal pressure, and water permeability.

Multiscale simulation of dynamic wetting

A sequential multiscale strategy that combines molecular dynamics (MD) with volume of fluid (VOF) simulations is proposed to study the spreading of droplets on surfaces. In this hybrid MD/VOF approach, VOF is applied everywhere in the domain with MD pre-simulations distributed along the wetted interface providing the crucial boundary information for the three-phase contact-line dynamics and the solid/liquid interfacial slip. For the latter, molecular shear flow simulations of the liquid in contact with the substrate are used to measure the local slip length to calibrate the Navier slip model. For the contact-line we use MD simulations of nanodroplets spreading on the substrate to calibrate the molecular kinetic theory (MKT) model for both the dynamic contact angle and the slip velocity. We validate this multiscale model for water nanodroplets spreading over a platinum surface by comparing our sequential hybrid simulations with the equivalent dynamics in a full MD treatment. We demonstrate that for nanodroplets spreading on surfaces, applying a dynamic contact angle model is not sufficient to pick up the molecular effects; we need to account for slip across the entire solid/liquid interface, in particular the large slip behaviour at the contact-line. We also demonstrate the application of this multiscale method to larger nanodroplets (up to ∼100nm diameters), where full MD simulation would be computationally intractable. We find that as the droplet size increases, the slip in the contact line region becomes less important. To simulate the full range of nano to macro droplets, an improved way of dealing with the VOF method is needed to reduce the overall number of grid points.

Hydrogen bonding and vibrational energy relaxation of interfacial water: A full DFT molecular dynamics simulation.

The air-water interface has been a subject of extensive theoretical and experimental studies due to its ubiquity in nature and its importance as a model system for aqueous hydrophobic interfaces. We report on the structure and vibrational energy transfer dynamics of this interfacial water system studied with equilibrium and non-equilibrium molecular dynamics simulations employing a density functional theory -based description of the system and the kinetic energy spectral density analysis. The interfacial water molecules are found to make fewer and weaker hydrogen (H)-bonds on average compared to those in the bulk. We also find that (i) the H-bonded OH groups conjugate to the free OH exhibit rather low vibrational frequencies (3000-3500 cm-1); (ii) the presence of a significant fraction (>10%) of free and randomly oriented water molecules at the interface ("labile water"), neither of whose OH groups are strong H-bond donors; (iii) the inertial rotation of free OH groups, especially from the labile water, contribute to the population decay of excited free OH groups with comparable rate and magnitude as intramolecular energy transfer between the OH groups. These results suggest that the labile water, which might not be easily detectable by the conventional vibrational sum frequency generation method, plays an important role in the surface water dynamics.

Effects of density on flow in a nano channel using a molecular-continuum hybrid method

A molecular-continuum hybrid method was developed to simulate micro- and nano-scale fluid flows that cannot be predicted using continuum fluidics. Molecular dynamics simulation was used near stationary solid surfaces, and Navier-Stokes equations were used in other regions. We carried out Couette flow simulation using this hybrid method and validated the results by comparing them with the analytical solution. We also studied the dependence of the velocity slip and slip length on the surface energy, liquid density, and roughness for a liquid channel flow with and without nano-structures on the solid surface. The behavior of the liquid near the solid wall changed with the surface energy as well as the liquid density. The variation of the velocity slip and slip length according to the surface energy also depended on the liquid density as well as the surface roughness. We compared the required computational time obtained from the molecular-continuum hybrid method with that obtained from full molecular dynamics simulation under the same computational condition, giving much shorter computational time for the case using the molecular-continuum hybrid method than that for full molecular dynamics simulation.