Calculations of a list of neighbors in Molecular Dynamics simulations

A new algorithm is proposed for the evaluation of nonbonded interactions in Molecular Dynamics simulations. The algorithm is based on a grid search and on partitioning of the atoms into boxes rather than on calculations of distances. The effort associated with the generation of the box list grows only linearly with the number of atoms. The algorithm is particularly advantageous for solvated systems. Test calculations show significant savings in CPU time and storage compared to commonly used algorithms for systems containing in excess of ca. 600 atoms.     

Macromolecular dynamics on a shared-memory multiprocessor

This is a study of the effect of data structures and algorithms on parallelism and vectorization for the molecular dynamics package GROMOS. About 80% of the computing time (for the GROMOS benchmark) is spent on the nonbonded interactions, and an important issue is the method that is used for finding all pairs of atoms (actually atom groups) which are within a given cutoff radius of each other. Favored by the authors of GROMOS is the use of the (Verlet) neighbor list with direct calculation of all distances. Substantial restructuring of their code for the nonbonded interactions has resulted in a speedup of a factor 6.9 for the eight-processor Alliant FX/8 through the use of concurrent subroutine calls. Vectorization gave a further improvement of 30%. However, the direct calculation of all distances does not scale up linearly with the number of atoms. (Nor would neighbor lists be needed if cutoffs were abandoned in favor of such techniques as the fast multipole method). Therefore an alternative GROMOS subroutine that implements the “gridcell plus linked list” approach without a neighbor list was also rewritten and tested. This resulted in the discovery of an apparent flaw in an early study of the relative merits of the two algorithms. It is observed that actually the neighbor list version of GROMOS is not much faster.     

Microscopic Symmetry Imposed by Rotational Symmetry Boundary Conditions in Molecular Dynamics Simulation

A large number of viral capsids, as well as other macromolecular assemblies, have icosahedral structure or structures with other rotational symmetries. This symmetry can be exploited during molecular dynamics (MD) to model in effect the full viral capsid using only a subset of primary atoms plus copies of image atoms generated from rotational symmetry boundary conditions (RSBC). A pure rotational symmetry operation results in both primary and image atoms at short range, and within nonbonded interaction distance of each other, so that nonbonded interactions can not be specified by the minimum image convention and explicit treatment of image atoms is required. As such an unavoidable consequence of RSBC is that the enumeration of nonbonded interactions in regions surrounding certain rotational axes must include both a primary atom and its copied image atom, thereby imposing microscopic symmetry for some forces. We examined the possibility of artifacts arising from this imposed microscopic symmetry of RSBC using two simulation systems: a water shell and human rhinovirus 14 (HRV14) capsid with explicit water. The primary unit was a pentamer of the icosahedron, which has the advantage of direct comparison of icosahedrally equivalent spatial regions, for example regions near a 2-fold symmetry axis with imposed symmetry and a 2-fold axis without imposed symmetry. Analysis of structural and dynamic properties of water molecules and protein atoms found similar behavior near symmetry axes with imposed symmetry and where the minimum image convention fails compared with that in other regions in the simulation system, even though an excluded volume effect was detected for water molecules near the axes with imposed symmetry. These results validate the use of RSBC for icosahedral viral capsids or other rotationally symmetric systems.     

Calculation of the minimum energy conformation of acetylcholine using a global optimization technique

The conformational energy of acetylcholine is minimized with respect to the distances between nonbonded atoms with the help of the Bremermann method of unconstrained global optimization. The set of distances for which the energy is the absolute minimum is then used to calculate the coordinates of all the atoms and hence the conformation of the molecule. The simplest type of potential function, namely the classical potential function is chosen for the calculation. The major advantages of this method are (i) that the starting point need not be close to the actual solution, (ii) that it gives the global minimum, irrespective of the starting point, (iii) that it is very general and can be used for any type of potential function, and (iv) that it does not require the computation of gradients. The results obtained are in very good agreement with those of other workers.     

