Analytic gradient for second order Møller-Plesset perturbation theory with the polarizable continuum model based on the fragment molecular orbital method

A new energy expression is proposed for the fragment molecular orbital method interfaced with the polarizable continuum model (FMO/PCM). The solvation free energy is shown to be more accurate on a set of representative polypeptides with neutral and charged residues, in comparison to the original formulation at the same level of the many-body expansion of the electrostatic potential determining the apparent surface charges. The analytic first derivative of the energy with respect to nuclear coordinates is formulated at the second-order M?ller-Plesset (MP2) perturbation theory level combined with PCM, for which we derived coupled perturbed Hartree-Fock equations. The accuracy of the analytic gradient is demonstrated on test calculations in comparison to numeric gradient. Geometry optimization of the small Trp-cage protein (PDB: 1L2Y) is performed with FMO/PCM/6-31(+)G(d) at the MP2 and restricted Hartree-Fock with empirical dispersion (RHF/D). The root mean square deviations between the FMO optimized and NMR experimental structure are found to be 0.414 and 0.426 A? for RHF/D and MP2, respectively. The details of the hydrogen bond network in the Trp-cage protein are revealed.     

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     

Extending the power of quantum chemistry to large systems with the fragment molecular orbital method

Following the brief review of the modern fragment-based methods and other approaches to perform quantum-mechanical calculations of large systems, the theoretical development of the fragment molecular orbital method (FMO) is covered in detail, with the emphasis on the physical properties, which can be computed with FMO. The FMO-based polarizable continuum model (PCM) for treating the solvent effects in large systems and the pair interaction energy decomposition analysis (PIEDA) are described in some detail, and a range of applications of FMO to biological studies is introduced. The factors determining the relative stability of polypeptide conformers (alpha-helix, beta-turn, and extended form) are elucidated using FMO/PCM and PIEDA, and the interactions in the Trp-cage miniprotein construct (PDB: 1L2Y) are analyzed using PIEDA.     

Second order Møller-Plesset perturbation theory based upon the fragment molecular orbital method

The fragment molecular orbital (FMO) method was combined with the second order M?ller-Plesset (MP2) perturbation theory. The accuracy of the method using the 6-31G(*) basis set was tested on (H(2)O)(n), n=16,32,64; alpha-helices and beta-strands of alanine n-mers, n=10,20,40; as well as on (H(2)O)(n), n=16,32,64 using the 6-31 + + G(**) basis set. Relative to the regular MP2 results that could be afforded, the FMO2-MP2 error in the correlation energy did not exceed 0.003 a.u., the error in the correlation energy gradient did not exceed 0.000 05 a.u./bohr and the error in the correlation contribution to dipole moment did not exceed 0.03 debye. An approximation reducing computational load based on fragment separation was introduced and tested. The FMO2-MP2 method demonstrated nearly linear scaling and drastically reduced the memory requirements of the regular MP2, making possible calculations with several thousands basis functions using small Pentium clusters. As an example, (H(2)O)(64) with the 6-31 + + G(**) basis set (1920 basis functions) can be run in 1 Gbyte RAM and it took 136 s on a 40-node Pentium4 cluster.     

Polarizable continuum model with the fragment molecular orbital-based time-dependent density functional theory

The polarizable continuum model (PCM) for describing the solvent effect was combined with the fragment molecular orbital-based time-dependent density functional theory (TDDFT). Several levels of the many-body expansion were implemented, and the importance of the many-body contributions to the singlet-excited states was discussed. To calibrate the accuracy, we performed a number of the model calculations using our method and the regular TDDFT in solution, applying them to phenol and polypeptides at the long-range corrected BLYP/6-31G* level. It was found that for systems up to 192 atoms the largest error in the excitation energy was 0.006 eV (vs. the regular TDDFT/PCM of the full system). The solvent shifts and the conformer effects were discussed, and the scaling was found to be nearly linear. Finally, we applied our method to the lowest singlet excitation of the photoactive yellow protein (PYP) in aqueous solution and determined the excitation energy to be in reasonable agreement with experiment. The excitation energy analysis provided the contributions of individual residues, and the main factors as well as their solvent shifts were determined.     

Change in a protein's electronic structure induced by an explicit solvent: an ab initio fragment molecular orbital study of ubiquitin

