Efficient lookup table using a linear function of inverse distance squared

The major bottleneck in molecular dynamics (MD) simulations of biomolecules exist in the calculation of pairwise nonbonded interactions like Lennard‐Jones and long‐range electrostatic interactions. Particle‐mesh Ewald (PME) method is able to evaluate long‐range electrostatic interactions accurately and quickly during MD simulation. However, the evaluation of energy and gradient includes time‐consuming inverse square roots and complementary error functions. To avoid such time‐consuming operations while keeping accuracy, we propose a new lookup table for short‐range interaction in PME by defining energy and gradient as a linear function of inverse distance squared. In our lookup table approach, densities of table points are inversely proportional to squared pair distances, enabling accurate evaluation of energy and gradient at small pair distances. Regardless of the inverse operation here, the new lookup table scheme allows fast pairwise nonbonded calculations owing to efficient usage of cache memory. © 2013 Wiley Periodicals, Inc.     

Microsecond simulations of DNA and ion transport in nanopores with novel ion–ion and ion–nucleotides effective potentials

We developed a novel scheme based on the grand‐canonical Monte Carlo/Brownian dynamics simulations and have extended it to studies of ion currents across three nanopores with the potential for single‐stranded DNA (ssDNA) sequencing: solid‐state nanopore Si₃N₄, α‐hemolysin, and E111N/M113Y/K147N mutant. To describe nucleotide‐specific ion dynamics compatible with ssDNA coarse‐grained model, we used the inverse Monte Carlo protocol, which maps the relevant ion–nucleotide distribution functions from all‐atom molecular dynamics (MD) simulations. Combined with the previously developed simulation platform for Brownian dynamics simulations of ion transport, it allows for microsecond‐ and millisecond‐long simulations of ssDNA dynamics in the nanopore with a conductance computation accuracy that equals or exceeds that of all‐atom MD simulations. In spite of the simplifications, the protocol produces results that agree with the results of previous studies on ion conductance across open channels and provide direct correlations with experimentally measured blockade currents and ion conductances that have been estimated from all‐atom MD simulations. © 2014 Wiley Periodicals, Inc.     

Computer simulation of trifluoromethane properties with ab initio force field

Intermolecular interaction potentials of the trifluoromethane dimer in 15 orientations have been calculated using the Hartree‐Fock (HF) self‐consistent theory and the second‐order Møller‐Plesset (MP2) perturbation theory. Single point energies at important geometries were also calibrated by the coupled cluster with single and double and perturbative triple excitation [CCSD(T)] calculations. We have employed Pople's medium size basis sets [up to 6‐311++G(3df,3pd)] and Dunning's correlation consistent basis sets (up to aug‐cc‐pVQZ). Basis set limit potential values were obtained through well‐studied extrapolation methods. The calculated MP2 potential data were employed to parameterize a 5‐site force field for molecular simulations. We performed molecular dynamics simulations using the constructed ab initio force field and compared the simulation results with experiments. Quantitative agreements for the atom‐wise radial distribution functions and the self‐diffusion coefficients over a wide range of experimental conditions can be obtained, thus validating the ab initio force field without using experimental data a priori. © 2011 Wiley Periodicals, Inc. J Comput Chem, 2011     

Force field dependence of phospholipid headgroup and acyl chain properties: Comparative molecular dynamics simulations of DMPC bilayers

The reliability of molecular simulations largely depends on the quality of the empirical force field parameters. Force fields used in lipid simulations continue to be improved to enhance the agreement with experiments for a number of different properties. In this work, we have carried out molecular dynamics simulations of neat DMPC bilayers using united‐atom Berger force field and three versions of all‐atom CHARMM force fields. Three different systems consisting of 48, 72, and 96 lipids were studied. Both particle mesh Ewald (PME) and spherical cut‐off schemes were used to evaluate the long‐range electrostatic interactions. In total, 21 simulations were carried out and analyzed to find out the dependence of lipid properties on force fields, system size, and schemes to calculate long‐range interactions. The acyl chain order parameters calculated from Berger and the recent versions of CHARMM simulations have shown generally good agreement with the experimental results. However, both sets of force fields deviate significantly from the experimentally observed P‐C dipolar coupling values for the carbon atoms that link the choline and glycerol groups with the phosphate groups. Significant differences are also observed in several headgroup parameters between CHARMM and Berger simulations. Our results demonstrate that when changes were introduced to improve CHARMM force field using PME scheme, all the headgroup parameters have not been reoptimized. The headgroup properties are likely to play a significant role in lipid–lipid, protein–lipid, and ligand–lipid interactions and hence headgroup parameters in phospholipids require refinement for both Berger and CHARMM force fields. © 2009 Wiley Periodicals, Inc.J Comput Chem, 2010     

