OOPSE: an object-oriented parallel simulation engine for molecular dynamics

OOPSE is a new molecular dynamics simulation program that is capable of efficiently integrating equations of motion for atom types with orientational degrees of freedom (e.g. "sticky" atoms and point dipoles). Transition metals can also be simulated using the embedded atom method (EAM) potential included in the code. Parallel simulations are carried out using the force-based decomposition method. Simulations are specified using a very simple C-based meta-data language. A number of advanced integrators are included, and the basic integrator for orientational dynamics provides substantial improvements over older quaternion-based schemes.     

MILCH SHAKE: An efficient method for constraint dynamics applied to alkanes

We describe a method to impose constraints in a molecular dynamics simulation. A technique developed to solve the special case of a linear topology (MILC SHAKE) is hybridized with the SHAKE algorithm. The methodology, which we term MILC‐hybridized SHAKE (or MILCH SHAKE), applies to more complex topologies. Here we consider the important case of all atom models of alkanes. Exploiting the mass difference between carbon and hydrogen we show that for higher alkanes MILCH SHAKE can be an order of magnitude faster than SHAKE. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2009     

金属Cu晶化与玻璃化行为的分子动力学模拟

A series of simulations of the crystallization and vitrification processes for metal Cu were carried out by means of the molecular dynamics technique. The radial distribution function, common neighbors, internal energy and volume of the system were recorded during the processes. The atomic internal energy, atomic Voronoi volume and atomic stress field of the relax system were analyzed at zero temperature. The interaction between atoms in the system is described using the embedded atom potential as proposed by Mishin. The simulation results show that crystalline and non-crystalline phases form at lower and higher cooling rates respectively. In comparison to nanocrystals, it is found that metallic glass has higher internal energy and larger volume. The intrinsic stress field is induced by distortion of the lattice.     

Solvated ensemble averaging in the calculation of partial atomic charges

In the calculation of partial atomic charges, for use in molecular mechanics or dynamics simulations, it is common practice to select only a single conformation for the molecule of interest. For molecules that contain rotatable bonds, it is preferable to compute the charges from several relevant conformations. We present here results from a charge derivation protocol that determines the partial charges by averaging charges computed for conformations selected from explicitly solvated MD simulations, performed under periodic boundary conditions. This approach leads to partial charges that are weighted by a realistic population of conformations and that are suitable for condensed phase simulations. This protocol can, in principle, be applied to any class of molecule and to nonaqueous solvation. Carbohydrates contain numerous hydroxyl groups that exist in an ensemble of orientations in solution, and in this report we apply ensemble averaging to a series of methyl glycosides. We report the extent to which ensemble averaging leads to charge convergence among the various monosaccharides and among the constituent atoms within a given monosaccharide. Due to the large number of conformations (200) in our ensembles, we are able to compute statistically relevant standard deviations for the partial charges. An analysis of the standard deviations allows us to assess the extent to which equivalent atom types may, nevertheless, require unique partial charges. The configurations of the hydroxyl groups exert considerable influence on internal energies, and the limits of ensemble averaged charges are discussed in terms of these properties.     

A comparison of neighbor search algorithms for large rigid molecules

Fast determination of neighboring atoms is an essential step in molecular dynamics simulations or Monte Carlo computations, and there exists a variety of algorithms to efficiently compute neighbor lists. However, most of these algorithms are general, and not specifically designed for a given type of application. As a result, although their average performance is satisfactory, they might be inappropriate in some specific application domains. In this article, we study the case of detecting neighbors between large rigid molecules, which has applications in, e.g., rigid body molecular docking, Monte Carlo simulations of molecular self-assembly or diffusion, and rigid body molecular dynamics simulations. More precisely, we compare the traditional grid-based algorithm to a series of hierarchy-based algorithms that use bounding volumes to rapidly eliminate large groups of irrelevant pairs of atoms during the neighbor search. We compare the performance of these algorithms based on several parameters: the size of the molecules, the average distance between them, the cutoff distance, as well as the type of bounding volume used in the culling hierarchy (AABB, OBB, wrapped, or layered spheres). We demonstrate that for relatively large systems (> 100,000 atoms) the algorithm based on the hierarchy of wrapped spheres shows the best results and the traditional grid-based algorithm gives the worst timings. For small systems, however, the grid-based algorithm and the one based on the wrapped sphere hierarchy are beneficial.     

YinYang atom: a simple combined ab initio quantum mechanical molecular mechanical model

