Accelerating chemical database searching using graphics processing units

The utility of chemoinformatics systems depends on the accurate computer representation and efficient manipulation of chemical compounds. In such systems, a small molecule is often digitized as a large fingerprint vector, where each element indicates the presence/absence or the number of occurrences of a particular structural feature. Since in theory the number of unique features can be exceedingly large, these fingerprint vectors are usually folded into much shorter ones using hashing and modulo operations, allowing fast "in-memory" manipulation and comparison of molecules. There is increasing evidence that lossless fingerprints can substantially improve retrieval performance in chemical database searching (substructure or similarity), which have led to the development of several lossless fingerprint compression algorithms. However, any gains in storage and retrieval afforded by compression need to be weighed against the extra computational burden required for decompression before these fingerprints can be compared. Here we demonstrate that graphics processing units (GPU) can greatly alleviate this problem, enabling the practical application of lossless fingerprints on large databases. More specifically, we show that, with the help of a ~$500 ordinary video card, the entire PubChem database of ~32 million compounds can be searched in ~0.2-2 s on average, which is 2 orders of magnitude faster than a conventional CPU. If multiple query patterns are processed in batch, the speedup is even more dramatic (less than 0.02-0.2 s/query for 1000 queries). In the present study, we use the Elias gamma compression algorithm, which results in a compression ratio as high as 0.097.     

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.     

Hybrid Monte Carlo: An efficient algorithm for condensed matter simulation

A detailed comparison has been made of the performance of molecular dynamics and hybrid Monte Carlo simulation algorithms for calculating thermodynamic properties of 2D Lennard-Jonesium. The hybrid Monte Carlo simulation required an order of magnitude fewer steps than the molecular dynamics simulation to calculate reproducible values of the specific heat. The ergodicity of the two algorithms was compared via the use of intermediate scattering functions. For classical systems the intermediate scattering functions should be real; however, a simple analysis demonstrates that this function will have a significant imaginary component when ergodicity breaks down. For q vectors near the zone boundary, the scattering functions are real for both algorithms. However, for q vectors near the zone center (i.e., harmonic, weakly coupled modes), the scattering function calculated via molecular dynamics had a significantly larger imaginary component than that calculated using hybrid Monte Carlo. Therefore, the hybrid Monte Carlo algorithm is more ergodic and samples phase space more efficiently than molecular dynamics for simulations of 2D Lennard-Jonesium. © 1994 by John Wiley & Sons, Inc.     

Efficient compression of molecular dynamics trajectory files

We investigate whether specific properties of molecular dynamics trajectory files can be exploited to achieve effective file compression. We explore two classes of lossy, quantized compression scheme: “interframe” predictors, which exploit temporal coherence between successive frames in a simulation, and more complex “intraframe” schemes, which compress each frame independently. Our interframe predictors are fast, memory‐efficient and well suited to on‐the‐fly compression of massive simulation data sets, and significantly outperform the benchmark BZip2 application. Our schemes are configurable: atomic positional accuracy can be sacrificed to achieve greater compression. For high fidelity compression, our linear interframe predictor gives the best results at very little computational cost: at moderate levels of approximation (12‐bit quantization, maximum error ≈ 10^?2 Å), we can compress a 1–2 fs trajectory file to 5–8% of its original size. For 200 fs time steps—typically used in fine grained water diffusion experiments—we can compress files to ～25% of their input size, still substantially better than BZip2. While compression performance degrades with high levels of quantization, the simulation error is typically much greater than the associated approximation error in such cases. © 2012 Wiley Periodicals, Inc.     

A new extension of classical molecular dynamics: An electron transfer algorithm