The effect of solvation on the electronic structure of the ubiquitin protein was analyzed using the ab initio fragment molecular orbital (FMO) method. FMO calculations were performed for the protein in vacuo, and the protein was immersed in an explicit solvent shell as thick as 12 A at the HF or MP2 level by using the 6-31G* basis set. The protein's physical properties examined were the net charge, the dipole moment, the internal energy, and the solvent interaction energy. Comparison of the computational results revealed the following changes in the protein upon solvation. First, the positively charged amino acid residues on the protein surface drew electrons from the solvent, while the negatively charged ones transfer electrons to the solvent. Second, the dipole moment of the protein was enhanced as a result of the polarization. Third, the internal energy of the protein was destabilized, but the destabilization was more than compensated for by the generation of a favorable protein-solvent interaction. Finally, the energetic changes were elicited both by the electron correlation effect of the first solvent shell and by the electrostatic effect of more distant solvent molecules. These findings were consistent with the picture of the solvated protein being a polarizable molecule dissolved in a dielectric media.     

The extension of the fragment molecular orbital method with the many-particle Green's function

By using the many-particle Green's function (GF) the extension of the fragment molecular orbital (FMO) method by Kitaura et al. [Chem. Phys. Lett. 313, 701 (1999)] is proposed. It is shown that the partial summation of the cluster expansion of GF reproduces the same extrapolation formula as that of FMO. Therefore we can determine the excitation energy, the transition moment, and the linear response of a molecule from GF approximated with the FMO procedure. It is also shown that no wave function exists which is consistent to the FMO results. The perturbation expansion in which the self-consistent charge approximation defines the unperturbed state is reported. By using it the three-body effects missing in the pair approximation of FMO are analyzed and the corrections to the energy and the reduced density matrices are proposed. In contrast to the previous works these new corrections are not expressed as the addition or the subtraction of the energies of fragments. They are size extensive and require only the quantities available by the FMO calculation. The accuracy of these corrections is validated with the extended Hubbard model and the several test molecules.     

Energy decomposition analysis in solution based on the fragment molecular orbital method

We develop the pair interaction energy decomposition analysis (PIEDA) in solution by combining the fragment molecular orbital (FMO) method with the polarizable continuum model (PCM). The solvent screening of the electrostatic interaction and the desolvation penalty in complex formation are described by this approach from ab initio calculations of fragments and their pairs. The applications to the complex of solvated sodium and chlorine ions, as well as to lysine and aspartic acid, show how the analysis helps reveal the physical picture. The PIEDA/PCM method is also applied to a small protein chignolin (PDB: 1UAO), and the solvent screening of the pair interactions is discussed.     

A study on orientation and absorption spectrum of interfacial molecules by using continuum model

In this work, a numerical procedure based on the continuum model is developed and applied to the solvation energy for ground state and the spectral shift against the position and the orientation of the interfacial molecule. The interface is described as a sharp boundary separating two bulk media. The polarizable continuum model (PCM) allows us to account for both electrostatic and nonelectrostatic solute-solvent interactions when we calculate the solvation energy. In this work we extend PCM to the interfacial system and the information about the position and orientation of the interfacial molecule can be obtained. Based on the developed expression of the electrostatic free energy of a nonequilibrium state, the numerical procedure has been implemented and used to deal with a series of test molecules. The time-dependent density functional theory (TDDFT) associated with PCM is used for the electron structure and the spectroscopy calculations of the test molecules in homogeneous solvents. With the charge distribution of the ground and excited states, the position- and orientation-dependencies of the solvation energy and the spectrum have been investigated for the interfacial systems, taking the electrostatic interaction, the cavitation energy, and the dispersion-repulsion interaction into account. The cavitation energy is paid particular attention, since the interface portion cut off by the occupation of the interfacial molecule contributes an extra part to the stabilization for the interfacial system. The embedding depth, the favorable orientational angle, and the spectral shift for the interfacial molecule have been investigated in detail. From the solvation energy calculations, an explanation has been given on why the interfacial molecule, even if symmetrical in structure, tends to take a tilting manner, rather than perpendicular to the interface.     

Covalent bond fragmentation suitable to describe solids in the fragment molecular orbital method

To improve the accuracy of the fragment molecular orbital method (FMO), we introduce a new fragmentation scheme based on using frozen orbitals to describe fractioned bonds. By applying this scheme to a set of polyalanine systems of up to 40 residues for the alpha-helix and beta-strand isomers, we established its accuracy, which is considerably improved compared to the original hybrid orbital projection method used for detaching bonds in FMO. For instance, at the two-body FMO expansion with the 6-311G* basis set, the error was typically reduced 2-4 times, and for 6-31G* the accuracy increase was even larger (10 times in terms of the maximum error). For the Trp-cage protein (PDB file 1L2Y) with many charged residues, a fairly large error was observed, which was shown to become small with a larger fragment size or at the three-body level. Consequently, we applied the new scheme to the adsorption of toluene and phenol on a faujasite zeolite, and we demonstrated that good accuracy can be achieved in reproducing ab initio results.     