A simple interface is proposed for combined quantum mechanical (QM) molecular mechanical (MM) calculations for the systems where the QM and MM regions are connected through covalent bonds. Within this model, the atom that connects the two regions, called YinYang atom here, serves as an ordinary MM atom to other MM atoms and as a hydrogen-like atom to other QM atoms. Only one new empirical parameter is introduced to adjust the length of the connecting bond and is calibrated with the molecule propanol. This model is tested with the computation of equilibrium geometries and protonation energies for dozens of molecules. Special attention is paid on the influence of MM point charges on optimized geometry and protonation energy, and it is found that it is important to maintain local charge-neutrality in the MM region in order for the accurate calculation of the protonation and deprotonation energies. Overall the simple YinYang atom model yields comparable results to some other QM/MM models.     

The absorption and diffusion of polyethylene chains on the carbon nanotube: The molecular dynamics study

Classical molecular dynamics simulations have been used to investigate the absorption and diffusion behavior of polyethylene (PE) chains on the surface of the side‐wall of the carbon nanotube (CNT). Different degrees of polymerization from 50 to 80 at separate temperatures of 300, 400, 500, and 600 K are considered. Through the simulation, it is examined that the PE chains are absorbed on the surface of CNT and form stable composites with the nanotube as capsules. It is found that the most probable distance between the CNT and the C atoms in backbone of PE molecules only attribute to the temperature, and at T = 300 K, this distance is about 3.8 Å. Furthermore, the pattern of the composites mainly depends on the temperature and the length matching of the chains and the CNT. In particular, the PE chains keep approximately linear conformation, and extend along the axis of the CNT at the room temperature. © 2007 Wiley Periodicals, Inc. J Polym Sci Part B: Polym Phys 46: 272–280, 2008     

投影解析法求直链烷烃内旋转时非键原子对间距离

