Linear Scaling Incremental Scheme for Correlation Energies with Embedding Generated Virtuals

A novel incremental scheme is presented including an incremental expansion of the virtual space for the calculation of electron correlation energies, which is compatible with any size-extensive correlation method and scales asymptotically linear for large molecules. The performance is studied for organic molecules, water clusters, and a La(III)-water complex, where the compatibility with pseudopotentials is also examined. The computational requirements are already reduced tremendously for medium-sized water clusters and hydrocarbons with respect to the canonical CCSD as well as the ordinary incremental scheme references. Correlation energies within chemical accuracy have been observed for all studied systems. The novelty of the method is that relatively small virtual spaces are used in combination with tuples of localized occupied spaces. The corresponding orthonormal occupied and virtual orbitals are obtained from QM/QM embedding calculations and can thus be used with standard quantum chemistry codes for correlation calculations. It is presented how relevant virtual spaces are selected and the correlation energies are linked in the new virtual space expansion.

Toward Accurate Prediction of Ion Mobility in Organic Semiconductors by Atomistic Simulation

A multiscale scheme (MLMS: Multi-Level Multi-Scale) to predict the ion mobility (μ) of amorphous organic semiconductors is proposed, which was successfully applied to the hole mobility predictions of 14 organic systems. An inverse relationship between μ and reorganization energy is observed due to local polaronic distortions. Another moderate inverse correlation between μ and distribution of site energy change exists, representing the effects of geometric flexibility. The former and the latter represent the intramolecular and intermolecular geometric effects, respectively. In addition, a linear correlation between transfer coupling and μ is observed, showing the importance of orbital overlaps between monomers. Especially, the highest hole mobility of C6-2TTN is due to its large transfer coupling. On the other hand, another high hole mobility of CBP turned out to come from the high first neighbor density (ρ_FND) of its first self-solvation, emphasizing the proper description of amorphous structural configurations with a sufficiently large number of monomers. In general, systems with either unusually high transfer coupling or high first neighbor density can potentially have high μ regardless of geometric effects. Especially, the newly suggested design parameter, ρ_FND, is pointing to a new direction as opposed to the traditional π-conjugation strategy.

Nonadiabatic Molecular Dynamics with Extended Density Functional Tight-Binding: Application to Nanocrystals and Periodic Solids

In this work, we report a new methodology for nonadiabatic molecular dynamics calculations within the extended tight-binding (xTB) framework. We demonstrate the applicability of the developed approach to finite and periodic systems with thousands of atoms by modeling "hot" electron relaxation dynamics in silicon nanocrystals and electron-hole recombination in both a graphitic carbon nitride monolayer and a titanium-based metal-organic framework (MOF). This work reports the nonadiabatic dynamic simulations in the largest Si nanocrystals studied so far by the xTB framework, with diameters up to 3.5 nm. For silicon nanocrystals, we find a non-monotonic dependence of "hot" electron relaxation rates on the nanocrystal size, in agreement with available experimental reports. We rationalize this relationship by a combination of decreasing nonadiabatic couplings related to system size and the increase of available coherent transfer pathways in systems with higher densities of states. We emphasize the importance of proper treatment of coherences for obtaining such non-monotonic dependences. We characterize the electron-hole recombination dynamics in the graphitic carbon nitride monolayer and the Ti-containing MOF. We demonstrate the importance of spin-adaptation and proper sampling of surface hopping trajectories in modeling such processes. We also assess several trajectory surface hopping schemes and highlight their distinct qualitative behavior in modeling the excited-state dynamics in superexchange-like models depending on how they handle coherences between nearly parallel states.

Local electronic structure analysis by ab initio elongation method: A benchmark using DNA block polymers

