Graphical processing units (GPUs) are specialized hardware for graphical renderings. Because video games and graphics renderings rely on simulations of mechanics and optics based on physical principles, it is not surprising that GPUs can also be used for general scientific simulations. Due to its fundamental differences from traditional processors, GPUs are best suited for executing a large group of similar tasks with minimal communication, and require completely different code design strategies. All of my methods have been optimized for GPUs, and it is now possible to use a single desktop workstation for calculations that used to require supercomputers. Such experience also prepares me for adapting quantum chemistry methods to other types of specialized novel hardware that may come out in the future.

This paper presents the polarizable quantum mechanics/molecular mechanics (QM/MM) embedding of state-averaged complete active space self-consistent field (SA-CASSCF) in the atomic multipole optimized energetics for biomolecular applications (AMOEBA) force field for the purpose of studying photoreactions in protein environments. We describe two extensions of our previous work that combines SA-CASSCF with AMOEBA water models allowing it to be generalized to AMOEBA models for proteins and other macromolecules. First, we discuss how our QM/MM model accounts for the discrepancy between the direct and polarization electric field that arises in the AMOEBA description of intramolecular polarization. A second improvement is the incorporation of double link-atom schemes to treat instances where the QM/MM boundary goes through covalent bonds. The new SA-CASSCF/AMOEBA method is applied to an elementary reaction step in NanoLuc, an artificial bioluminescence luciferase. We will show how the reaction mechanism is different when calculated in the gas phase, in polarizable continuum medium (PCM), versus in protein AMOEBA models. 

This paper introduces a spin-free formulation of the supporting subspace factorization, enabling a reduction in the computational scaling of extended multi-state complete active space second-order perturbation method (XMS-CASPT2) for arbitrary spins. We present a new derivation of the idea with the purpose of understanding its physical interpretation. The sources that make CASPT2 difficult are seperated into the “same-site interactions” and “inter-site interactions”. We first use the Kronecker sum to remove the same-site interactions in the absence of inter-site interactions, leading to MP2 energy in dressed orbitals. We then recover the inter-site interactions exactly using Lowdin partition, where the supporting subspace concept will naturally arise. 

This paper presents the state-averaged complete active space self-consistent field (SA-CASSCF) in AMOEBA polarizable water model, which enables rigorous simulation of non-adiabatic molecular dynamics with nonequilibrium solvation effects. When comparing the newly implemented AMOEBA with other solvent models (including PCM and fixed charge forcefield) for simulating photodynamics, AMOEBA is the only model that can capture a wide range of solvent effects, including the nonequilibrium electrostatic effects, hydrogen bonding, friction and steric hindrance. The paper also presents non-adiabatic molecular dynamics simulations of the model retinal protonated Schiff base molecule in water, where ultrafast transition to the ground state is observed.

In this paper, a new theoretical method is developed that enables simulating polarization effects of solvents over excited state dynamics for the first time. When applying the new method to dipolar merocyanine dye, the theoretical simulations revealed how the barrier heights of photoreaction pathways can be tuned with polarity of solvents, leading to systematic control of excited state lifetime. In another application to green fluorescent protein chromophore, the theoretical simulations show that the geometry and the shape of the conical intersection in solvents are qualitatively changed compared to in the gas phase, suggesting non-trivial solvation effects over excited state dynamics in photoreactions. 

This work combines automatic differentiation from computer algebra and diagrammatic notation from quantum computing, and showed that analytical gradients for complicated electronic structure methods can be automatically developed. The cost of the gradient calculation is shown to be roughly equal to 3 energy calculations. The tools developed in this work can significantly simplify future developments of analytical gradients and non-adiabatic couplings for new excited state methods. 

This work shows that by using supporting subspace method, the CASPT2 analytical gradients can be formulated as MP2 gradients and Fock gradients, thus significantly reduces the computational cost. The key is that the computational bottleneck in gradient calculations, i.e. the "lambda equation", takes the same form as the Bloch equation for solving energy.  This allows the supporting subspace factorization originally proposed for energy to be naturally applied in the gradient calculations. We present results from both ab initio molecular dynamics simulations as well as reaction pathway studies, which demonstrate the potential applications enabled by the new method.  

This work applied the supporting subspace formulation to multi-state CASPT2, and enabled accurate and efficient calculations of excited state energies. By combining this with QM/MM, the new method was applied to studies of green fluorescent proteins including a few hundred atoms in QM region on a single desktop. This work lays the foundation for future theoretical  studies of photochemical properties for large molecules. 

This work was the first introduction to the supporting subspace method, which enables the reformulation of multi-reference perturbation theory in an exact way, but is much simpler and faster. This work is also the first application of tensor hyper-contraction in the multi-reference context. In the new formulation,  the computational cost of CASPT2 has been reduced from O(N⁶) to O(N⁴).

This work was the first report of analytical gradients for tensor hyper-contraction. Using MP2 as an example, I showed that THC can also reduce the formal scaling for analytical gradient calculations of correlation methods to O(N⁴), in addition to reducing the scaling of energy calculations as already shown in previous works. 

I also proposed a new variant of THC in this work, called "dual grid THC", which is particularly appealing for gradient calculations due to its simple expressions and fast performance. 

I proposed a new THC variant called "local tensor hyper-contraction" in this work, which avoids the cubic scaling metric matrix inversion bottleneck in previously reported least-squares THC. In particular, in this work I discussed the close connection between THC and signal processing. One of the major achievements in this work is that we can now treat small proteins (e.g. ubiquitin) at SOS-MP2 level on a single desktop machine within a couple of hours. 

This work was the first report of GPU implementation of tensor hyper-contraction. In addition, this work discussed in details how spatial sparsity between atomic orbitals and gridpoints can be used to reduce the computational cost. This work was the first step towards making THC practical for large molecule calculations. 

This work introduced an Automated Code Engine (ACE), which can automatically generate optimized GPU kernels for molecular integrals starting from analytical expressions. The main new concept is the decision graph, which is motivated by graphs used in compiler theory. The decision graph is a directed colored graph: the nodes represent variables, while the connectivity and colors encode different operation-storage tradeoffs for computing the same integral expression. The optimal code variant is selected through empirical timing, which can naturally account for the differences between hardware generations.  

This work reported an efficient GPU implementation for ECP integrals required for descriptions of metal centers. I introduced a new upper-bound formula to estimate the sizes of the integrals, where integrals with upper-bound below a certain threshold will be screened and not computed. I also introduced a sorting strategy to balance the loads between GPU threads. The new implementation allows evaluating ECP integrals for 10,000 basis functions within a few minutes on a single GPU.