Three pillars for achieving quantum mechanical molecular dynamics simulations of huge systems: Divide‐and‐conquer,density‐functional tight‐binding,and massively parallel computation

The linear‐scaling divide‐and‐conquer (DC) quantum chemical methodology is applied to the density‐functional tight‐binding (DFTB) theory to develop a massively parallel program that achieves on‐the‐fly molecular reaction dynamics simulations of huge systems from scratch. The functions to perform large scale geometry optimization and molecular dynamics with DC‐DFTB potential energy surface are implemented to the program called DC‐DFTB‐K. A novel interpolation‐based algorithm is developed for parallelizing the determination of the Fermi level in the DC method. The performance of the DC‐DFTB‐K program is assessed using a laboratory computer and the K computer. Numerical tests show the high efficiency of the DC‐DFTB‐K program, a single‐point energy gradient calculation of a one‐million‐atom system is completed within 60 s using 7290 nodes of the K computer. © 2016 Wiley Periodicals, Inc.     

Performance of density‐functional tight‐binding models in describing hydrogen‐bonded anionic‐water clusters

Density‐functional tight‐binding (DFTB) models are computationally efficient approximations to density‐functional theory that have been shown to predict reliable structural and energetic properties for various systems. In this work, the reliability and accuracy of the self‐consistent‐charge DFTB model and its recent extension(s) in predicting the structures, binding energies, charge distributions, and vibrational frequencies of small water clusters containing polyatomic anions of the Hofmeister series (carbonate, sulfate, hydrogen phosphate, acetate, nitrate, perchlorate, and thiocyanate) have been carefully and systematically evaluated on the basis of high‐level ab initio quantum‐chemistry [MP2/aug‐cc‐pVTZ and CCSD(T)/aug‐cc‐pVQZ] reference data. Comparison with available experimental data has also been made for further validation. The self‐consistent‐charge DFTB model, and even more so its recent extensions, are shown to properly account for the structural properties, energetics, intermolecular polarization, and spectral signature of hydrogen‐bonding in anionic water clusters at a fraction of the computational cost of ab initio quantum‐chemistry methods. This makes DFTB models candidates of choice for investigating much larger systems such as seeded water droplets, their structural properties, formation thermodynamics, and infrared spectra. © 2014 Wiley Periodicals, Inc.     

Rapid QM/MM approach for biomolecular systems under periodic boundary conditions: Combination of the density‐functional tight‐binding theory and particle mesh Ewald method

A quantum mechanical/molecular mechanical (QM/MM) approach based on the density‐functional tight‐binding (DFTB) theory is a useful tool for analyzing chemical reaction systems in detail. In this study, an efficient QM/MM method is developed by the combination of the DFTB/MM and particle mesh Ewald (PME) methods. Because the Fock matrix, which is required in the DFTB calculation, is analytically obtained by the PME method, the Coulomb energy is accurately and rapidly computed. For assessing the performance of this method, DFTB/MM calculations and molecular dynamics simulation are conducted for a system consisting of two amyloid‐β(1‐16) peptides and a zinc ion in explicit water under periodic boundary conditions. As compared with that of the conventional Ewald summation method, the computational cost of the Coulomb energy by utilizing the present approach is drastically reduced, i.e., 166.5 times faster. Furthermore, the deviation of the electronic energy is less than . © 2016 Wiley Periodicals, Inc.     

Energy gradients in combined fragment molecular orbital and polarizable continuum model (FMO/PCM) calculation

The analytic energy gradients for the combined fragment molecular orbital and polarizable continuum model (FMO/PCM) method are derived and implemented. Applications of FMO/PCM geometry optimization to polyalanine show that the structures obtained with the FMO/PCM method are very close to those obtained with the corresponding full ab initio PCM methods. FMO/PCM (RHF/6‐31G* level) is used to optimize the solution structure of the 304‐atom Trp‐cage miniprotein and the result is in agreement with NMR experiments. The key factors determining the relative stability of the α‐helix, β‐turn and the extended form in solution are elucidated for polyalanine. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2010     