The molecular dynamics is one of the most widely used methods for the simulation of the properties corresponding to ionic motion. Unfortunately, classical molecular dynamics cannot be applied for electron transfer simulation. Suggested modification of the molecular dynamics allows performing the electron transfer from one particle to another during simulation runtime. All additional data structure and the corresponding algorithms are presented in this article. The method can be applied to the systems with pair Van der Waals and Coulomb interactions. Moreover, it may be extended for many‐bodied interatomic interactions. In addition, an algorithm of transference numbers calculation has been designed. This extension is not an independent method but it can be useful for simulating the systems with high concentration of electron donors and acceptors. © 2017 Wiley Periodicals, Inc.     

Routine Microsecond Molecular Dynamics Simulations with AMBER on GPUs. 1. Generalized Born

We present an implementation of generalized Born implicit solvent all-atom classical molecular dynamics (MD) within the AMBER program package that runs entirely on CUDA enabled NVIDIA graphics processing units (GPUs). We discuss the algorithms that are used to exploit the processing power of the GPUs and show the performance that can be achieved in comparison to simulations on conventional CPU clusters. The implementation supports three different precision models in which the contributions to the forces are calculated in single precision floating point arithmetic but accumulated in double precision (SPDP), or everything is computed in single precision (SPSP) or double precision (DPDP). In addition to performance, we have focused on understanding the implications of the different precision models on the outcome of implicit solvent MD simulations. We show results for a range of tests including the accuracy of single point force evaluations and energy conservation as well as structural properties pertainining to protein dynamics. The numerical noise due to rounding errors within the SPSP precision model is sufficiently large to lead to an accumulation of errors which can result in unphysical trajectories for long time scale simulations. We recommend the use of the mixed-precision SPDP model since the numerical results obtained are comparable with those of the full double precision DPDP model and the reference double precision CPU implementation but at significantly reduced computational cost. Our implementation provides performance for GB simulations on a single desktop that is on par with, and in some cases exceeds, that of traditional supercomputers.     

Parallel replica dynamics with a heterogeneous distribution of barriers: application to n-hexadecane pyrolysis

Parallel replica dynamics simulation methods appropriate for the simulation of chemical reactions in molecular systems with many conformational degrees of freedom have been developed and applied to study the microsecond-scale pyrolysis of n-hexadecane in the temperature range of 2100-2500 K. The algorithm uses a transition detection scheme that is based on molecular topology, rather than energetic basins. This algorithm allows efficient parallelization of small systems even when using more processors than particles (in contrast to more traditional parallelization algorithms), and even when there are frequent conformational transitions (in contrast to previous implementations of the parallel replica algorithm). The parallel efficiency for pyrolysis initiation reactions was over 90% on 61 processors for this 50-atom system. The parallel replica dynamics technique results in reaction probabilities that are statistically indistinguishable from those obtained from direct molecular dynamics, under conditions where both are feasible, but allows simulations at temperatures as much as 1000 K lower than direct molecular dynamics simulations. The rate of initiation displayed Arrhenius behavior over the entire temperature range, with an activation energy and frequency factor of E(a) = 79.7 kcal/mol and log A/s(-1) = 14.8, respectively, in reasonable agreement with experiment and empirical kinetic models. Several interesting unimolecular reaction mechanisms were observed in simulations of the chain propagation reactions above 2000 K, which are not included in most coarse-grained kinetic models. More studies are needed in order to determine whether these mechanisms are experimentally relevant, or specific to the potential energy surface used.     

Reactive Bond-Order Simulations Using Both Spatial and Temporal Approaches to Parallelism

We describe two different approaches to exploiting parallel computing architecture that have been used successfully for reactive molecular simulation using bond-order potentials. These potentials are based on the Tersoff bond-order formalism, and allow accurate treatement of covalent bonding reactions in the framework of a classical potential. They include both Brenner's reactive empirical bond order (REBO) potential and our adaptive intermolecular version of this potential (AIREBO). Traditional spatial and atom-based parallel decompositioon techniques have been employed in the RMD-CE program developed for parallel molecular dynamics simulations with a variety of reactive potentials. Key features of this implementation, including the object-oriented approach and novel algorithms for the integrator and neighbor lists, are discussed. The resulting code provides efficient scaling down to system sizes of 400 atoms per processor, and has been applied successfully to systems of as many as half a million atoms. For smaller systems, the parallel replica dynamics algorithm has been successfully applied to take advantage of parallelism in the time domain for rare-event systems. This approach takes advantage of the independence of different parts of a dynamics trajectory, and provides excellent parallel efficiencies for systems as small as tens of atoms, where other parallel simulation techniques are not applicable. This technique has been used to model the pyrolysis of hexadecane on the microsecond timescale, at more realistic temperatures than are achievable with other simulation methods.     

