Compressing molecular dynamics trajectories: Breaking the one‐bit‐per‐sample barrier

Molecular dynamics simulations yield large amounts of trajectory data. For their durable storage and accessibility an efficient compression algorithm is paramount. State of the art domain‐specific algorithms combine quantization, Huffman encoding and occasionally domain knowledge. We propose the high resolution trajectory compression scheme (HRTC) that relies on piecewise linear functions to approximate quantized trajectories. By splitting the error budget between quantization and approximation, our approach beats the current state of the art by several orders of magnitude given the same error tolerance. It allows storing samples at far less than one bit per sample. It is simple and fast enough to be integrated into the inner simulation loop, store every time step, and become the primary representation of trajectory data. © 2016 Wiley Periodicals, Inc.     

New ways to boost molecular dynamics simulations

We describe a set of algorithms that allow to simulate dihydrofolate reductase (DHFR, a common benchmark) with the AMBER all‐atom force field at 160 nanoseconds/day on a single Intel Core i7 5960X CPU (no graphics processing unit (GPU), 23,786 atoms, particle mesh Ewald (PME), 8.0 Å cutoff, correct atom masses, reproducible trajectory, CPU with 3.6 GHz, no turbo boost, 8 AVX registers). The new features include a mixed multiple time‐step algorithm (reaching 5 fs), a tuned version of LINCS to constrain bond angles, the fusion of pair list creation and force calculation, pressure coupling with a “densostat,” and exploitation of new CPU instruction sets like AVX2. The impact of Intel's new transactional memory, atomic instructions, and sloppy pair lists is also analyzed. The algorithms map well to GPUs and can automatically handle most Protein Data Bank (PDB) files including ligands. An implementation is available as part of the YASARA molecular modeling and simulation program from www.YASARA.org . © 2015 The Authors Journal of Computational Chemistry Published by Wiley Periodicals, Inc.     

An efficient and extensible format,library, and API for binary trajectory data from molecular simulations

Molecular dynamics simulations is an important application in theoretical chemistry, and with the large high‐performance computing resources available today the programs also generate huge amounts of output data. In particular in life sciences, with complex biomolecules such as proteins, simulation projects regularly deal with several terabytes of data. Apart from the need for more cost‐efficient storage, it is increasingly important to be able to archive data, secure the integrity against disk or file transfer errors, to provide rapid access, and facilitate exchange of data through open interfaces. There is already a whole range of different formats used, but few if any of them (including our previous ones) fulfill all these goals. To address these shortcomings, we present “Trajectory Next Generation” (TNG)—a flexible but highly optimized and efficient file format designed with interoperability in mind. TNG both provides state‐of‐the‐art multiframe compression as well as a container framework that will make it possible to extend it with new compression algorithms without modifications in programs using it. TNG will be the new file format in the next major release of the GROMACS package, but it has been implemented as a separate library and API with liberal licensing to enable wide adoption both in academic and commercial codes. © 2013 Wiley Periodicals, Inc.     

Energy extrapolation schemes for adaptive multi-scale molecular dynamics simulations

This paper evaluates simple schemes to extrapolate potential energy values using the set of energies and forces extracted from a molecular dynamics trajectory. In general, such a scheme affords the maximum amount of information about a molecular system at minimal computational cost. More specifically, schemes like this are very important in the field of adaptive multi-scale molecular dynamics simulations. In this field, often the computation of potential energy values at certain trajectory points is not required for the simulation itself, but solely for the a posteriori analysis of the simulation data. Extrapolating the values at these points from the available data can save considerable computational time. A set of extrapolation schemes are employed based on Taylor series and central finite difference approximations. The schemes are first tested on the trajectories of molecular systems of varying sizes, obtained at MM and QM level using velocity-Verlet integration with standard simulation time steps. Remarkably good accuracy was obtained with some of the approximations, while the failure of others can be explained in terms of the distinct features of a molecular dynamics trajectory. We have found that, for a Taylor expansion of the potential energy, both a first and a second order truncation exhibit errors that grow with system size. In contrast, the second order central finite difference approximation displays an accuracy that is independent of the size of the system, while giving a very good estimate of the energy, and costing as little as a first order truncation of the Taylor series. A fourth order central finite difference approximation requires more input data, which is not always available in adaptive multi-scale simulations. Furthermore, this approximation gives errors of similar magnitude or larger than its second order counterpart, at standard simulation time steps. This leads to the conclusion that a second order central finite difference approximation is the optimal choice for energy extrapolation from molecular dynamics trajectories. This finding is confirmed in a final application to the analysis of an adaptive multi-scale simulation.     