用旋转投影解析几何的方法,计算了链烃分子链上相邻原子上的非键原子对间的距离r(φ_i),及相间原子上的非键原子对间的距离r(φ_i,φ₁₊₁.所得公式适用于链上不含双键的链状分子.     

Atom and Bond Fukui Functions and Matrices: A Hirshfeld‐I Atoms‐in‐Molecule Approach

The Fukui function is often used in its atom‐condensed form by isolating it from the molecular Fukui function using a chosen weight function for the atom in the molecule. Recently, Fukui functions and matrices for both atoms and bonds separately were introduced for semiempirical and ab initio levels of theory using Hückel and Mulliken atoms‐in‐molecule models. In this work, a double partitioning method of the Fukui matrix is proposed within the Hirshfeld‐I atoms‐in‐molecule framework. Diagonalizing the resulting atomic and bond matrices gives eigenvalues and eigenvectors (Fukui orbitals) describing the reactivity of atoms and bonds. The Fukui function is the diagonal element of the Fukui matrix and may be resolved in atom and bond contributions. The extra information contained in the atom and bond resolution of the Fukui matrices and functions is highlighted. The effect of the choice of weight function arising from the Hirshfeld‐I approach to obtain atom‐ and bond‐condensed Fukui functions is studied. A comparison of the results with those generated by using the Mulliken atoms‐in‐molecule approach shows low correlation between the two partitioning schemes.     

Cluster distance geometry of polypeptide chains

Distance geometry has been a broadly useful tool for dealing with conformational calculations. Customarily each atom is represented as a point, constraints on the distances between some atoms are obtained from experimental or theoretical sources, and then a random sampling of conformations can be calculated that are consistent with the constraints. Although these methods can be applied to small proteins having on the order of 1000 atoms, for some purposes it is advantageous to view the problem at lower resolution. Here distance geometry is generalized to deal with distances between sets of points. In the end, much of the same techniques produce a sampling of different configurations of these sets of points subject to distance constraints, but now the radii of gyration of the different sets play an important role. A simple example is given of how the packing constraints for polypeptide chains combine with loose distance constraints to give good calculated protein conformers at a very low resolution.     

Balancing an accurate representation of the molecular surface in generalized born formalisms with integrator stability in molecular dynamics simulations

Different integrator time steps in NVT and NVE simulations of protein and nucleic acid systems are tested with the GBMV (Generalized Born using Molecular Volume) and GBSW (Generalized Born with simple SWitching) methods. The simulation stability and energy conservation is investigated in relation to the agreement with the Poisson theory. It is found that very close agreement between generalized Born methods and the Poisson theory based on the commonly used sharp molecular surface definition results in energy drift and simulation artifacts in molecular dynamics simulation protocols with standard 2-fs time steps. New parameters are proposed for the GBMV method, which maintains very good agreement with the Poisson theory while providing energy conservation and stable simulations at time steps of 1 to 1.5 fs.     

Fragmentation at sp2 carbon atoms in fragment molecular orbital method

In the fragment molecular orbital (FMO) method, a given molecular system is usually fragmented at sp³ carbon atoms. However, fragmentation at different sites sometimes becomes necessary. Hence, we propose fragmentation at sp² carbon atoms in the FMO method. Projection operators are constructed using sp² local orbitals. To maintain practical accuracy, it is essential to consider the three-body effect. In order to suppress the corresponding increase of computational cost, we propose approximate models considering local trimers. Numerical verification shows that the present models are as accurate as or better than the standard FMO2 method in total energy with fragmentation at sp³ carbon atoms.     

Strategies to model the near-solute solvent molecular density/polarization

The solvent molecular distribution significantly affects the behavior of the solute molecules and is thus important in studying many biological phenomena. It can be described by the solvent molecular density distribution, g, and the solvent electric dipole distribution, p. The g and p can be computed directly by counting the number of solvent molecules/dipoles in a microscopic volume centered at r during a simulation or indirectly from the mean force F and electrostatic field E acting on the solvent molecule at r, respectively. However, it is not clear how the g and p derived from simulations depend on the solvent molecular center or the solute charge and if the g(F) and p(E) computed from the mean force and electric field acting on the solvent molecule, respectively, could reproduce the corresponding g and p obtained by direct counting. Hence, we have computed g, p, g(F), and p(E) using different water centers from simulations of a solute atom of varying charge solvated in TIP3P water. The results show that g(F) and p(E) can reproduce the g and p obtained using a given count center. This implies that rather than solving the coordinates of each water molecule by MD simulations, the distribution of water molecules could be indirectly obtained from analytical formulas for the mean force F and electrostatic field E acting on the solvent molecule at r. Furthermore, the dependence of the g and p distributions on the solute charge revealed provides an estimate of the change in g and p surrounding a biomolecule upon a change in its conformation.     

Reaching multi‐nanosecond timescales in combined QM/MM molecular dynamics simulations through parallel horsetail sampling

We report an enhanced sampling technique that allows to reach the multi‐nanosecond timescale in quantum mechanics/molecular mechanics molecular dynamics simulations. The proposed technique, called horsetail sampling, is a specific type of multiple molecular dynamics approach exhibiting high parallel efficiency. It couples a main simulation with a large number of shorter trajectories launched on independent processors at periodic time intervals. The technique is applied to study hydrogen peroxide at the water liquid–vapor interface, a system of considerable atmospheric relevance. A total simulation time of a little more than 6 ns has been attained for a total CPU time of 5.1 years representing only about 20 days of wall‐clock time. The discussion of the results highlights the strong influence of the solvation effects at the interface on the structure and the electronic properties of the solute. © 2017 Wiley Periodicals, Inc.     

Structural transformations in silver nanoparticles

A molecular dynamics simulation was performed for silver clusters of 147, 309, and 561 atoms with the initial cuboctahedral habit in the temperature range 0–1000 K with an embedded atom potential for silver. Structural transitions of the silver clusters to complex twins (icosahedral habit) with coherent (111)/(111) boundaries over all edges of icosahedra were found, which started at temperatures of 50 K, 350 K, and 700 K, respectively. To analyze the structural transformations in nanoparticles, an algorithm is proposed based on a simplicial Delaunay decomposition (Delaunay triangulation). It was found that after the transition of silver nanoparticles to complex twins, the atomic motion becomes vibrational; the atoms vibrate around the sites that correspond to the vertices of the regular polyhedra. In the case of the 147-atom silver nanoparticle, the polyhedra are arranged in the following sequence, starting from the center of mass: icosahedron (12 atoms), icosododecahedron (30 atoms), icosahedron (12 atoms), dodecahedron (20 atoms), truncated icosahedron (60 atoms, isostructural with fullerene C₆₀), icosahedron (12 atoms), and one atom at the center of mass.     

Applications of weighting and chirality strategies for distance geometry algorithms to an enterotoxin peptide analog

Significant improvements are made to a recent algorithm that finds molecular conformations using distance geometry on nuclear magnetic resonance data. Weighting factors for the nearest approximation of the distance matrix to a data matrix are allowed to vary between iterations of the algorithm. These changes are proportional to the error of the distance between atoms in the configuration and the nuclear magnetic resonance data bounds. The weight changes increase the rate of convergence by an order of magnitude. Penalty functions are proposed to ensure the correct chirality. Numerical results for these modifications and subsequent energy calculations using CHARMm are given for an analog of the heat stable (ST) enterotoxin peptide STh produced byE. coli in humans.     

A method for conformational sampling of loops in proteins based on adiabatic decoupling and temperature or force scaling

A method for conformational Boltzmann sampling of loops in proteins in aqueous solution is presented that is based on adiabatic decoupling molecular dynamics (MD) simulation with temperature or force scaling. To illustrate the enhanced sampling, the loop from residues 33 to 43 in the bovine protein ribonuclease A is adiabatically decoupled from the rest of the protein and the solvent with a mass scaling factor s(m) =1000 and the sampling is enhanced with a scaling of the temperature using s(T) =2 or of the force using s(V) =0.667. Over 5 ns of simulation the secondary structure of the protein remains unaltered while a combined dihedral-angle conformational cluster analysis shows an increase of conformations outside the first most populated cluster of loop conformations for adiabatic decoupling MD with temperature scaling using s(T) =2 or force scaling using s(V) =0.667 compared to the standard MD simulation. The atom-positional root-mean-square fluctuations of the C(α) atoms of the loop show an increase in the movement of the loop as well, indicating that adiabatic decoupling MD with upscaling of the temperature or downscaling of the force is a promising method for conformational Boltzmann sampling.     

Molecular design of novel indacenodithiophene-based organic dyes for efficient dye-sensitized solar cells applications

Indacenodithiophene (IDT)-based high-efficiency photovoltaics have received increasing attention recently. This paper reports a density functional theory investigation of the electronic and optical properties of three IDT-based organic dyes together with the dye/(TiO₂)₄₆ interface. In order to enhance the photoelectric properties of IDT dyes, this paper considers two methods for the structure modification of the experimentally reported dye DPInDT (J. Org. Chem. 2011, 76, 8977): the extension of the conjugation length by dithienothiophene as well as the heteroatom substitution of the bridging atoms by electron-rich nitrogen atoms. Our calculations show that both methods obviously affect the distributions of the molecular orbitals and notably red shift the absorption peaks of around 20 nm, with the former method demonstrating enhanced light harvesting efficiency. The structure modifications proposed also enhance the emission spectrum properties for IDT-based organic dyes. The calculated ultrafast injection time of electrons from the excited state of IDT dyes to the (TiO₂)₄₆ belongs to the femtosecond order of magnitude, and is ideal for efficient photoelectric conversion process in dye-sensitized solar cells (DSSCs) applications. The IDT dyes designed in this paper have good electronic and spectroscopic properties. This study is expected to provide useful guidance for the development of novel IDT dyes for applications in DSSCs.     

聚丙烯腈PAN分子链的第一性原理研究

The equilibrium geometries of the polyacrylonitrile (PAN) chain was theoretically studied using the Hartree-Fock method at the STO-3G levels. As for the optimized structures, the average distance of the C atom couple in the main chain is 155.6 pm; the average distance of the C atom couple in the branch chain is 149.7 pm; the average distance between a C atom bonding with N atom is 115.5 pm. For the charge distribution, because of the influence of a N atom with its comparatively larger negative charge, the C atoms in the main chain are different in their charge distribution. Finally, the vibration models of the chain have been analyzed to clarify the reaction sequence of dehydrogenation and cyclization during pre oxidation and carbonization of the polyacrylonitrile.     

Density Functional Theory Study of Pd Aggregation on a Pyridine-Terminated Self-Assembled Monolayer

By using density functional theory calculations, the initial steps towards Pd metal cluster formation on a pyridine-terminated self-assembled monolayer (SAM) consisting of 3-(4-(pyridine-4-yl)phenyl)propane-1-thiol on an Au(1 1 1) surface are investigated. Theoretical modelling allows the investigation of structural details of the SAM surface and the metal/SAM interface at the atomic level, which is essential for elucidating the nature of Pd–SAM and Pd–Pd interactions at the liquid/solid interface and gaining insight into the mechanism of metal nucleation in the initial stage of electrodeposition. The structural flexibility of SAM molecules was studied first and the most stable conformation was identified, planar molecules in a herringbone packing, as the model for Pd adsorption. Two binding sites are found for Pd atoms on the pyridine end group of the SAM. The strong interaction between Pd atoms and pyridines illustrates the importance of SAM functionalisation in the metal nucleation process. Consistent with an energetic driving force of approximately −0.3 eV per Pd atom towards Pd aggregation suggested by static calculations, a spontaneous Pd dimerisation is observed in ab initio molecular dynamic studies of the system. Nudged elastic band calculations suggest a potential route with a low energy barrier of 0.10 eV for the Pd atom diffusion and then dimerisation on top of the SAM layer.