A theoretical study on the origin of cooperativity in the formation of 3(10)- and alpha-helices.

By using a simple repeating unit method, we have conducted a theoretical study which delineates the preferences for beta-strand, 2(7)-ribbon, 3(10)-helix, and alpha-helix formation for a series of polyglycine models up to 14 amino acid residues (Ac-(Gly)(n), n = 0, 1, 2,., 14). Interactions among residues, which result in cooperativity, are clearly indicated by variations in calculated energies of the residues. Whereas no cooperativity is found in the formation of beta-strands and 2(7)-ribbons, there is a significant cooperativity in the formation of 3(10)- and alpha-helices, especially for the latter. In the case of alpha-helices, the 14th residue is more stable than the 3rd by about 3 kcal/mol. A good correlation between calculated residue energy and residue dipole moment was uncovered, indicating the importance of long-range electrostatic interactions to the cooperativity. The results of our calculations are compared with those of the AMBER and PM3 methods, and indicate that both methods, AMBER and PM3, need further development in the cooperative view of electrostatic interactions. The result should be of importance in providing insight into protein folding and formation of helical structures in a variety of polymeric compounds. This also suggests a strategy for the development of more consistent molecular mechanics force fields.     

Ab initio fragment molecular orbital (FMO) method applied to analysis of the ligand-protein interaction in a pheromone-binding protein

Full quantum computation of the electronic state of proteins has recently become possible by the advent of the ab initio fragment molecular orbital (FMO) method. We applied this method to the analysis of the interaction between the Bombyx mori pheromone-binding protein and its ligand, bombykol. The protein–ligand interaction of this molecular complex was minutely analyzed by the FMO method, and the analysis revealed several important interactions between the ligand and amino acid residues.     

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.     

Proton Solvation in Nonaqueous HF as a Main Factor That Determines Its Superacidic Properties

The Hartree–Fock method and electron correlation at the MP2 level in the 6-31++G** basis set is used together with the combined method for taking into account solvation (the discrete model of microsolvation and the polarizable continuum model (PCM)) to study the following molecular systems in the medium of liquid HF: [HFH]⁺, [HFHFH]⁺, [HF...HFH...FH]⁺, [HF...HFHFH...FH]⁺, and [HF...HF...HFH...FH...FH]⁺. The results of calculation make it possible to explain the anomalously high activity of protons in nonaqueous HF by a lower solvation energy of protons in liquid HF (by 48 kcal/mol) than in water. A conclusion is drawn that the combined method used in this work (the discrete model of microsolvation and PCM) is an efficient tool in estimating the contribution from the solvation to the thermodynamic parameters of the reaction in liquid HF.     

Fe(II) hydrolysis in aqueous solution: A DFT study

Species arising from Fe(II) hydrolysis in aqueous solution have been investigated using density-functional methods (DFT). The different tautomers and multiplicities of each species have been calculated. The solvation energy has been estimated using the UAHF–PCM method. The hydrolysis free energies have been estimated and compared with the available experimental data. The different hydrolysis species have distinct geometries and electronic structures. The estimated ionization potential of the hydrolyzed species is linearly dependent to the number of hydroxyls present in the complex. The estimated Fe(II)/Fe(III) oxidation potential is in good agreement with previously published results about 0.29 V larger than the experimental value. The results highlight the importance of the chemical speciation in describing electron transfer processes at a molecular level. The PBE/TZVP/UAHF–PCM method has been found to describe correctly the hydrolysis free energies of Fe(II) with an average error about 5 kcal mol⁻¹ from the experimental values.     

Excitation energy transfer (EET) between molecules in condensed matter: a novel application of the polarizable continuum model (PCM)

We present a quantum-mechanical theory to study excitation energy transfers between molecular systems in solution. The model is developed within the time-dependent (TD) density-functional theory and the solvent effects are introduced in terms of the polarizable continuum model (PCM). Unique characteristic of this model is that both "reaction field" and screening effects are included in a coherent and self-consistent way. This is obtained by introducing proper solvent-specific operators in the Kohn-Sham equations and in the corresponding TD scheme. The solvation model exploits the integral equation formalism (IEF) version of PCM and it defines the solvent operators on a molecular cavity modeled on the real three-dimensional (3D) structure of the solute systems. Applications to EET in dimers of ethylene and naphtalene are presented and discussed.     