Adsorption of carbon monoxide on Pd low‐index surfaces and (110) reconstructed surface

Adsorption and diffusion of carbon monoxide on Pd low‐index surfaces and missing‐row Pd (110) reconstructed surface have been investigated by the extended London–Eyring–Polyani–Sato (LEPS) method constructed by means of a five‐parameter Morse potential. All critical characteristics, such as adsorption site, adsorption geometry, binding energy, CO vibrational frequency have been obtained and compared with the experimental and theoretical data. On these surfaces, the stable adsorption sites of CO are changed with increasing CO coverage. On the missing‐row Pd (110) reconstructed surface, there are five stable adsorption sites: H₁, H₂ (H₁ and H₂ are threefold hollow sites on (111) subsurface), B (bridge site on the second layer), SB (short‐bridge site), and T (top site). Copyright © 2010 John Wiley & Sons, Ltd.     

Computational design of protein–ligand binding: Modifying the specificity of asparaginyl‐tRNA synthetase

A method for computational design of protein–ligand interactions is implemented and tested on the asparaginyl‐ and aspartyl‐tRNA synthetase enzymes (AsnRS, AspRS). The substrate specificity of these enzymes is crucial for the accurate translation of the genetic code. The method relies on a molecular mechanics energy function and a simple, continuum electrostatic, implicit solvent model. As test calculations, we first compute AspRS‐substrate binding free energy changes due to nine point mutations, for which experimental data are available; we also perform large‐scale redesign of the entire active site of each enzyme (40 amino acids) and compare to experimental sequences. We then apply the method to engineer an increased binding of aspartyl‐adenylate (AspAMP) into AsnRS. Mutants are obtained using several directed evolution protocols, where four or five amino acid positions in the active site are randomized. Promising mutants are subjected to molecular dynamics simulations; Poisson‐Boltzmann calculations provide an estimate of the corresponding, AspAMP, binding free energy changes, relative to the native AsnRS. Several of the mutants are predicted to have an inverted binding specificity, preferring to bind AspAMP rather than the natural substrate, AsnAMP. The computed binding affinities are significantly weaker than the native, AsnRS:AsnAMP affinity, and in most cases, the active site structure is significantly changed, compared to the native complex. This almost certainly precludes catalytic activity. One of the designed sequences has a higher affinity and more native‐like structure and may represent a valid candidate for Asp activity. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2010     

Computer Simulation of Long Side‐Chain Substituted Poly(phenylene vinylene) Polymers

A molecular dynamics (MD) simulation was employed to investigate structure features and segment orientation of four poly(phenylene vinylene) (PPV) derivatives with long, flexible side chains at room temperature. In the simulations, the main chains of the polymers were found to be semirigid and exhibit a tendency to coil into ellipsoidal helices or form zigzag conformations of limited regularity. The simulations show that continuous quasi‐coplanar segments along the backbone are in a range of ≈2–4 repeat units. The ordered orientation and coupling distance of interchain aromatic rings can be correlated with optical properties of materials. A simplified quantum‐mechanical method was developed to investigate optical properties based on MD trajectories. The method was tested to simulate the absorption spectra of four PPV derivatives. The absorption maxima of the calculated spectra are in reasonable agreement with experimental data. This work implies that long‐range electron transfer along the backbones of these polymers may not occur, but may be mediated by interchain interactions.     