The ab initio elongation (ELG) method based on a polymerization concept is a feasible way to perform linear-scaling electronic structure calculations for huge aperiodic molecules while maintaining computational accuracy. In the method, the electronic structures are sequentially elongated by repeating (1) the conversion of canonical molecular orbitals (CMOs) to region-localized MOs (RLMOs), that is, active RLMOs localized onto a region close to an attacking monomer or frozen RLMOs localized onto the remaining region, and the subsequent (2) partial self-consistent-field calculations for an interaction space composed of the active RLMOs and the attacking monomer. For each ELG process, one can obtain local CMOs for the interaction space and the corresponding local orbital energies. Local site information, such as the local highest-occupied/lowest-unoccupied MOs, can be acquired with linear-scaling efficiency by correctly including electronic effects from the frozen region. In this study, we performed a local electronic structure analysis using the ELG method for various DNA block polymers with different sequential patterns. This benchmark aimed to confirm the effectiveness of the method toward the efficient detection of a singular local electronic structure in unknown systems as a future practical application. We discussed the high-throughput efficiency of our method and proposed a strategy to detect singular electronic structures by combining with a machine learning technique.

Towards a quantitative description of excitonic couplings in photosynthetic pigment-protein complexes: quantum chemistry driven multiscale approaches

A structure-based quantitative calculation of excitonic couplings between photosynthetic pigments has to describe the dynamical polarization of the protein/solvent environment of the pigments, giving rise to reaction field and screening effects. Here, this challenging problem is approached by combining the fragment molecular orbital (FMO) method with the polarizable continuum model (PCM). The method is applied to compute excitonic couplings between chlorophyll a (Chl a) pigments of the water-soluble chlorophyll-binding protein (WSCP). By calibrating the vacuum dipole strength of the 0–0 transition of the Chl a chromophores according to experimental data, an excellent agreement between calculated and experimental linear absorption and circular dichroism spectra of WSCP is obtained. The effect of the mutual polarization of the pigment ground states is calculated to be very small. The simple Poisson-Transition-charge-from-Electrostatic-potential (Poisson-TrEsp) method is found to accurately describe the screening part of the excitonic coupling, obtained with FMO/PCM. Taking into account that the reaction field effects of the latter method can be described by a scalar constant leads to an improvement of Poisson-TrEsp that is expected to provide the basis for simple and realistic calculations of optical spectra and energy transfer in photosynthetic light-harvesting complexes. In addition, we present an expression for the estimation of Huang–Rhys factors of high-frequency pigment vibrations from experimental fluorescence line-narrowing spectra that takes into account the redistribution of oscillator strength by the interpigment excitonic coupling. Application to WSCP results in corrected Huang–Rhys factors that are less than one third of the original values obtained by the standard electronic two-state analysis that neglects the above redistribution. These factors are important for the estimation of the dipole strength of the 0–0 transition of the chromophores and for the development of calculation schemes for the spectral density of the exciton-vibrational coupling.

The importance of reaction field and screening effects on the excitonic couplings is demonstrated, and from quantum-chemical calculations a single scaling factor is derived that can be used to improve simple models based on the Poisson equation.

Computational and data driven molecular material design assisted by low scaling quantum mechanics calculations and machine learning

Electronic structure methods based on quantum mechanics (QM) are widely employed in the computational predictions of the molecular properties and optoelectronic properties of molecular materials. The computational costs of these QM methods, ranging from density functional theory (DFT) or time-dependent DFT (TDDFT) to wave-function theory (WFT), usually increase sharply with the system size, causing the curse of dimensionality and hindering the QM calculations for large sized systems such as long polymer oligomers and complex molecular aggregates. In such cases, in recent years low scaling QM methods and machine learning (ML) techniques have been adopted to reduce the computational costs and thus assist computational and data driven molecular material design. In this review, we illustrated low scaling ground-state and excited-state QM approaches and their applications to long oligomers, self-assembled supramolecular complexes, stimuli-responsive materials, mechanically interlocked molecules, and excited state processes in molecular aggregates. Variable electrostatic parameters were also introduced in the modified force fields with the polarization model. On the basis of QM computational or experimental datasets, several ML algorithms, including explainable models, deep learning, and on-line learning methods, have been employed to predict the molecular energies, forces, electronic structure properties, and optical or electrical properties of materials. It can be conceived that low scaling algorithms with periodic boundary conditions are expected to be further applicable to functional materials, perhaps in combination with machine learning to fast predict the lattice energy, crystal structures, and spectroscopic properties of periodic functional materials.