Fragment molecular orbital study of the electronic excitations in the photosynthetic reaction center of Blastochloris viridis

All electron calculations were performed on the photosynthetic reaction center of Blastochloris viridis, using the fragment molecular orbital (FMO) method. The protein complex of 20,581 atoms and 77,754 electrons was divided into 1398 fragments, and the two‐body expansion of FMO/6‐31G* was applied to calculate the ground state. The excited electronic states of the embedded electron transfer system were separately calculated by the configuration interaction singles approach with the multilayer FMO method. Despite the structural symmetry of the system, asymmetric excitation energies were observed, especially on the bacteriopheophytin molecules. The asymmetry was attributed to electrostatic interaction with the surrounding proteins, in which the cytoplasmic side plays a major role. © 2009 Wiley Periodicals, Inc. J Comput Chem 2010     

Time-dependent density functional theory based upon the fragment molecular orbital method

Time-dependent density functional theory (TDDFT) was combined with the two-body fragment molecular orbital method (FMO2). In this FMO2-TDDFT scheme, the system is divided into fragments, and the electron density for fragments is determined self-consistently. Consequently, only one main fragment of interest and several fragment pairs including it are calculated by TDDFT. To demonstrate the accuracy of FMO2-TDDFT, we computed several low-lying singlet and triplet excited states of solvated phenol and polyalanine using our method and the standard TDDFT for the full system. The BLYP functional with the long-range correction (LC-BLYP) was employed with the 6-31G(*) basis set (some tests were also performed with 6-311G(*), as well as with B3LYP and time-dependent Hartree-Fock). Typically, FMO2-TDDFT reproduced the full TDDFT excitation energies within 0.1 eV, and for one excited state the error was about 0.2 eV. Beside the accurate reproduction of the TDDFT excitation energies, we also automatically get an excitation energy decomposition analysis, which provides the contributions of individual fragments. Finally, the efficiency of our approach was exemplified on the LC-BLYP6-31G(*) calculation of the lowest singlet excitation of the photoactive yellow protein which consists of 1931 atoms, and the obtained value of 3.1 eV is in agreement with the experimental value of 2.8 eV.     

Analytic gradient and molecular dynamics simulations using the fragment molecular orbital method combined with effective potentials

The completely analytic energy gradients are derived and implemented for the two-body fragment molecular orbital (FMO2) method combined with the model core potentials (MCP) and effective fragment potentials (EFP). The many-body terms in EFP require solving coupled-perturbed Hartree-Fock equations, which are derived and implemented. The molecular dynamics (MD) simulations are performed using the FMO2/MCP method for the capped alanine decamer and with the FMO2/EFP method for the zwitterionic conformer of glycine tetramer immersed in the water layer of 6.0 Å (135 water molecules). The results of the MD simulations using the FMO2/EFP and FMO2/MCP gradients show that the total energy is conserved at the time steps less than 1 fs.     

Study of multipole contributions to the structure of water around ions in solution using the soft sticky dipole-quadrupole-octupole (SSDQO) model of water

The solvation of ions in the soft sticky dipole-quadrupole-octupole (SSDQO) model for liquid water is presented here. This new potential energy function for liquid water describes water-water interactions by a Lennard-Jones term plus a sticky potential consisting of an approximate moment expansion with point dipole, quadrupole, and octupole moments. The SSDQO potential energy function using the moments from extended simple point charge (SPC/E), TIP3P, or TIP5P reproduces the pair potential energy functions and radial distribution functions of the respective multipoint model but it is much faster than even the three-point models. Here, the solvation of ions in SSDQO water is studied using ion-water potential energy functions consisting of moment expansions up to the charge-quadrupole term, up to the charge-octupole term, and up to an approximate charge-hexadecapole term using the moments of SPC/E water. The radial distributions from Monte Carlo simulations show the best agreement with the results for ions in SPC/E water for the expansion up to the charge-hexadecapole term. Thus, the best results are obtained when the water-water and ion-water potentials are exact up to the 1r(4) term and also contain an approximate 1r(5) term. Overall, the simplicity, efficiency, and accuracy of the SSDQO potential make it potentially very useful for computer simulations of aqueous solvation.