On Coarse-Graining by the Inverse Monte Carlo Method: Dissipative Particle Dynamics Simulations Made to a Precise Tool in Soft Matter Modeling

We present a promising coarse-graining strategy for linking micro- and mesoscales of soft matter systems. The approach is based on effective pairwise interaction potentials obtained from detailed atomistic molecular dynamics (MD) simulations, which are then used in coarse-grained dissipative particle dynamics (DPD) simulations. Here, the effective potentials were obtained by applying the inverse Monte Carlo method [Lyubartsev and Laaksonen, Phys. Rev. E. 52, 3730 (1995)] on a chosen subset of degrees of freedom described in terms of radial distribution functions. In our first application of the method, the effective potentials were used in DPD simulations of aqueous NaCl solutions. With the same computational effort we were able to simulate systems of one order of magnitude larger than the MD simulations. The results from the MD and DPD simulations are in excellent agreement.     

Correlation between Dynamic Heterogeneity and Local Structure in a Room‐Temperature Ionic Liquid: A Molecular Dynamics Study of [bmim][PF6]

The complex dynamics of a room‐temperature ionic liquid, 1‐n‐butyl‐3‐methylimidazolium hexafluorophosphate ([bmim][PF₆]), is studied using equilibrium classical molecular dynamics simulations in the temperature range of 250–450 K. The activation energies for the self‐diffusion of ions are around 30–34 kJ mol^?1, with that of the anion a little higher than that for the cation. The electrical conductivity of the liquid is calculated and good agreement with experiments is obtained. Structural relaxation is studied through the decay of coherent (total density–density correlation) and incoherent (self part of density–density correlation) intermediate scattering functions over a range of temperatures and wave vectors relevant to the system. The relaxation data are used to identify and characterize two processes, α and β. The dependence of the two relaxation times on temperature and wave vector is obtained. The dynamical heterogeneity of the ions determined through the non‐Gaussian parameter indicates the motion of the cation to be more heterogeneous than that of the anion. The faster ones among the cations are coordinated to faster anions, while slower cations are surrounded predominantly by slower anions. Thus, the dynamical heterogeneity in this ionic liquid is shown to have structural signatures.     

Infrared Fingerprint Engineering: A Molecular‐Design Approach to Long‐Wave Infrared Transparency with Polymeric Materials

Optical technologies in the long‐wave infrared (LWIR) spectrum (7–14 μm) offer important advantages for high‐resolution thermal imaging in near or complete darkness. The use of polymeric transmissive materials for IR imaging offers numerous cost and processing advantages but suffers from inferior optical properties in the LWIR spectrum. A major challenge in the design of LWIR‐transparent organic materials is that nearly all organic molecules absorb in this spectral window which lies within the so‐called IR‐fingerprint region. We report on a new molecular‐design approach to prepare high refractive index polymers with enhanced LWIR transparency. Computational methods were used to accelerate the design of novel molecules and polymers. Using this approach, we have prepared chalcogenide hybrid inorganic/organic polymers (CHIPs) with enhanced LWIR transparency and thermomechanical properties via inverse vulcanization of elemental sulfur with new organic co‐monomers.     

Pair potentials from diffraction data on liquids: a neural network solution

The inverse theorem of liquids states a one to one correspondence between classical mechanical pair potentials and structural functions. Molecular-dynamics and Monte Carlo simulations provide exact structural functions for known pair interactions. There is no exact or widespread method in the opposite direction, where the pair interactions are to be determined from a priori known pair-correlation functions or structure factors. The methods based on the integral equation theories of liquids are approximate and the iterative refinements of pair potentials with simulations take a long time. We applied artificial neural networks to get pair interactions from known structure factors in this study. We performed molecular-dynamics simulations on one-component systems with different pair potentials and the structure factors were calculated. To optimize (train) the weights of neural networks 2000 pair interaction-structure factor pairs were used. The performance of the method was tested on further 200 data pairs. The method provided reasonable potentials for the majority of the systems opening a "quick and dirty" method to determine pair interactions.     

Systematic implicit solvent coarse graining of dimyristoylphosphatidylcholine lipids