Parallel implementation of efficient charge–charge interaction evaluation scheme in periodic divide‐and‐conquer density‐functional tight‐binding calculations

A low‐computational‐cost algorithm and its parallel implementation for periodic divide‐and‐conquer density‐functional tight‐binding (DC‐DFTB) calculations are presented. The developed algorithm enables rapid computation of the interaction between atomic partial charges, which is the bottleneck for applications to large systems, by means of multipole‐ and interpolation‐based approaches for long‐ and short‐range contributions. The numerical errors of energy and forces with respect to the conventional Ewald‐based technique can be under the control of the multipole expansion order, level of unit cell replication, and interpolation grid size. The parallel performance of four different evaluation schemes combining previous approaches and the proposed one are assessed using test calculations of a cubic water box on the K computer. The largest benchmark system consisted of 3,295,500 atoms. DC‐DFTB energy and forces for this system were obtained in only a few minutes when the proposed algorithm was activated and parallelized over 16,000 nodes in the K computer. The high performance using a single node workstation was also confirmed. In addition to liquid water systems, the feasibility of the present method was examined by testing solid systems such as diamond form of carbon, face‐centered cubic form of copper, and rock salt form of sodium chloride. © 2017 Wiley Periodicals, Inc.     

Assessment and acceleration of binding energy calculations for protein–ligand complexes by the fragment molecular orbital method

In the field of drug discovery, it is important to accurately predict the binding affinities between target proteins and drug applicant molecules. Many of the computational methods available for evaluating binding affinities have adopted molecular mechanics‐based force fields, although they cannot fully describe protein–ligand interactions. A noteworthy computational method in development involves large‐scale electronic structure calculations. Fragment molecular orbital (FMO) method, which is one of such large‐scale calculation techniques, is applied in this study for calculating the binding energies between proteins and ligands. By testing the effects of specific FMO calculation conditions (including fragmentation size, basis sets, electron correlation, exchange‐correlation functionals, and solvation effects) on the binding energies of the FK506‐binding protein and 10 ligand complex molecule, we have found that the standard FMO calculation condition, FMO2‐MP2/6‐31G(d), is suitable for evaluating the protein–ligand interactions. The correlation coefficient between the binding energies calculated with this FMO calculation condition and experimental values is determined to be R = 0.77. Based on these results, we also propose a practical scheme for predicting binding affinities by combining the FMO method with the quantitative structure–activity relationship (QSAR) model. The results of this combined method can be directly compared with experimental binding affinities. The FMO and QSAR combined scheme shows a higher correlation with experimental data (R = 0.91). Furthermore, we propose an acceleration scheme for the binding energy calculations using a multilayer FMO method focusing on the protein–ligand interaction distance. Our acceleration scheme, which uses FMO2‐HF/STO‐3G:MP2/6‐31G(d) at R_int = 7.0 Å, reduces computational costs, while maintaining accuracy in the evaluation of binding energy. © 2015 Wiley Periodicals, Inc.     

Analytic first and second derivatives of the energy in the fragment molecular orbital method combined with molecular mechanics

Analytic first and second derivatives of the energy are developed for the fragment molecular orbital method interfaced with molecular mechanics in the electrostatic embedding scheme at the level of Hartree-Fock and density functional theory. The importance of the orbital response terms is demonstrated. The role of electrostatic embedding upon molecular vibrations is analyzed, comparing force field and quantum mechanical treatments for an ionic liquid and a solvated protein. The method is applied for 100 protein conformations sampled in molecular dynamics (MD) to take into account the complexity of a flexible protein structure in solution, and a good agreement with experimental data is obtained: Frequencies from an experimental infrared (IR) spectrum are reproduced within 17 cm⁻¹ .     

Moment expansion of the linear density‐density response function