Quantum-Mechanical Molecular Dynamics Simulations on Secondary Proton Transfer in Bacteriorhodopsin Using Realistic Models

Bacteriorhodopsin (BR) transports a proton from intracellular to extracellular (EC) sites through five proton transfers. The second proton transfer is the release of an excess proton stored in BR into the EC medium, and an atomistic understanding of this whole process has remained unexplored due to its ubiquitous environment. Here, fully quantum mechanical (QM) molecular dynamics (MD) and metadynamics (MTD) simulations for this process were performed at the divide-and-conquer density-functional tight-binding level using realistic models (∼50000 and ∼20000 atoms) based on the time-resolved photointermediate structures from an X-ray free electron laser. Regarding the proton storage process, the QM-MD/MTD simulations confirmed the Glu-shared mechanism, in which an excess proton is stored between Glu194 and Glu204, and clarified that the activation occurs by localizing the proton at Glu204 in the photocycle. Furthermore, the QM-MD/MTD simulations elucidated a release pathway from Glu204 through Ser193 to the EC water molecules and clarified that the proton release starts at ∼250 μs. In the ubiquitous proton diffusion in the EC medium, the transient proton receptors predicted experimentally were assigned to carboxylates in Glu9 and Glu74. Large-scale QM-MD/MTD simulations beyond the conventional sizes, which provided the above findings and confirmations, were possible by adopting our Dcdftbmd program.

Ab initio multi-level layered elongation method and its application to local interaction analysis between DNA bulge and ligand molecules

A multi-level layered elongation method was developed for efficiently analyzing the electronic states of local structures in large bio/nano-systems at the full ab initio level of theory. The original elongation method developed during the last three decades in our group has focused on the system in one direction from one terminal to the other terminal to sequentially construct the electronic states of a polymer, called a theoretical synthesis of polymers. In this study, an important region termed the central (C) part is targeted in a large polymer and the remainder are terminal (T) parts. The electronic structures along with polymer elongation are calculated repeatedly from both end T parts to the C central part at the same time. The important C part is treated with large basis sets (high level) and the other regions are treated with small basis sets (low level) in the ab initio theoretical framework. The electronic structures besides the C part can be reused for other systems with different structures at the C part, which renders the method computationally efficient. This multi-level layered elongation method was applied to the investigation on DNA single bulge recognition of small molecules (ligands). The reliability and validity of our approach were examined in comparison with the results obtained by direct calculations using a conventional quantum chemical method for the entire system. Furthermore, stabilization energies by the formation of the complex of bulge DNA and a ligand were estimated with basis set superposition error corrections incorporated into the elongation method.

Trajectory Surface Hopping Approach to Condensed-Phase Nonradiative Relaxation Dynamics Using Divide-and-Conquer Spin-Flip Time-Dependent Density-Functional Tight Binding

Nonradiative relaxation of excited molecules is central to many crucial issues in photochemistry. Condensed phases are typical contexts in which such problems are considered, and the nonradiative relaxation dynamics are expected to be significantly affected by interactions with the environment, for example, a solvent. We developed a nonadiabatic molecular dynamics simulation technique that can treat the nonradiative relaxation and explicitly include the environment in the calculations without a heavy computational burden. Specifically, we combined trajectory surface hopping with Tully's fewest-switches algorithm, a tight-binding approximated version of spin-flip time-dependent density-functional theory, and divide-and-conquer (DC) spatial fragmentation scheme. Numerical results showed that this method can treat systems with thousands of atoms within reasonable computational resources, and the error arising from DC fragmentation is negligibly small. Using this method, we obtained molecular insights into the solvent dependence of the photoexcited-state dynamics of trans-azobenzene, which demonstrate the importance of the environment for condensed-phase nonradiative relaxation.