We have used systematic structure‐based coarse graining to derive effective site–site potentials for a 10‐site coarse‐grained dimyristoylphosphatidylcholine (DMPC) lipid model and investigated their state point dependence. The potentials provide for the coarse‐grained model the same site–site radial distribution functions, bond and angle distributions as those computed in atomistic simulations carried out at four different lipid–water molar ratios. It was shown that there is a non‐negligible dependence of the effective potentials on the concentration at which they were generated, which is also manifested in the properties of the lipid bilayers simulated using these potentials. Thus, effective potentials computed at low lipid concentration favor to more condensed and ordered structure of the bilayer with lower average area per lipid, while potentials obtained at higher lipid concentrations provide more fluid‐like structure. The best agreement with the reference data and experiment was achieved using the set of potentials derived from atomistic simulations at 1:30 lipid:water molar ratio providing fully saturated hydration of DMPC lipids. Despite theoretical limitations of pairwise coarse‐grained potentials expressed in their state point dependence, all the resulting potentials provide a stable bilayer structure with correct partitioning of different lipid groups across the bilayer as well as acceptable values of the average lipid area, compressibility and orientational ordering. In addition to bilayer simulations, the model has proven its robustness in modeling of self‐aggregation of lipids from randomly dispersed solution to ordered bilayer structures, bicelles, and vesicles. © 2014 Wiley Periodicals, Inc.     

Molecular dynamics simulations of fluid methane properties using ab initio intermolecular interaction potentials

Intermolecular interaction energy data for the methane dimer have been calculated at a spectroscopic accuracy and employed to construct an ab initio potential energy surface (PES) for molecular dynamics (MD) simulations of fluid methane properties. The full potential curves of the methane dimer at 12 symmetric conformations were calculated by the supermolecule counterpoise‐corrected second‐order Møller‐Plesset (MP2) perturbation theory. Single‐point coupled cluster with single and double and perturbative triple excitations [CCSD(T)] calculations were also carried out to calibrate the MP2 potentials. We employed Pople's medium size basis sets [up to 6‐311++G(3df, 3pd)] and Dunning's correlation consistent basis sets (cc‐pVXZ and aug‐cc‐pVXZ, X = D, T, Q). For each conformer, the intermolecular carbon–carbon separation was sampled in a step 0.1 Å for a range of 3–9 Å, resulting in a total of 732 configuration points calculated. The MP2 binding curves display significant anisotropy with respect to the relative orientations of the dimer. The potential curves at the complete basis set (CBS) limit were estimated using well‐established analytical extrapolation schemes. A 4‐site potential model with sites located at the hydrogen atoms was used to fit the ab initio potential data. This model stems from a hydrogen–hydrogen repulsion mechanism to explain the stability of the dimer structure. MD simulations using the ab initio PES show quantitative agreements on both the atom‐wise radial distribution functions and the self‐diffusion coefficients over a wide range of experimental conditions. © 2008 Wiley Periodicals, Inc. J Comput Chem 2009     

Structural Organization of Water‐Containing Nafion: The Integral Equation Theory

Using the mutually consistent variant of the integral equation polymer reference interaction site model (RISM) theory and the chemically realistic rotational isomeric state (RIS) model, we perform a molecular level modeling of the specific structural organization of perfluorosulfonic acid ionomer (Nafion) with certain amount of physisorbed water. Our results establish molecular scale information necessary for understanding the equilibrium structure and thermodynamics of water‐containing Nafion as well as water distribution and ionic (molecular) transport phenomena in hydrated Nafion membranes. As a first step in this direction, semi‐empirical quantum mechanical calculations on the molecular structures and energies of the polymeric backbone and pendant chains of Nafion with and without additional water molecules are carried out. These data are used in the single‐polymer RIS Monte Carlo simulation, in which the short‐range intramolecular interactions are taken into account via appropriate matrices of statistical weights. The local structure and morphology of the entire bicomponent water‐containing system is calculated consistently on the basis of the RISM theory at different contents of absorbed water. We find that the addition of even a relatively small amount of water leads to the strong intensification of aggregation processes observed for polar sulfonic acid (SO₃H) groups. Water molecules and polar SO₃H groups form mixed aggregates with a three‐layer structure. Incorporation of water molecules inside the aggregates results in an increase of their stability and leads to the increase of the number of associating groups in a stable aggregate. With increasing the number of incorporated solvent molecules, the average size of mixed aggregates increases. The results obtained in the present study support the concept of an irregularly shaped cluster surfaces. Such geometries are favorable to the formation of long channels of connected water‐containing aggregates providing the unique permeability characteristics of Nafion membranes. In addition, we investigate the clustering and continuum percolation in the water phase, using the formalism of pair‐connectedness correlation functions. The mean cluster size found theoretically is in reasonable agreement with the corresponding experimental estimates existing in the literature.     