We present a low rank moment expansion of the linear density‐density response function. The general interacting (fully nonlocal) density‐density response function is calculated by means of its spectral decomposition via an iterative Lanczos diagonalization technique within linear density functional perturbation theory. We derive a unitary transformation in the space of the eigenfunctions yielding subspaces with well‐defined moments. This transformation generates the irreducible representations of the density‐density response function with respect to rotations within SO(3). This allows to separate the contributions to the electronic response density from different multipole moments of the perturbation. Our representation maximally condenses the physically relevant information of the density‐density response function required for intermolecular interactions, yielding a considerable reduction in dimensionality. We illustrate the performance and accuracy of our scheme by computing the electronic response density of a water molecule to a complex interaction potential. © 2015 Wiley Periodicals, Inc.     

Specific interactions and binding free energies between thermolysin and dipeptides: molecular simulations combined with ab initio molecular orbital and classical vibrational analysis

Thermolysin (TLN) is a metalloprotease widely used as a nonspecific protease for sequencing peptide and synthesizing many useful chemical compounds by the chemical industry. It was experimentally shown that the activity and functions of TLN are inhibited by the binding of many types of amino acid dipeptides. However, the binding mechanisms between TLN and dipeptides have not been clarified at the atomic and electronic levels. In this study, we investigated the binding mechanisms between TLN and four dipeptides. Specific interactions and binding free energies (BFEs) between TLN and the dipeptides were calculated using molecular simulations based on classical molecular dynamics and ab initio fragment molecular orbital (FMO) methods. The molecular systems were embedded in solvating water molecules during calculations. The calculated BFEs were qualitatively consistent with the trend of the experimentally observed inhibition of TLN activity by binding of the dipeptides. In addition, the specific interactions between the dipeptides and each amino acid residue of TLN or solvating water molecules were elucidated by the FMO calculations.     

The self‐consistent charge density functional tight‐binding theory study of carbon adatoms using tuned Hubbard U parameters

The self‐consistent charge density functional tight‐binding (DFTB) theory is a useful tool for realizing the electronic structures of large molecular complex systems. In this study, the electronic structure of C₆₁ formed by fullerene C₆₀ with a carbon adatom is analyzed, using the fully localized limit and pseudo self‐interaction correction methods of DFTB to adjust the Hubbard U parameter (DFTB + U). The results show that both the methods used to adjust U can significantly reduce the molecular orbital energy of occupied states localized on the defect carbon atom and improve the gap between highest occupied molecular orbital(HOMO) and lowest unoccupied molecular orbital(LUMO) of C₆₁. This work will provide a methodological reference point for future DFTB calculations of the electronic structures of carbon materials.     

Shift‐and‐invert parallel spectral transformation eigensolver: Massively parallel performance for density‐functional based tight‐binding

The Shift‐and‐invert parallel spectral transformations (SIPs), a computational approach to solve sparse eigenvalue problems, is developed for massively parallel architectures with exceptional parallel scalability and robustness. The capabilities of SIPs are demonstrated by diagonalization of density‐functional based tight‐binding (DFTB) Hamiltonian and overlap matrices for single‐wall metallic carbon nanotubes, diamond nanowires, and bulk diamond crystals. The largest (smallest) example studied is a 128,000 (2000) atom nanotube for which ～330,000 (～5600) eigenvalues and eigenfunctions are obtained in ～190 (～5) seconds when parallelized over 266,144 (16,384) Blue Gene/Q cores. Weak scaling and strong scaling of SIPs are analyzed and the performance of SIPs is compared with other novel methods. Different matrix ordering methods are investigated to reduce the cost of the factorization step, which dominates the time‐to‐solution at the strong scaling limit. A parallel implementation of assembling the density matrix from the distributed eigenvectors is demonstrated. © 2015 Wiley Periodicals, Inc.     

Protein‐specific force field derived from the fragment molecular orbital method can improve protein–ligand binding interactions