Investigation of the distinction between van der Waals interaction and chemical bonding based on the PAEM‐MO diagram

In recent years, the basic problem of understanding chemical bonding, nonbonded, and/or van der Waals interactions has been intensively debated in terms of various theoretical methods. We propose and construct the potential acting on one electron in a molecule‐molecular orbital (PAEM‐MO) diagram, which draws the PAEM inserted the MO energy levels with their major atomic orbital components. PAEM‐MO diagram is able to show clear distinction of chemical bonding from nonbonded and/or vdW interactions. The rule for this is as follows. Along the line connecting two atoms in a molecule or a complex, the existence of chemical bonding between these two atoms needs to satisfy two conditions: (a) a critical point of PAEM exists and (b) PAEM barrier between the two atoms is lower in energy than the occupied major valence‐shell bonding MO which contains in‐phase atomic components (positive overlap) of the two considered atoms. In contrast to the chemical bonding, for a nonbonded interaction or van der Waals interaction between two atoms, both conditions (a) and (b) do not be satisfied at the same time. This is demonstrated and discussed by various typical cases, particularly those related to helium atom and H? H bonding in phenanthrene. There are helium bonds in HHeF and HeBeO molecules, whereas no H? H bonding in phenanthrene. The validity and limitation for this rule is demonstrated through the investigations of the curves of the PAEM barrier top and MO energies versus the internuclear distances for He₂, H₂, and He₂⁺ systems. © 2014 Wiley Periodicals, Inc.     

Linear-scaling soft-core scheme for alchemical free energy calculations

Alchemical free energy calculations involving the removal or insertion of atoms into condensed phase systems generally make use of soft-core scaling of nonbonded interactions, designed to circumvent numerical instabilities that arise from weakly interacting "hard" atoms in close proximity. Current methods model soft-core atoms by introducing a nonlinear dependence between the shape of the interaction potential and the strength of the interaction. In this article, we propose a soft-core method that avoids introducing such a nonlinear dependence, through the application of a smooth flattening of the potential energy only in a region that is energetically accessible under normal conditions. We discuss the benefits that this entails and explore a selection of applications, including enhanced methods for the estimation of free energy differences and for the automated optimization of the placement of intermediate states in multistage alchemical calculations.     

Conformational energetics of 1,3-dichloropropane as predicted by several calculations methods

Two semiempirical methods (MNDO and AM1), a molecular mechanics technique (MM2) and two ab initio approaches (6–31G* full optimization and 3–21G/6–31G*) were used to calculate the ordering of and energy difference between conformers in 1,3-dichloropropane. The semiempirical methods did not order the conformers properly or predict correct energy differences. Both ab initio methods ordered the conformers and predicted energy differences correctly, with the 6–31G* full optimization performing slightly better. The MM2 results were presented for calculations involving a force field with no hydrogens and a full force field of all atoms. The full force field properly ordered the conformers but did not correctly predict the energy differences. The nonhydrogen field ordered the conformers based on the Cl…Cl nonbonded distance. The data show that conformer stability is not a simple matter of maximizing the Cl…Cl nonbonded distance, but is also related to some other stabilizing interaction(s).     

Rigid-CLL: avoiding constant-distance computations in cell linked-lists algorithms

Many of the existing molecular simulation tools require the efficient identification of the set of nonbonded interacting atoms. This is necessary, for instance, to compute the energy values or the steric contacts between atoms. Cell linked-lists can be used to determine the pairs of atoms closer than a given cutoff distance in asymptotically optimal time. Despite this long-term optimality, many spurious distances are anyway computed with this method. Therefore, several improvements have been proposed, most of them aiming to refine the volume of influence for each atom. Here, we suggest a different improvement strategy based on avoiding to fill cells with those atoms that are always at a constant distance of a given atom. This technique is particularly effective when large groups of the particles in the simulation behave as rigid bodies as it is the case in simplified models considering only few of the degrees of freedom of the molecule. In these cases, the proposed technique can reduce the number of distance computations by more than one order of magnitude, as compared with the standard cell linked-list technique. The benefits of this technique are obtained without incurring in additional computation costs, because it carries out the same operations as the standard cell linked-list algorithm, although in a different order. Since the focus of the technique is the order of the operations, it might be combined with existing improvements based on bounding the volume of influence for each atom.     