Lossless compression of chemical fingerprints using integer entropy codes improves storage and retrieval

Many modern chemoinformatics systems for small molecules rely on large fingerprint vector representations, where the components of the vector record the presence or number of occurrences in the molecular graphs of particular combinatorial features, such as labeled paths or labeled trees. These large fingerprint vectors are often compressed to much shorter fingerprint vectors using a lossy compression scheme based on a simple modulo procedure. Here, we combine statistical models of fingerprints with integer entropy codes, such as Golomb and Elias codes, to encode the indices or the run lengths of the fingerprints. After reordering the fingerprint components by decreasing frequency order, the indices are monotone-increasing and the run lengths are quasi-monotone-increasing, and both exhibit power-law distribution trends. We take advantage of these statistical properties to derive new efficient, lossless, compression algorithms for monotone integer sequences: monotone value (MOV) coding and monotone length (MOL) coding. In contrast to lossy systems that use 1024 or more bits of storage per molecule, we can achieve lossless compression of long chemical fingerprints based on circular substructures in slightly over 300 bits per molecule, close to the Shannon entropy limit, using a MOL Elias Gamma code for run lengths. The improvement in storage comes at a modest computational cost. Furthermore, because the compression is lossless, uncompressed similarity (e.g., Tanimoto) between molecules can be computed exactly from their compressed representations, leading to significant improvements in retrival performance, as shown on six benchmark data sets of druglike molecules.     

Exploring the role of decoherence in condensed-phase nonadiabatic dynamics: a comparison of different mixed quantum/classical simulation algorithms for the excited hydrated electron

Mixed quantum/classical (MQC) molecular dynamics simulation has become the method of choice for simulating the dynamics of quantum mechanical objects that interact with condensed-phase systems. There are many MQC algorithms available, however, and in cases where nonadiabatic coupling is important, different algorithms may lead to different results. Thus, it has been difficult to reach definitive conclusions about relaxation dynamics using nonadiabatic MQC methods because one is never certain whether any given algorithm includes enough of the necessary physics. In this paper, we explore the physics underlying different nonadiabatic MQC algorithms by comparing and contrasting the excited-state relaxation dynamics of the prototypical condensed-phase MQC system, the hydrated electron, calculated using different algorithms, including: fewest-switches surface hopping, stationary-phase surface hopping, and mean-field dynamics with surface hopping. We also describe in detail how a new nonadiabatic algorithm, mean-field dynamics with stochastic decoherence (MF-SD), is to be implemented for condensed-phase problems, and we apply MF-SD to the excited-state relaxation of the hydrated electron. Our discussion emphasizes the different ways quantum decoherence is treated in each algorithm and the resulting implications for hydrated-electron relaxation dynamics. We find that for three MQC methods that use Tully's fewest-switches criterion to determine surface hopping probabilities, the excited-state lifetime of the electron is the same. Moreover, the nonequilibrium solvent response function of the excited hydrated electron is the same with all of the nonadiabatic MQC algorithms discussed here, so that all of the algorithms would produce similar agreement with experiment. Despite the identical solvent response predicted by each MQC algorithm, we find that MF-SD allows much more mixing of multiple basis states into the quantum wave function than do other methods. This leads to an excited-state lifetime that is longer with MF-SD than with any method that incorporates nonadiabatic effects with the fewest-switches surface hopping criterion.     