Accurate computational estimate of the protein–ligand binding affinity is of central importance in rational drug design. To improve accuracy of the molecular mechanics (MM) force field (FF) for protein–ligand simulations, we use a protein‐specific FF derived by the fragment molecular orbital (FMO) method and by the restrained electrostatic potential (RESP) method. Applying this FMO‐RESP method to two proteins, dodecin, and lysozyme, we found that protein‐specific partial charges tend to differ more significantly from the standard AMBER charges for isolated charged atoms. We did not see the dependence of partial charges on the secondary structure. Computing the binding affinities of dodecin with five ligands by MM PBSA protocol with the FMO‐RESP charge set as well as with the standard AMBER charges, we found that the former gives better correlation with experimental affinities than the latter. While, for lysozyme with five ligands, both charge sets gave similar and relatively accurate estimates of binding affinities. © 2013 Wiley Periodicals, Inc.     

Possible mechanism of BN fullerene formation from a boron cluster: Density‐functional tight‐binding molecular dynamics simulations

We simulate the formation of a BN fullerene from an amorphous B cluster at 2000 K by quantum mechanical molecular dynamics based on the density‐functional tight‐binding method. We run 30 trajectories 200 ps in length, where N atoms are supplied around the target cluster, which is initially an amorphous B₃₆ cluster. Most of the incident N atoms are promptly incorporated into the target cluster to form B‐N‐B bridges or NB₃ pyramidal local substructures. BN fullerene formation is initiated by alternating BN ring condensation. Spontaneous atomic rearrangement and N₂ dissociation lead to the construction of an sp² single‐shelled structure, during which the BN cluster undergoes a transition from a liquid‐like to a solid‐like state. Continual atomic rearrangement and sporadic N₂ dissociation decrease the number of defective rings in the BN cluster and increase the number of six‐membered rings, forming a more regular shell structure. The number of four‐membered rings tends to remain constant, and contributes to more ordered isolated‐tetragon‐rule ring placement. © 2016 Wiley Periodicals, Inc.     

Ab initio quantum mechanical study of the binding energies of human estrogen receptor alpha with its ligands: an application of fragment molecular orbital method

We have theoretically examined the relative binding affinities (RBA) of typical ligands, 17beta-estradiol (EST), 17alpha-estradiol (ESTA), genistein (GEN), raloxifene (RAL), 4-hydroxytamoxifen (OHT), tamoxifen (TAM), clomifene (CLO), 4-hydroxyclomifene (OHC), diethylstilbestrol (DES), bisphenol A (BISA), and bisphenol F (BISF), to the alpha-subtype of the human estrogen receptor ligand-binding domain (hERalpha LBD), by calculating their binding energies. The ab initio fragment molecular orbital (FMO) method, which we have recently proposed for the calculations of macromolecules such as proteins, was applied at the HF/STO-3G level. The receptor protein was primarily modeled by 50 amino acid residues surrounding the ligand. The number of atoms in these model complexes is about 850, including hydrogen atoms. For the complexes with EST, RAL, OHT, and DES, the binding energies were calculated again with the entire ERalphaLBD consisting of 241 residues or about 4000 atoms. No significant difference was found in the calculated binding energies between the model and the real protein complexes. This indicates that the binding between the protein and its ligands is well characterized by the model protein with the 50 residues. The calculated binding energies relative to EST were very well correlated with the experimental RBA (the correlation coefficient r=0.837) for the ligands studied in this work. We also found that the charge transfer between ER and ligands is significant on ER-ligand binding. To our knowledge, this is the first achievement of ab initio quantum mechanical calculations of large molecules such as the entire ERalphaLBD protein.     

Theoretical study of the prion protein based on the fragment molecular orbital method