On searching neighbors in computer simulations of macromolecular systems

Various algorithms for evaluating nonbonded interactions in molecular dynamics (MD ) simulations of macromolecular systems are considered, and a combination of two techniques using a space grid for finding neighbor atoms is proposed. The application of grid search techniques to nonrectangular periodic systems is discussed. Finally, the computing time required by different algorithms is compared on a Cray-1 vector processing computer as well as on a CDC Cyber 170/760 sequential computer. Neighbor list techniques turn out to be faster than the grid search techniques for the systems considered here. However, storage requirements may exclude the use of neighbor list techniques for large systems. Finally, MD of complex macromolecular systems turned out to be about 4–14 times faster on a Cray-1 than on a Cyber 170/760.     

Hydrogen-hydrogen bonding in planar biphenyl, predicted by atoms-in-molecules theory, does not exist

Based on an Atoms-in-Molecules (AIM) analysis, Matta et al. recently claimed evidence for the existence of hydrogen-hydrogen bonding between ortho-hydrogen atoms, pointing towards each other from adjacent phenyl groups in planar biphenyl. This AIM result is opposed to the classical view that nonbonded steric repulsion between the ortho-hydrogen atoms is responsible for the higher energy of the planar as compared to the twisted geometry of biphenyl. In the present work, we address the question if hydrogen-hydrogen bonding in biphenyl exists, as suggested by AIM, or not. To this end, we have analyzed the potential energy surface for internal rotation of biphenyl in terms of two interacting phenyl radicals using density functional theory (DFT) at BP86/TZ2P. A detailed analysis of the bonding mechanism and a quantitative bond energy decomposition in the framework of Kohn-Sham DFT show that Pauli (or overlap) repulsion, mainly between C(ortho)--H(ortho) phenyl MOs, prevents biphenyl from being planar and forces it to adopt a twisted equilibrium geometry. Furthermore, a derivative of biphenyl in which all four ortho-hydrogen atoms have been removed does adopt a planar equilibrium geometry. Thus, our results confirm the classical view of steric repulsion between ortho-hydrogen atoms in biphenyl and they falsify the hypothesis of hydrogen-hydrogen bonding.     

Conformational sensitivity of polyether macrocycles to electrostatic potential: Partial atomic charges,molecular mechanics,and conformational prediction

Molecular recognition (whether by enzymes, the immune system, or chelating ligands) depends critically on molecular conformation. Molecular mechanics predicts energetically favorable molecular conformations by locating low energy conformations using an empirical fit of molecular potential energy as a function of internal coordinates. Molecular mechanics analysis of 18-crown-6 demonstrates that the nonbonded term (primarily the electrostatic part) is the largest contributor to the conformational energy. Nevertheless, common methods of treating the electrostatic interaction for 18-crown-6 yield inconsistent values for conformational energies partly because partial charges assigned to each atom can change with conformation due to through-space inductive effects which are not considered in most molecular mechanics programs. Similar findings from several other groups are reviewed to support our conclusions. We argue for care and caution in predicting conformational preferences of molecules with two or more highly polar atoms. We also discuss the desirability of using an empirical method of partial charge determination such as the charge equilibration algorithm of Rappé and Goddard (or a suitable generalization which includes polarization) as a method of including these effects in molecular mechanics and molecular dynamics calculations.     

Short nonbonded contact distances in organic molecules and their use as atom-clash criteria in conformer validation and searching