Nonadiabatic excited-state molecular dynamics: numerical tests of convergence and parameters

Nonadiabatic molecular dynamics simulations, involving multiple Born-Oppenheimer potential energy surfaces, often require a large number of independent trajectories in order to achieve the desired convergence of the results, and simulation relies on different parameters that should be tested and compared. In addition to influencing the speed of the simulation, the chosen parameters combined with the frequently reduced number of trajectories can sometimes lead to unanticipated changes in the accuracy of the simulated dynamics. We have previously developed a nonadiabatic excited state molecular dynamics methodology employing Tully's fewest switches surface hopping algorithm. In this study, we seek to investigate the impact of the number of trajectories and the various parameters on the simulation of the photoinduced dynamics of distyrylbenzene (a small oligomer of polyphenylene vinylene) within our developed framework. Various user-defined parameters are analyzed: classical and quantum integration time steps, the value of the friction coefficient for Langevin dynamics, and the initial seed used for stochastic thermostat and hopping algorithms. Common approximations such as reduced number of nonadiabatic coupling terms and the classical path approximation are also investigated. Our analysis shows that, at least for the considered molecular system, a minimum of ~400 independent trajectories should be calculated in order to achieve statistical averaging necessary for convergence of the calculated relaxation timescales.     

Molecular dynamics simulation of configurational ensembles compatible with experimental FRET efficiency data through a restraint on instantaneous FRET efficiencies

Förster resonance energy transfer (FRET) measurements are widely used to investigate (bio)molecular interactions or/and association. FRET efficiencies, the primary data obtained from this method, give, in combination with the common assumption of isotropic chromophore orientation, detailed insight into the lengthscale of molecular phenomena. This study illustrates the application of a FRET efficiency restraint during classical atomistic molecular dynamics simulations of a mutant mastoparan X peptide in either water or 7 M aqueous urea. The restraint forces acting on the donor and acceptor chromophores ensure that the sampled peptide configurational ensemble satisfies the experimental primary data by modifying interchromophore separation and chromophore transition dipole moment orientations. By means of a conformational cluster analysis, it is seen that indeed different configurational ensembles may be sampled without and with application of the restraint. In particular, while the FRET efficiency and interchromophore distances monitored in an unrestrained simulation may differ from the experimentally‐determined values, they can be brought in agreement with experimental data through usage of the FRET efficiency restraining potential. Furthermore, the present results suggest that the assumption of isotropic chromophore orientation is not always justified. The FRET efficiency restraint allows the generation of configurational ensembles that may not be accessible with unrestrained simulations, and thereby supports a meaningful interpretation of experimental FRET results in terms of the underlying molecular degrees of freedom. Thus, it offers an additional tool to connect the realms of computer and wet‐lab experimentation. © 2014 Wiley Periodicals, Inc.     

Ab initio based force field and molecular dynamics simulations of crystalline TATB

An all-atom force field for 1,3,5-triamino-2,4,6-trinitrobenzene (TATB) is presented. The classical intermolecular interaction potential for TATB is based on single-point energies determined from high-level ab initio calculations of TATB dimers. The newly developed potential function is used to examine bulk crystalline TATB via molecular dynamics simulations. The isobaric thermal expansion and isothermal compression under hydrostatic pressures obtained from the molecular dynamics simulations are in good agreement with experiment. The calculated volume-temperature expansion is almost one dimensional along the c crystallographic axis, whereas under compression, all three unit cell axes participate, albeit unequally.     

Chemical reaction optimization for solving shortest common supersequence problem