CHARMM‐GUI ligand reader and modeler for CHARMM force field generation of small molecules

Reading ligand structures into any simulation program is often nontrivial and time consuming, especially when the force field parameters and/or structure files of the corresponding molecules are not available. To address this problem, we have developed Ligand Reader & Modeler in CHARMM‐GUI. Users can upload ligand structure information in various forms (using PDB ID, ligand ID, SMILES, MOL/MOL2/SDF file, or PDB/mmCIF file), and the uploaded structure is displayed on a sketchpad for verification and further modification. Based on the displayed structure, Ligand Reader & Modeler generates the ligand force field parameters and necessary structure files by searching for the ligand in the CHARMM force field library or using the CHARMM general force field (CGenFF). In addition, users can define chemical substitution sites and draw substituents in each site on the sketchpad to generate a set of combinatorial structure files and corresponding force field parameters for throughput or alchemical free energy simulations. Finally, the output from Ligand Reader & Modeler can be used in other CHARMM‐GUI modules to build a protein‐ligand simulation system for all supported simulation programs, such as CHARMM, NAMD, GROMACS, AMBER, GENESIS, LAMMPS, Desmond, OpenMM, and CHARMM/OpenMM. Ligand Reader & Modeler is available as a functional module of CHARMM‐GUI at http://www.charmm-gui.org/input/ligandrm . © 2017 Wiley Periodicals, Inc.     

New method for approximate Hartree-Fock calculations using density approximations and coulomb field corrections. II

The HF approximation method that was outlined in Paper I is tested with respect to several molecular properties. Three different levels of approximation a, b, and c are considered. Satisfactory results—compared to corresponding “exact” HF calculations—are obtained with the STO -3G basis and the approximation level a. At this level the error in the binding energy is 0.001–0.025 a.u. for all considered molecules which contain up to six first-row atoms as, e.g., cyclopentanone (C₅OH₈). The error in the reaction energies considered here is about 4 kcal/mol (the maximal error is 9 kcal/mol). Orbital energies, dipole moments, gross charges, equilibrium geometries, and barriers to internal rotation are well reproduced by the approximation method at all three levels.     

Development of an internal searching algorithm for parameterization of the MM2/MM3 force fields

Application of Allinger's MM2/MM3 force fields to molecules of real interest is frequently hindered by the lack of parameters for various heterocyclic systems and for poly-functionalized molecules. A common approach to this problem is to manually choose missing parameters “by analogy” with those that are part of the force field's internal parameter set. Naturally, this is generally attempted only by those possessing extensive experience with force fields. In order to use the MM2/MM3 force fields to study herbicides, an algorithm has been developed to automate this process for the non MM2 specialist. Using a set of “relative cost” criteria for atom type replacement, the algorithm searches the force field parameter set and selects the most appropriate parameters for a given molecule whose MM2 output file contains “missing parameter” errors. The program selects parameter error messages from a standard MM2 output file, finds analogous parameters, asks the user to verify their appropriateness and creates a standard MM2 parameter deck for the molecule of interest.     

Chemical Syntheses and Chemical Biology of Carboxyl Polyether Ionophores: Recent Highlights

A central goal of chemical biology is to develop molecular probes that enable fundamental studies of cellular systems. In the hierarchy of bioactive molecules, the so‐called ionophore class occupies an unflattering position in the lower branches, with typical labels being “non‐specific” and “toxic”. In fact, the mere possibility that a candidate molecule possesses “ionophore activity” typically prompts its removal from further studies; ionophores—from a chemical genetics perspective—are molecular outlaws. In stark contrast to this overall poor reputation of ionophores, synthetic chemistry owes some of its most amazing achievements to studies of ionophore natural products, in particular the carboxyl polyethers renowned for their intricate molecular structures. These compounds have for decades been academic battlegrounds where new synthetic methodology is tested and retrosynthetic tactics perfected. Herein, we review the most exciting recent advances in carboxyl polyether ionophore (CPI) synthesis and in addition discuss the burgeoning field of CPI chemical biology.     

Response function theory of electron correlation