Short, intramolecular nonbonded contact distances from a large sample of organic molecules retrieved from the Cambridge Structural Database (CSD) have been analyzed. With the exception of the element pairs N···S, O···P, O···S, and S···S, the first percentiles of X···Y distance distributions (X, Y = C, Br, Cl, F, N, O, P, S) are well estimated by ∑vdw -0.5 ?, where ∑vdw is the sum of the Bondi van der Waals radii. The 0.1th percentiles are typically a further 0.1 ? shorter. Some 99% of well-refined organic molecules from the CSD have no nonbonded contacts shorter than the 0.1th percentile and no more than two contacts shorter than the first percentile. This can be used as the basis of an atom-clash test for validating less precise crystal-structure data, such as the geometries of protein ligands. In principle, the same test can be used in molecular modeling to identify and filter out unacceptable conformations generated in a conformational search. This is complicated by the fact that conformer generation is usually performed on molecular models with standard bond angles that are not relaxed during the search. In consequence, conformations often appear to contain untenable nonbonded contacts, which would, however, be removed by bond-angle relaxation. This is particularly likely for molecules containing conjugating substituents bonded to adjacent atoms of an aromatic ring, or on the same side of a double bond. Other molecules particularly likely to be affected are those containing rings or other bulky groups separated by a single-atom linkage, and those with the capacity to form intramolecular hydrogen bonds. The problem is greatly ameliorated by the fact that there are many ways to approximate a true conformation, leading to an increased probability that at least one of the approximations will satisfy atom-clash criteria.     

A novel decomposition of torsional potentials into pairwise interactions: A study of energy second derivatives

A general method of analyzing intramolecular torsional potentials in terms of energy second derivatives that couple the rotating atoms is presented. The method offers a rigorous decomposition of the total torsional potential into pairwise (dihedral) interactions and enables one to derive nonbonded torsional interactions between 1–4 atoms as well as between more distant atoms and sites. The method is demonstrated on ethane, propane and acetaldehyde. It is shown that the 1–4 H…H dihedral potentials in ethane and propane are very similar, thereby supporting the notion of transferable force field potential functions. However, the dihedral potentials that are obtained differ from 1–4 potentials that are used in current force fields. Intramolecular three body effects are clearly seen in this method and are found to be relatively large for the dihedral interactions, although in the one case studied (propane) the overall effect on the methyl-methyl interaction is negligible due to cancellation of terms. The analysis explicitly shows that the barrier in acetaldehyde is due mainly to the dihedral H…H interaction.     

Hydrogen-hydrogen bonding: a stabilizing interaction in molecules and crystals

Bond paths linking two bonded hydrogen atoms that bear identical or similar charges are found between the ortho-hydrogen atoms in planar biphenyl, between the hydrogen atoms bonded to the C1-C4 carbon atoms in phenanthrene and other angular polybenzenoids, and between the methyl hydrogen atoms in the cyclobutadiene, tetrahedrane and indacene molecules corseted with tertiary-tetra-butyl groups. It is shown that each such H-H interaction, rather than denoting the presence of "nonbonded steric repulsions", makes a stabilizing contribution of up to 10 kcal mol(-1) to the energy of the molecule in which it occurs. The quantum theory of atoms in molecules-the physics of an open system-demonstrates that while the approach of two bonded hydrogen atoms to a separation less than the sum of their van der Waals radii does result in an increase in the repulsive contributions to their energies, these changes are dominated by an increase in the magnitude of the attractive interaction of the protons with the electron density distribution, and the net result is a stabilizing change in the energy. The surface virial that determines the contribution to the total energy decrease resulting from the formation of the H-H interatomic surface is shown to account for the resulting stability. It is pointed out that H-H interactions must be ubiquitous, their stabilization energies contributing to the sublimation energies of hydrocarbon molecular crystals, as well as solid hydrogen. H-H bonding is shown to be distinct from "dihydrogen bonding", a form of hydrogen bonding with a hydridic hydrogen in the role of the base atom.     

Systematic drug receptor mapping: A new approach to the analysis of conformational energy calculations of flexible molecules with application to dopaminergic and adrenergic agonists