Computational Methods for Biochemical Simulations Implemented in GAMESS

Computational methods for modeling biochemical processes implemented in GAMESS package are reviewed; in particular, quantum mechanics combined with molecular mechanics (QM/MM), semi-empirical, and fragmentation approaches. A detailed summary of capabilities is provided for the QM/MM implementation in QuanPol program and the fragment molecular orbital (FMO) method. Molecular modeling and visualization packages useful for biochemical simulations with GAMESS are described. GAMESS capabilities with corresponding references are tabulated for reader's convenience.

Density-Functional Tight-Binding Study of Carbonaceous Species Diffusion on the (100)-γ-Al2O3 Surface

Carbonaceous or oxy-carbon species are intermediates formed during C _x H _y combustion on a Pt _n /Al₂O₃ catalyst, which contain carbon, hydrogen, and oxygen atoms. The accumulation of the carbonaceous species, arguably, leads to catalytic deactivation; therefore, their removal is of importance. As the diffusion process is occasionally the rate-determining step in the growth of carbonaceous species, the present study aims to reveal the diffusion mechanisms. The free energy barriers of acetate, formate, and methoxy diffusion on the (100)-γ-Al₂O₃ surface were evaluated through extensive metadynamics simulations at the density-functional tight-binding level. The present work deduces that each adopted carbonaceous species exhibits different diffusion mechanisms and supports experimental evidence that the acetate species exhibits the slowest diffusivity among the adopted carbonaceous species.

Quantum mechanical molecular dynamics simulations of polaron formation in methylammonium lead iodide perovskite

Polaron formation in a halide perovskite is analyzed via nanometre-scale quantum mechanical molecular dynamics simulations.

Solvent Screening in Zwitterions Analyzed with the Fragment Molecular Orbital Method

Based on induced solvent charges, a new model of solvent screening is developed in the framework of the fragment molecular orbital combined with the polarizable continuum model. The developed model is applied to analyze interactions in a prototypical zwitterionic system, sodium chloride in water, and it is shown that the large underestimation of the interaction in the original solvent screening based on local charges is successfully corrected. The model is also applied to a complex of the Trp-cage (PDB: 1L2Y ) miniprotein with an anionic ligand, and the physical factors determined protein-ligand binding in solution are unraveled.

Combining Accuracy and Efficiency: An Incremental Focal-Point Method Based on Pair Natural Orbitals

In this work, we present a new pair natural orbitals (PNO)-based incremental scheme to calculate CCSD(T) and CCSD(T0) reaction, interaction, and binding energies. We perform an extensive analysis, which shows small incremental errors similar to previous non-PNO calculations. Furthermore, slight PNO errors are obtained by using T_PNO = T_TNO with appropriate values of 10^-7 to 10^-8 for reactions and 10^-8 for interaction or binding energies. The combination with the efficient MP2 focal-point approach yields chemical accuracy relative to the complete basis-set (CBS) limit. In this method, small basis sets (cc-pVDZ, def2-TZVP) for the CCSD(T) part are sufficient in case of reactions or interactions, while some larger ones (e.g., (aug)-cc-pVTZ) are necessary for molecular clusters. For these larger basis sets, we show the very high efficiency of our scheme. We obtain not only tremendous decreases of the wall times (i.e., factors >10²) due to the parallelization of the increment calculations as well as of the total times due to the application of PNOs (i.e., compared to the normal incremental scheme) but also smaller total times with respect to the standard PNO method. That way, our new method features a perfect applicability by combining an excellent accuracy with a very high efficiency as well as the accessibility to larger systems due to the separation of the full computation into several small increments.