The full perturbation expansion for the response (or density—density correlation) function is examined in order to provide a useful general theory of excitation energies, oscillator strengths, dynamic polarizabilities, etc., that is more accurate than the random phase approximation. It is first shown how the formal partition of the diagrammatic version of the perturbation expansion into reducible and irreducible diagrams is generally useless as the latter category contains all the difficult terms which have heretofore resisted analysis in all but a haphazard form. It is then shown how the diagram for the response function can be partitioned into “correlated” and “uncorrelated” subsets. Restricting attention to the particle—hole blocks of the full response function, the “uncorrelated” diagrams desecribe the propagation of a particle—hole pair in an N-electron system where the particle and hole are each interacting with the remaining electrons but they are not interacting with each other. The “correlated” diagrams are those containing the hole—particle interactions, and, by defining a new class of reducible and irreducible diagrams, these are all summed to provide a perturbation expansion of the effective two-body hole—particle interaction that appears in the inverse of the response function. The “uncorrelated” diagrams are further partitioned into two sets, one of which is summed to all orders, while the other set is inverted in an order by order fashion. The final result presents a perturbation expansion for the inverse of the response function that is analogous to the Dyson equation for one-electron Green functions. Maintaining the perturbation expansion through first order for the inverse of the response function yields the eigenvalue equation of the familiar random phase approximation, while truncation at second order provides the most advanced theories that have been generated by the equations-of-motion method.     

Integral data compression for FPS 64-BIT processors: Improved I/O and reduced storage

Compression of two-electron integral data is used to reduce integral storage and I/O requirements with FPS M64 Series (formerly FPS-X64) processors. Schemes are developed and implemented in assembly language to compress floating-point values to a fixed-point accuracy, and unsigned integer numbers. The floating-point scheme stores only the significant bits of the mantissa and a short, biased exponent. The unsigned integers are packed into fixed-length fields just long enough to hold the largest value. The packing procedures are tested on FPS-164 and FPS-264 processors (since renamed M64/145 and M64/60 by FPS) and incorporated into HONDO to compress two-electron integral files. Reduction factors of 0.2–0.4 are obtained for floating-point compression and 0.3–0.5 for index packing, with typical overall factors around one-third. The advantages of improved I/O and storage efficiency are accompanied by a small increase in processor time to perform the packing and unpacking. Timing changes for HONDO are presented, and both packing schemes dramatically reduce SCF elapsed times with FPS-264 processors. It is concluded that compression effectively extends external storage capacities, improves I/O efficiency, and can reduce the elapsed time of I/O bound calculations.     

Molecular dynamics of 300 keV (D2O)100 clusters on Ti-D

We use the molecular dynamics code DAMSEL to predict the velocity distributions for beam and lattice atoms after bombardment of Ti-D “foils” of thickness 20.86 Å by 300 keV (D₂O)₁₀₀ cluster ions. From these distributions we estimate the D-D nuclear fusion yield. We find that cluster bombardment reduces the overall energy deposition of the beam in the lattice compared to that of the individual beam atoms of the same velocity. However, a small portion of lattice atoms (<1%), and a larger percentage of beam atoms (～30%), have energies above the maximum present in the case of bombardment by individual D or O atoms. The folding of the standard D-D fusion cross sections over the relative velocity distributions produced by beam and lattice deuterons produce a fusion yield estimate of ～1×10^?21 fusions per cluster, with the high-energy distributions of beam deuterons playing the most important role. This is nine orders of magnitude lower that the data of Beuhler et al. While transient (～10 fs) atom densities 50% higher than that of the initial lattice are recorded in the course of the simulation, the average energy transferred per lattice atom — 23 eV — is insufficient to support any “collision spike” explanation of the observed fusion yield.     

Differential Behavior of Amino–Imino Constitutional Isomers in Nonlinear Optical Processes

A detailed study of the “blocked” amino–imino tautomers derived from N‐acridine‐substituted 2‐aminobenzothiazole—and their effect on the nonlinear optical response—is presented. The synthesis, characterization, and nonlinear optical properties of these frozen tautomers, namely, N‐methyl‐N‐(2‐nitroacridin‐6‐yl)‐2‐aminobenzothia‐zole and 3‐methyl‐N‐(7‐nitroacridin‐3‐yl)‐2‐iminobenzothiazole, are reported. A theoretical model based on valence–bond theory is also proposed and used to analyze the effects of the nuclear configuration corresponding to each frozen tautomer structure. In the present case, the aromatic form and the allylic‐anion‐like system of the ? N? C? N? group inherent to each isomer are crucial for understanding and analyzing the different responses of each “blocked” tautomer.     

Adaptative finite element simulation of currents at microelectrodes to a guaranteed accuracy. Application to a simple model problem