The computer program PRODIS is used to find low energy conformations of flexible molecules by searching the potential energy surface(s) of one or more torsion angles via rigid rotation. The n-dimensional grid of energy versus torsion angles is then converted to a Boltzman probability distribution, with the probability being represented not as a function of torsion angle, but rather a distance between two atoms. These atoms are chosen by comparison with a known, active analogue in which certain atoms have previously been determined as requirements for drug activity. PRODIS produces a list of low energy conformations, their corresponding interatomic distances and the Boltzman probability for each distance ±0.125, as well as the total probability for each conformation. The user also specifies a target interatomic distance and range (usually derived from a more rigid analogue) for which PRODIS lists all conformations and their Boltzman probability that meet this distance.     

Conformational and energetic effects of truncating nonbonded interactions in an aqueous protein dynamics simulation

In a previous aqueous protein dynamics study, we compared the rms deviation relative to the crystal structure for distance-dependent and constant dielectric models with and without a nonbonded cutoff. The structures obtained from a constant dielectric simulation with a cutoff were substantially different from the structures obtained from a distance-dependent dielectric simulation, with and without cutoff, and a constant dielectric model without a cutoff. In fact, structures from the distance-dependent dielectric simulations were insensitive to the nonbonded cutoff and in good agreement with the structures generated from the constant dielectric simulation without a cutoff. In addition, the solute-solvent temperature differential and solvent evaporation artifacts, characteristic of the constant dielectric simulation with a cutoff, were not present for the distance-dependent dielectric simulations. In this current work, we explore whether this dielectric-dependent cutoff-sensitive behavior for a constant dielectric model arises from the discontinuities in the forces at the nonbonded cutoff or from neglecting the structure-stabilizing interactions beyond the nonbonded cutoff. We also examine the origin of the dielectric-dependent artifacts, and its potential influence on the structural disparity. Several protocols for protein dynamics simulations are compared using both constant and distance-dependent dielectric models, including implementation of a switching function and a nonbonded cutoff and two different temperature coupling algorithms. We show that the distance-dependent dielectric model conserves energy in the SPASMS molecular mechanics and dynamics software for the time steps and nonbonded cutoffs commonly used in macromolecule simulations. Although the switching function simulation also conserved energy over a range of commonly used cutoffs, the constant dielectric model with a switching function yielded conformational results more similar to a constant dielectric simulation without a switching function than to a constant dielectric model without a nonbonded cutoff. Therefore, the conformational disparity between the dielectric models arises from neglecting important structure-stabilizing interactions beyond the cutoff, rather than differences in energy conservation. © 1993 John Wiley & Sons, Inc.     

Investigating detailed mechanism of hydrogen molecules adsorbing on single-wall carbon nanotubes using fitted force field parameters containing carbon–carbon interactions

The nonbonded and bonded force field parameters for carbon atoms in single-wall carbon nanotubes (SWNT) are fitted by means of quantum chemistry calculations with considering the periodic boundary conditions. The nonbonded parameters between carbon atoms and hydrogen atoms are fitted as well. All the fitted parameters are verified by comparing to quantum chemistry results and by calculating Young's modulus. Adsorption of Hydrogen molecules are then carried out on a bundle of self-assembled SWNTs. The adsorption isotherms are consistent to the Freundlich equation. Both hydrogen molecules adsorbed outside and inside the SWNTs are counted. According to our result, hydrogen molecules adsorbed inside the SWNTs are more stable at a relatively high temperature and are playing an important part in total amount of the adsorbed molecules. While C(10,10) have the highest adsorption capacities in most of the temperatures, hydrogen molecules inside C(5,5) are the most stable of all the four kinds of SWNTs. Thus, balancing adsorption capacities and strength of interaction can be important in choosing SWNT for gas adsorption. Besides, we deduct an equation that can describe the relation between hydrogen pressure and amount of SWNTs based on our simulation results. The hydrogen pressure may decrease by adding SWNTs in the system. The fitting method in our system is valid to SWNTs and can be tested in further studies of similar systems. © 2018 Wiley Periodicals, Inc.     