Nonlinear shear of entangled polymers from nonequilibrium molecular dynamics

This study aims to use molecular dynamics (MD) simulations of Kremer–Grest (KG) chains to inform future developments of models of entangled polymer dynamics. We perform nonequilibrium MD simulations, under shear flow, for well‐entangled KG chains. We study chains of 512 and 1000 KG beads, corresponding to 8 and 15 entanglements, respectively. We compute the linear rheological properties from equilibrium simulations of the stress autocorrelation and obtain from these data the tube model parameters. Under nonlinear shear flow, we compute the shear viscosity, the first and second normal stress differences, and chain contour length. For chains of 512 monomers, we obtain agreement with the results of Cao and Likhtman (ACS Macro. Lett. 2015, 4, 1376). We also compare our nonlinear results with the Graham, Likhtman and Milner‐McLeish (GLaMM) model. We identify some systematic disagreement that becomes larger for the longer chains. We made a comparison of the transient shear stress maximum from our simulations, two nonlinear models and experiments on a wide range of melts and solutions, including polystyrene (PS), polybutadiene, and styrene–butadiene rubber. This comparison establishes that the PS melt data show markedly different behavior to all other melts and solutions and KG simulations reproduce the PS data more closely than either the GLaMM or Xie and Schweizer models. We discuss the performance of these models against the data and simulations. Finally, by imposing a rapid reversing flow, we produce a method to extract the recoverable strain from MD simulations, valid for sufficiently entangled monodisperse polymers. We explore how the resulting data can probe the melt state just before the reversing flow. © 2019 The Authors. Journal of Polymer Science Part B: Polymer Physics published by Wiley Periodicals, Inc. J. Polym. Sci., Part B: Polym. Phys. 2019 , 57, 1692–1704     

Numerical interpretation of molecular surface field in dielectric modeling of solvation

Continuum solvent models, particularly those based on the Poisson‐Boltzmann equation (PBE), are widely used in the studies of biomolecular structures and functions. Existing PBE developments have been mainly focused on how to obtain more accurate and/or more efficient numerical potentials and energies. However to adopt the PBE models for molecular dynamics simulations, a difficulty is how to interpret dielectric boundary forces accurately and efficiently for robust dynamics simulations. This study documents the implementation and analysis of a range of standard fitting schemes, including both one‐sided and two‐sided methods with both first‐order and second‐order Taylor expansions, to calculate molecular surface electric fields to facilitate the numerical calculation of dielectric boundary forces. These efforts prompted us to develop an efficient approximated one‐dimensional method, which is to fit the surface field one dimension at a time, for biomolecular applications without much compromise in accuracy. We also developed a surface‐to‐atom force partition scheme given a level set representation of analytical molecular surfaces to facilitate their applications to molecular simulations. Testing of these fitting methods in the dielectric boundary force calculations shows that the second‐order methods, including the one‐dimensional method, consistently perform among the best in the molecular test cases. Finally, the timing analysis shows the approximated one‐dimensional method is far more efficient than standard second‐order methods in the PBE force calculations. © 2017 Wiley Periodicals, Inc.     

Prediction of protein–protein complexes using replica exchange with repulsive scaling