Molecular Dipole Moments within the Incremental Scheme Using the Domain-Specific Basis-Set Approach

We present the first implementation of the fully automated incremental scheme for CCSD unrelaxed dipole moments using the domain-specific basis-set approach. Truncation parameters are varied, and the accuracy of the method is statistically analyzed for a test set of 20 molecules. The local approximations introduce small errors at second order and negligible ones at third order. For a third-order incremental CCSD expansion with a CC2 error correction, a cc-pVDZ/SV domain-specific basis set (tmain = 3.5 Bohr), and the truncation parameter f = 30 Bohr, we obtain a mean error of 0.00 mau (-0.20 mau) and a standard deviation of 1.95 mau (2.17 mau) for the total dipole moments (Cartesian components of the dipole vectors). By analyzing incremental CCSD energies, we demonstrate that the MP2 and CC2 error correction schemes are an exclusive correction for the domain-specific basis-set error. Our implementation of the incremental scheme provides fully automated computations of highly accurate dipole moments at reduced computational cost and is fully parallelized in terms of the calculation of the increments. Therefore, one can utilize the incremental scheme, on the same hardware, to extend the basis set in comparison to standard CCSD and thus obtain a better total accuracy.

Subsystem Analysis for the Fragment Molecular Orbital Method and Its Application to Protein-Ligand Binding in Solution

A subsystem analysis is derived incorporating interfragment interactions into the fragment properties, such as energies or charges. The relative stabilities of three alanine isomers, the α-helix, the β-turn, and the extended form are studied and the differences in fragment properties are elucidated. The analysis is further elaborated for studies of binding energies. The binding of the Trp-cage protein (PDB: 1L2Y ) to two ligands is studied in detail. Binding energies defined for each fragment can be used as a convenient descriptor for analyzing contributions to binding in solution.

The fragment molecular orbital method combined with density-functional tight-binding and the polarizable continuum model

The electronic gap in proteins is analyzed in detail, and it is shown that FMO-DFTB/PCM is efficient and accurate in describing the molecular structure of proteins in solution.

First UHF Implementation of the Incremental Scheme for Open-Shell Systems

The incremental scheme makes it possible to compute CCSD(T) correlation energies to high accuracy for large systems. We present the first extension of this fully automated black-box approach to open-shell systems using an Unrestricted Hartree-Fock (UHF) wave function, extending the efficient domain-specific basis set approach to handle open-shell references. We test our approach on a set of organic and metal organic structures and molecular clusters and demonstrate standard deviations from canonical CCSD(T) values of only 1.35 kJ/mol using a triple ζ basis set. We find that the incremental scheme is significantly more cost-effective than the canonical implementation even for relatively small systems and that the ease of parallelization makes it possible to perform high-level calculations on large systems in a few hours on inexpensive computers. We show that the approximations that make our approach widely applicable are significantly smaller than both the basis set incompleteness error and the intrinsic error of the CCSD(T) method, and we further demonstrate that incremental energies can be reliably used in extrapolation schemes to obtain near complete basis set limit CCSD(T) reaction energies for large systems.

Electron-correlated fragment-molecular-orbital calculations for biomolecular and nano systems

Recent developments in the fragment molecular orbital (FMO) method for theoretical formulation, implementation, and application to nano and biomolecular systems are reviewed. The FMO method has enabled ab initio quantum-mechanical calculations for large molecular systems such as protein-ligand complexes at a reasonable computational cost in a parallelized way. There have been a wealth of application outcomes from the FMO method in the fields of biochemistry, medicinal chemistry and nanotechnology, in which the electron correlation effects play vital roles. With the aid of the advances in high-performance computing, the FMO method promises larger, faster, and more accurate simulations of biomolecular and related systems, including the descriptions of dynamical behaviors in solvent environments. The current status and future prospects of the FMO scheme are addressed in these contexts.