We performed fragment molecular orbital (FMO) calculations to examine the molecular interactions between the prion protein (PrP) and GN8, which is a potential curative agent for prion diseases. This study has the following novel aspects: we introduced the counterpoise method into the FMO scheme to eliminate the basis set superposition error and examined the influence of geometrical fluctuation on the interaction energies, thereby enabling rigorous analysis of the molecular interaction between PrP and GN8. This analysis could provide information on key amino acid residues of PrP as well as key units of GN8 involved in the molecular interaction between the two molecules. The present FMO calculations were performed using an original program developed in our laboratory, called “Parallelized ab initio calculation system based on FMO (PAICS)”. © 2009 Wiley Periodicals, Inc. J Comput Chem 2009     

Theoretical analysis of binding specificity of influenza viral hemagglutinin to avian and human receptors based on the fragment molecular orbital method

The hemagglutinin (HA) protein of the influenza virus binds to the host cell receptor in the early stage of viral infection. A change in binding specificity from avian 2-3 to human 2-6 receptor is essential for optimal human-to-human transmission and pandemics. Therefore, it is important to reveal the key factors governing the binding affinity of HA-receptor complex at the molecular level for the understanding and prediction of influenza pandemics. In this work, on the basis of ab initio fragment molecular orbital (FMO) method, we have carried out the interaction energy analysis of HA-receptor complexes to quantitatively elucidate the binding specificity of HAs to avian and human receptors. To discuss the binding property of influenza HA comprehensively, a number of HAs from human H1, swine H1, avian H3 and avian H5 viruses were analyzed. We performed detailed investigations about the interaction patterns of complexes of various HAs and receptor analogues, and revealed that intra-molecular interactions between conserved residues in HA play an important role for HA-receptor binding. These results may provide a hint to understand the role of conserved acidic residues at the receptor binding site which are destabilized by the electrostatic repulsion with sialic acid. The calculated binding energies and interaction patterns between receptor and HAs are consistent with the binding specificities of each HA and thus explain the receptor binding mechanism. The calculated results in the present analysis have provided a number of viewpoints regarding the models for the HA-receptor binding specificity associated with mutated residues. Examples include the role of Glu190 and Gln226 for the binding specificity of H5 HA. Since H5 HA has not yet been adapted to human receptor and the mechanism of the specificity change is unknown, this result is helpful for the prediction of the change in receptor specificity associated with forthcoming possible pandemics.     

1,2,3-三氮杂苯－(水)3复合物多体相互作用

The interaction between 1,2,3-triazine and three water molecules was studied using density functional theory B3LYP method at 6-31-t＋＋G^＊＊ basis set. Various structures for 1,2,3-triazine-（water）n （n= 1, 2, 3） complex were investigated and the different lower energy structures were reported. Many-body analysis was also carded out to obtain relaxation energy and many-body interaction energy （two, three, and four-body）, and the most stable conformer has the basis set superposition error corrected interaction energy of -- 102.61 kJ/mol. The relaxation energy, two- and three-body interactions have significant contribution to the total interaction energy whereas four-body interaction was very small for 1,2,3-triazine-（water）3 complex.     

Rapid and accurate assessment of GPCR–ligand interactions Using the fragment molecular orbital‐based density‐functional tight‐binding method

The reliable and precise evaluation of receptor–ligand interactions and pair‐interaction energy is an essential element of rational drug design. While quantum mechanical (QM) methods have been a promising means by which to achieve this, traditional QM is not applicable for large biological systems due to its high computational cost. Here, the fragment molecular orbital (FMO) method has been used to accelerate QM calculations, and by combining FMO with the density‐functional tight‐binding (DFTB) method we are able to decrease computational cost 1000 times, achieving results in seconds, instead of hours. We have applied FMO‐DFTB to three different GPCR–ligand systems. Our results correlate well with site directed mutagenesis data and findings presented in the published literature, demonstrating that FMO‐DFTB is a rapid and accurate means of GPCR–ligand interactions. © 2017 Authors. Journal of Computational Chemistry Published by Wiley Periodicals, Inc.