The realistic prediction of protein–protein complex structures is import to ultimately model the interaction of all proteins in a cell and for the design of new protein–protein interactions. In principle, molecular dynamics (MD) simulations allow one to follow the association process under realistic conditions including full partner flexibility and surrounding solvent. However, due to the many local binding energy minima at the surface of protein partners, MD simulations are frequently trapped for long times in transient association states. We have designed a replica-exchange based scheme employing different levels of a repulsive biasing between partners in each replica simulation. The bias acts only on intermolecular interactions based on an increase in effective pairwise van der Waals radii (repulsive scaling (RS)-REMD) without affecting interactions within each protein or with the solvent. For a set of five protein test cases (out of six) the RS-REMD technique allowed the sampling of near-native complex structures even when starting from the opposide site with respect to the native binding site for one partner. Using the same start structures and same computational demand regular MD simulations sampled near native complex structures only for one case. The method showed also improved results for the refinement of docked structures in the vicinity of the native binding geometry compared to regular MD refinement.     

Molecular dynamics study of the interface between water and 2-nitrophenyl octyl ether

We present results of molecular dynamics simulations of the interface between water and 2-nitrophenyl octyl ether (NPOE). This system is analyzed in detail using a procedure to calculate intrinsic profiles of several important properties (density, radial distribution functions, hydrogen bonds, molecular orientation, self-diffusion). The interface was found to be molecularly sharp but corrugated by thermal fluctuations. Using a method based on capillary wave theory, we have estimated the interfacial tension and obtained good agreement with values calculated from the virial route. The results were compared to simulations of the water/nitrobenzene interface. The presence of an alkyl chain in NPOE introduces an added degree of hydrophobicity, which causes an increase in the interfacial tension. Furthermore, interfacial NPOE molecules are less organized than nitrobenzene and show a distinct dynamic response. These results shed light on the observed differences between these two organic liquids in electrochemical studies.     

Length‐Dependent Conductance of Molecular Wires and Contact Resistance in Metal–Molecule–Metal Junctions

Molecular wires are covalently bonded to gold electrodes—to form metal–molecule–metal junctions—by functionalizing each end with a ? SH group. The conductance of a wide variety of molecular junctions is studied theoretically by using first‐principles density functional theory (DFT) combined with the nonequilibrium Green′s function (NEGF) formalism. Based on the chain‐length‐dependent conductance of the series of molecular wires, the attenuation factor β is obtained and compared with the experimental data. The β value is quantitatively correlated to the molecular HOMO–LUMO gap. Coupling between the metallic electrode and the molecular bridge plays an important role in electron transport. A contact resistance of 6.0±2.0 KΩ is obtained by extrapolating the molecular‐bridge length to zero. This value is of the same magnitude as the quantum resistance.     

Reaction energetics on long‐range corrected density functional theory: Diels–Alder reactions

The possibility of quantitative reaction analysis on the orbital energies of long‐range corrected density functional theory (LC‐DFT) is presented. First, we calculated the Diels–Alder reaction enthalpies that have been poorly given by conventional functionals including B3LYP functional. As a result, it is found that the long‐range correction drastically improves the reaction enthalpies. The barrier height energies were also computed for these reactions. Consequently, we found that dispersion correlation correction is also crucial to give accurate barrier height energies. It is, therefore, concluded that both long‐range exchange interactions and dispersion correlations are essentially required in conventional functionals to investigate Diels–Alder reactions quantitatively. After confirming that LC‐DFT accurately reproduces the orbital energies of the reactant and product molecules of the Diels–Alder reactions, the global hardness responses, the halves of highest occupied molecular orbital (HOMO)‐lowest unoccupied molecular orbital (LUMO) energy gaps, along the intrinsic reaction coordinates of two Diels–Alder reactions were computed. We noticed that LC‐DFT results satisfy the maximum hardness rule for overall reaction paths while conventional functionals violate this rule on the reaction pathways. Furthermore, our results also show that the HOMO‐LUMO gap variations are close to the reaction enthalpies for these Diels–Alder reactions. Based on these results, we foresee quantitative reaction analysis on the orbital energies. © 2012 Wiley Periodicals, Inc.