Incremental evaluation of coupled cluster dipole polarizabilities

In this work we present the first implementation of the incremental scheme for coupled cluster linear-response frequency-dependent dipole polarizabilities.

Divide-and-conquer quantum mechanical material simulations with exascale supercomputers

Recent developments in large-scale materials science simulations, especially under the divide-and-conquer method, are reviewed. The pros and cons of the divide-and-conquer method are discussed. It is argued that the divide-and-conquer method, such as the linear-scaling 3D fragment method, is an ideal approach to take advantage of the heterogeneous architectures of modern-day supercomputers despite their relatively large prefactors among linear-scaling methods. Some developments in graphics processing unit (GPU) electronic structure calculations are also reviewed. The accelerators like GPU could be an essential part for the future exascale supercomputing.

The use of many-body expansions and geometry optimizations in fragment-based methods

Conspectus Chemists routinely work with complex molecular systems: solutions, biochemical molecules, and amorphous and composite materials provide some typical examples. The questions one often asks are what are the driving forces for a chemical phenomenon? How reasonable are our views of chemical systems in terms of subunits, such as functional groups and individual molecules? How can one quantify the difference in physicochemical properties of functional units found in a different chemical environment? Are various effects on functional units in molecular systems additive? Can they be represented by pairwise potentials? Are there effects that cannot be represented in a simple picture of pairwise interactions? How can we obtain quantitative values for these effects? Many of these questions can be formulated in the language of many-body effects. They quantify the properties of subunits (fragments), referred to as one-body properties, pairwise interactions (two-body properties), couplings of two-body interactions described by three-body properties, and so on. By introducing the notion of fragments in the framework of quantum chemistry, one obtains two immense benefits: (a) chemists can finally relate to quantum chemistry, which now speaks their language, by discussing chemically interesting subunits and their interactions and (b) calculations become much faster due to a reduced computational scaling. For instance, the somewhat academic sounding question of the importance of three-body effects in water clusters is actually another way of asking how two hydrogen bonds affect each other, when they involve three water molecules. One aspect of this is the many-body charge transfer (CT), because the charge transfers in the two hydrogen bonds are coupled to each other (not independent). In this work, we provide a generalized view on the use of many-body expansions in fragment-based methods, focusing on the general aspects of the property expansion and a contraction of a many-body expansion in a formally two-body series, as exemplified in the development of the fragment molecular orbital (FMO) method. Fragment-based methods have been very successful in delivering the properties of fragments, as well as the fragment interactions, providing insights into complex chemical processes in large molecular systems. We briefly review geometry optimizations performed with fragment-based methods and present an efficient geometry optimization method based on the combination of FMO with molecular mechanics (MM), applied to the complex of a subunit of protein kinase 2 (CK2) with a ligand. FMO results are discussed in comparison with experimental and MM-optimized structures.

Extension of accompanying coordinate expansion and recurrence relation method for general-contraction basis sets

An algorithm of the accompanying coordinate expansion and recurrence relation (ACE-RR), which is used for the rapid evaluation of the electron repulsion integral (ERI), has been extended to the general-contraction (GC) scheme. The present algorithm, denoted by GC-ACE-RR, is designed for molecular calculations including heavy elements, whose orbitals consist of many primitive functions with and without higher angular momentum such as d- and f-orbitals. The performance of GC-ACE-RR was assessed for (ss|ss)-, (pp|pp)-, (dd|dd)-, and (ff|ff)-type ERIs in terms of contraction length and the number of GC orbitals. The present algorithm was found to reduce the central processing unit time compared with the ACE-RR algorithm, especially for higher angular momentum and highly contracted orbitals. Compared with HONDOPLUS and GAMESS program packages, GC-ACE-RR computations for ERIs of three-dimensional gold clusters Aun (n = 1, 2, …, 10, 15, 20, and 25) are more than 10 times faster.