In this series of papers we consider the general problem of numerical simulation of the currents at microelectrodes using an adaptive finite element approach. Microelectrodes typically consist of an electrode embedded (or recessed) in an insulating material. For all such electrodes, numerical simulation is made difficult by the presence of a boundary singularity at the electrode edge (where the electrode meets the insulator), manifested by the large increase in the current density at this point, often referred to as the ‘edge-effect’. Our approach to overcoming this problem involves the derivation of an a posteriori bound on the error in the numerical approximation for the current that can be used to drive an adaptive mesh-generation algorithm. This allows us to calculate the current to within a prescribed tolerance. We begin by demonstrating the power of the method for a simple model problem — an E reaction mechanism at a microdisc electrode — for which the analytical solution is known. In this paper we give the background to the problem, and show how an a posteriori error bound can be used to drive an adaptive mesh-generation algorithm. We then use the algorithm to solve our model problem and obtain very accurate results on comparatively coarse meshes in minimal computing time. We give the technical details of the background theory and the derivation of the error bound in the accompanying paper.     

Comparison of radii sets,entropy, QM methods,and sampling on MM‐PBSA,MM‐GBSA,and QM/MM‐GBSA ligand binding energies of F. tularensis enoyl‐ACP reductase (FabI)

To validate a method for predicting the binding affinities of FabI inhibitors, three implicit solvent methods, MM‐PBSA, MM‐GBSA, and QM/MM‐GBSA were carefully compared using 16 benzimidazole inhibitors in complex with Francisella tularensis FabI. The data suggests that the prediction results are sensitive to radii sets, GB methods, QM Hamiltonians, sampling protocols, and simulation length, if only one simulation trajectory is used for each ligand. In this case, QM/MM‐GBSA using 6 ns MD simulation trajectories together with GB^neck2, PM3, and the mbondi2 radii set, generate the closest agreement with experimental values (r² = 0.88). However, if the three implicit solvent methods are averaged from six 1 ns MD simulations for each ligand (called “multiple independent sampling”), the prediction results are relatively insensitive to all the tested parameters. Moreover, MM/GBSA together with GB^HCT and mbondi, using 600 frames extracted evenly from six 0.25 ns MD simulations, can also provide accurate prediction to experimental values (r² = 0.84). Therefore, the multiple independent sampling method can be more efficient than a single, long simulation method. Since future scaffold expansions may significantly change the benzimidazole's physiochemical properties (charges, etc.) and possibly binding modes, which may affect the sensitivities of various parameters, the relatively insensitive “multiple independent sampling method” may avoid the need of an entirely new validation study. Moreover, due to large fluctuating entropy values, (QM/)MM‐P(G)BSA were limited to inhibitors’ relative affinity prediction, but not the absolute affinity. The developed protocol will support an ongoing benzimidazole lead optimization program. © 2015 Wiley Periodicals, Inc.     

Chiral Atropisomeric Indenocorannulene Bowls: Critique of the Cahn–Ingold–Prelog Conception of Molecular Chirality

Chiral corannulenes abound, but suffer generally from configurational lability associated with bowl‐to‐bowl inversion, 1 thus obviating questions of stereogenicity and stereoelement construction. 2 In contrast, peri‐annulated corannulenes show greatly increased barriers for bowl‐to‐bowl inversion; specifically indenocorannulenes invert on a time scale too slow to observe by normal NMR methods and raise the possibility of creating chiral atropisomeric bowl‐shaped aromatics. 3 Two methods for preparing indenocorannulene from simple 2‐haloarylcorannulenes—silyl cation C–F activation, 4 and Pd‐mediated C–Cl activation^[5]—enable the synthesis of an array of such chiral atropisomeric indenocorannulenes. 6 Resolution of the enantiomers by high‐performance liquid chromatography over chiral support phases motivates the study of chiroptical properties, the assignment of absolute “Cartesian” configuration, and the assessment of configurational stability. 7 These studies bring into question any systematic assignment of nontrivial stereoelements (i.e. not the molecule in its entirety) and refute any assertion of congruence between “Cahn–Ingold–Prelog elements” and the physical or “Cartesian” basis of chirality.     

Rapid sampling of all‐atom peptides using a library‐based polymer‐growth approach

We adapted existing polymer growth strategies for equilibrium sampling of peptides described by modern atomistic forcefields with a simple uniform dielectric solvent. The main novel feature of our approach is the use of precalculated statistical libraries of molecular fragments. A molecule is sampled by combining fragment configurations—of single residues in this study—which are stored in the libraries. Ensembles generated from the independent libraries are reweighted to conform with the Boltzmann‐factor distribution of the forcefield describing the full molecule. In this way, high‐quality equilibrium sampling of small peptides (4–8 residues) typically requires less than one hour of single‐processor wallclock time and can be significantly faster than Langevin simulations. Furthermore, approximate, clash‐free ensembles can be generated for larger peptides (up to 32 residues in this study) in less than a minute of single‐processor computing. We discuss possible applications of our growth procedure to free energy calculation, fragment assembly protein‐structure prediction protocols, and to “multi‐resolution” sampling. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2011