Shortest common supersequence (SCS) is a classical NP-hard problem, where a string to be constructed that is the supersequence of a given string set. The SCS problem has an enormous application of data compression, query optimization in the database and different bioinformatics activities. Due to NP-hardness, the exact algorithms fail to compute SCS for larger instances. Many heuristics and meta-heuristics approaches were proposed to solve this problem. In this paper, we propose a meta-heuristics approach based on chemical reaction optimization, CRO_SCS that is designed inspired by the nature of the chemical reactions. For different optimization problems like 0-1 knapsack, quadratic assignment, global numeric optimization problems CRO algorithm shows very good performance. We have redesigned the reaction operators and a new reform function to solve the SCS problem. The outcomes of the proposed CRO_SCS algorithm are compared with those of the enhanced beam search (IBS_SCS), deposition and reduction (DR), ant colony optimization (ACO) and artificial bee colony (ABC) algorithms. The length of supersequence, execution time and standard deviation of all related algorithms show that CRO_SCS gives better results on the average than all other algorithms.     

Image analysis techniques for automatic evaluation of two-dimensional electrophoresis.

Techniques for automatic analysis of two-dimensional electrophoresis gels by computer-aided image analysis are described. Original gels or photographic films are scanned using a laser scanner and the files are transferred to a microcomputer. The program package first performs a compression and preevaluation of the files. Spot identification and quantification is performed by the chain code algorithm after appropriate zooming and cutting. Labeling facilitates spot identification and quantification in numerical and graphical (pseudocolor) representation on peripheral devices for camera ready output. Interpolation between measured basepoints is performed by cubic spline algorithms which are automatically switched on and off, depending on the need by the program. High speed analysis and graphic representation is achieved using fast Assembler language routines rather than high level languages. One-dimensional gels can be analyzed using the same software. Spot matching between parallel two-dimensional gels has not yet been implemented.     

Time-reversible ab initio molecular dynamics

Time-reversible ab initio molecular dynamics based on a lossless multichannel decomposition for the integration of the electronic degrees of freedom [Phys. Rev. Lett. 97, 123001 (2006)] is explored. The authors present a lossless time-reversible density matrix molecular dynamics scheme. This approach often allows for stable Hartree-Fock simulations using only one single self-consistent field cycle per time step. They also present a generalization, introducing an additional "forcing" term, that in a special case includes a hybrid Lagrangian, i.e., Car-Parrinello-type, method, which can systematically be constrained to the Born-Oppenheimer potential energy surface by using an increasing number of self-consistency cycles in the nuclear force calculations. Furthermore, in analog to the reversible and symplectic leapfrog or velocity Verlet schemes, where not only the position but also the velocity is propagated, the authors propose a Verlet-type density velocity formalism for time-reversible Born-Oppenheimer molecular dynamics.     

Modeling signal transduction networks: a comparison of two stochastic kinetic simulation algorithms

Computational efficiency of stochastic kinetic algorithms depend on factors such as the overall species population, the total number of reactions, and the average number of nodal interactions or connectivity in a network. These size measures of the network model can have a significant impact on computational efficiency. In this study, two scalable biological networks are used to compare the size scaling efficiencies of two popular and conceptually distinct stochastic kinetic simulation algorithms--the random substrate method of Firth and Bray (FB), and the Gillespie algorithm as implemented using the Gibson-Bruck method (GGB). The arithmetic computational efficiencies of these two algorithms, respectively, scale with the square of the total species population and the logarithm of the total number of active reactions. The two scalable models considered are the size scalable model (SSM), a four compartment reaction model for a signal transduction network involving receptors with single phosphorylation binding sites, and the variable connectivity model (VCM), a single compartment model where receptors possess multiple phosphorylation binding sites. The SSM has fixed species connectivity while the connectivity between species in VCM increases with the number of phosphorylation sites. For SSM, we find that, as the total species population is increased over four orders of magnitude, the GGB algorithm performs significantly better than FB for all three SSM compartment models considered. In contrast, for VCM, we find that as the overall species population decreases while the number of phosphorylation sites increases (implying an increase in network linkage) there exists a crossover point where the computational demands of the GGB method exceed that of the FB.     

Rigid-body dynamics in the isothermal-isobaric ensemble: a test on the accuracy and computational efficiency