Nature of nonbonded Se...O interactions characterized by 17O NMR spectroscopy and NBO and AIM analyses

To investigate the nature of nonbonded Se...O interactions, three series of 2-substituted benzeneselenenyl derivatives [2-(CHO)C6H4SeX (1), 2-(CH2OH)C6H4SeX, (2), 2-(CH2OiPr)C6H4SeX (3); X = Cl, Br, CN, SPh, SeAr, Me] were synthesized. The 17O NMR absorption observed for 17O-enriched aldehydes 1 appeared upfield relative to benzaldehyde (PhCHO), while the opposite downfield shifts relative to benzyl alcohol (PhCH2OH) were observed for 17O-enriched alcohols 2 and ethers 3. The magnitude of both the upfield and the downfield shifts became larger as the electron-withdrawing ability of a substituent X increased. Quantum chemical calculations at the B3LYP level revealed that for all model compounds the most stable conformer has an intramolecular nonbonded Se.O interaction. Thus, the relative 17O NMR chemical shifts (DeltadeltaO) for 1-3 would reflect the strengths of the Se...O interactions. The natural bond orbital (NBO) analysis demonstrated that the stabilization energy due to an nO --> sigma Se-X orbital interaction (ESe...O) correlates with the Se...O atomic distance on a single curve irrespective of the type of the O atom. On the other hand, the atoms in molecules (AIM) analysis showed that the nonbonded Se...O interactions can be characterized by the presence of a bond critical point, the total energy density (HSe...O) of which decreases with strengthening of the interaction. The results suggested that Se...O interactions have a dominant covalent character rather than an electrostatic one.     

Parallelizing a molecular dynamics algorithm on a multiprocessor workstation using OpenMP

The atomistic molecular dynamics program YASP has been parallelized for shared-memory computer architectures. Parallelization was restricted to the most CPU-time-consuming parts: neighbor-list construction, calculation of nonbonded, angle and dihedral forces, and constraints. Most of the sequential FORTRAN code was kept; parallel constructs were inserted as compiler directives using the OpenMP standard. Only in the case of the neighbor list did the data structure have to be changed. The parallel code achieves a useful speedup over the sequential version for systems of several thousand atoms and above. On an IBM Regatta p690+, the throughput increases with the number of processors up to a maximum of 12-16 processors depending on the characteristics of the simulated systems. On dual-processor Xeon systems, the speedup is about 1.7.     

PVM-AMBER: A parallel implementation of the AMBER molecular mechanics package for workstation clusters

A parallel version of the popular molecular mechanics package AMBER suitable for execution on workstation clusters has been developed. Computer-intensive portions of molecular dynamics or free-energy perturbation computations, such as nonbonded pair list generation or calculation of nonbonded energies and forces, are distributed across a collection of Unix workstations linked by Ethernet or FDDI connections. This parallel implementation utilizes the message-passing software PVM (Parallel Virtual Machine) from Oak Ridge National Laboratory to coordinate data exchange and processor synchronization. Test simulations performed for solvated peptide, protein, and lipid bilayer systems indicate that reasonable parallel efficiency (70–90%) and computational speedup (2–5 × serial computer runtimes) can be achieved with small workstation clusters (typically six to eight machines) for typical biomolecular simulation problems. PVM-AMBER is also easily and rapidly portable to different hardware platforms due to the availability of PVM for numerous computers. The current version of PVM-AMBER has been tested successfully on Silicon Graphics, IBM RS6000, DEC ALPHA, and HP 735 workstation clusters and heterogeneous clusters of these machines, as well as on CRAY T3D and Kendall Square KSR2 parallel supercomputers. Thus, PVM-AMBER provides a simple and cost-effective mechanism for parallel molecular dynamics simulations on readily available hardware platforms. Factors that affect the efficiency of this approach are discussed. © 1995 by John Wiley & Sons, Inc.