We have developed a time-reversible rigid-body (rRB) molecular dynamics algorithm in the isothermal-isobaric (NPT) ensemble. The algorithm is an extension of rigid-body dynamics [Matubayasi and Nakahara, J Chem Phys 1999, 110, 3291] to the NPT ensemble on the basis of non-Hamiltonian statistical mechanics [Martyna, G. J. et al., J Chem Phys 1994, 101, 4177]. A series of MD simulations of water as well as fully hydrated lipid bilayer systems have been undertaken to investigate the accuracy and efficiency of the algorithm. The rRB algorithm was shown to be superior to the state-of-the-art constraint-dynamics algorithm SHAKE/RATTLE/ROLL, with respect to computational efficiency. However, it was revealed that both algorithms produced accurate trajectories of molecules in the NPT as well as NVT ensembles, as long as a reasonably short time step was used. A couple of multiple time-step (MTS) integration schemes were also examined. The advantage of the rRB algorithm for computational efficiency increased when the MD simulation was carried out using MTS on parallel processing computer systems; total computer time for MTS-MD of a lipid bilayer using 64 processors was reduced by about 40% using rRB instead of SHAKE/RATTLE/ROLL.     

Rapid grid-based construction of the molecular surface and the use of induced surface charge to calculate reaction field energies: applications to the molecular systems and geometric objects

This article describes a number of algorithms that are designed to improve both the efficiency and accuracy of finite difference solutions to the Poisson-Boltzmann equation (the FDPB method) and to extend its range of application. The algorithms are incorporated in the DelPhi program. The first algorithm involves an efficient and accurate semianalytical method to map the molecular surface of a molecule onto a three-dimensional lattice. This method constitutes a significant improvement over existing methods in terms of its combination of speed and accuracy. The DelPhi program has also been expanded to allow the definition of geometrical objects such as spheres, cylinders, cones, and parallelepipeds, which can be used to describe a system that may also include a standard atomic level depiction of molecules. Each object can have a different dielectric constant and a different surface or volume charge distribution. The improved definition of the surface leads to increased precision in the numerical solutions of the PB equation that are obtained. A further improvement in the precision of solvation energy calculations is obtained from a procedure that calculates induced surface charges from the FDPB solutions and then uses these charges in the calculation of reaction field energies. The program allows for finite difference grids of large dimension; currently a maximum of 571(3) can be used on molecules containing several thousand atoms and charges. As described elsewhere, DelPhi can also treat mixed salt systems containing mono- and divalent ions and provide electrostatic free energies as defined by the nonlinear PB equation.     

A novel reversible pH-triggered release immobilized enzyme system

A novel immobilized enzyme system supported by poly(acrylic acid/N,N’-methylene-bisacryl-amide) hydrogel microspheres was prepared. This system exhibited characteristics of reversible pH-triggered release. The morphology, size, and chemical structure were examined through optical microscopy, particle size analyzer, and Fourier transform infrared spectrometer. Immobilization and release features were further investigated under different conditions, including pH, time, and microsphere quantity. Results showed the microspheres were regularly spherical with 3.8 ~ 6.6 μm diameter. Loading efficiencies of bovine serum albumin immobilized by gel entrapment and adsorption methods were 93.9% and 56.2%, respectively. The pH-triggered protein release of the system occurred when medium pH was above 6.0, while it was hardly detected when medium pH was below 6.0. Release efficiencies of entrapped and adsorbed protein were 6.38% and 95.0%, respectively. Hence, adsorption method was used to immobilize trypsin. Loading efficiency of 77.2% was achieved at pH 4.0 in 1 h. Release efficiency of 91.6% was obtained under optimum pH catalysis condition set at 8.0 and trypsin was free in solutions with retention activity of 63.3%. And 51.5% of released trypsin could be reloaded in 10 min. The results indicate this kind of immobilized enzyme system offers a promising alternative for enzyme recovery in biotechnology.