To obtain accurate electronic structures in large systems, we need to scale up high-level electronic structure calculations. Kohn-Sham density functional theory (DFT) is exact with the exact exchange-correlation (XC) potential. We developed a simple, yet effective method to direct construct accurate XC potentials in a system. The method consists of two steps: (a) the system partitioning and (b) the XC potential patching. We developed two schemes to partition the system: (1) partition the system’s density and (2) partition the system’s density matrix. Once the partitioning is finished, an embedded cluster is defined and its XC potential can be computed by inverting its electron density with advanced, orbital-based XC functionals, or correlated-wave function methods. The cluster’s XC potential is projected to its central atom, and the system’s XC potential is obtained by patching it atom by atom. We demonstrate the performance of this exchange-correlation potential patching (XCPP) method of patching the potential of a fully nonlocal XC functional: the exact-exchange (EX) and the correlation based on the random phase approximation (RPA). XCPP-RPA is applied to two one-dimensional systems: a H20 chain and a H10Li8 chain. In both tests, XCPP-RPA results converge to the benchmark as we increase the cluster sizes. We observed an effective error cancellation between the patched EX and RPA energies when the density partitioning is employed. The patched EX+RPA potentials agree well with the benchmarks. This work serves the first step toward self-consistent RPA simulations of large systems within the framework of XCPP.

Physisorption is an effective route to meet hydrogen gas (H2) storage and delivery requirements for transportation because it is fast and fully reversible under mild conditions. However, most current candidates have too small binding enthalpies to H2 which leads to volumetric capacity less than 10 g/L compared to that of the system target of 40 g/L at 298 K. Accurate quantum mechanical (QM) methods were used to determine the H2 binding enthalpy of 5 linkers which
were chelated with 11 different transition metals (Tm), including abundant first-row Tm (Sc through Cu), totaling 60 molecular compounds with more than 4 configurations related to the different number of H2 that interact with the molecular compound. It was found that first-row Tm gave similar and sometimes superior van der Waals interactions with H2 than precious Tm. Based on these linkers, 30 new covalent organic frameworks (COFs) were constructed. The H2 uptakes of these new COFs were determined using quantum mechanics (QM)-based force fields and grand canonical Monte Carlo (GCMC) simulations. For the first time, the range for the adsorption pressure was explored for 0-700 bar and 298 K. It was determined that Co-, Ni-, and Fe-based COFs can give high H2 uptake and delivery when compared to bulk H2 on this unexplored range of pressure.

The CASSCF approach enables an accurate description of the electronic structure of many-electron systems where non-dynamical correlation effects are important. Variational two-electron reduced-density-matrix (v2RDM) methods provide a route to polynomial-scaling implementations of CASSCF enabling the description of active-spaces larger than those that can be considered using configuration interaction-(CI-) based CASSCF. In this work, we present an implementation of analytic energy gradients for (v2RDM)-driven CASSCF. Expressions for analytic gradients are simplified by the fact that the Lagrangian for the active-space energy is stationary with respect to variations in the active-space reduced-density matrices. We assess the relative performance of v2RDM-CASSCF relative to CI-based CASSCF in the geometry optimization of 20 molecules. For these molecules, bond lengths from geometry optimization with v2RDM-driven CASSCF are in good agreement with those from (CI)-driven CASSCF. When enforcing two-particle N-representability conditions, v2RDM-CASSCF-optimized bond lengths display a mean unsigned error of 0.006 Å and a maximum unsigned error of 0.0265 Å, relative to those obtained with CI-CASSCF. When enforcing partial three-particle N-representability conditions, the mean and maximum errors are reduced to 0.0006 Å and 0.0054 Å, respectively.

The protonation/oxidation state of proteins and other biomolecules can be related to their structure and function, and it can affect properties like stability, ligand binding, catalysis, absorption spectrum, among others. This happens because pH and/or the redox potential affect the charge distribution on the biomolecules due to changes in the predominant protonation/oxidation state of the relevant groups. Also, ligands can change their protonation/oxidation state upon protein binding process or during an enzymatic reaction. Therefore, theoretical methods that can correctly describe the protonation/reduction state at constant pH and/or at constant redox potential are very important.
On this presentation, we show the implementation on AMBER of an extension to constant redox potential of codes already implemented: constant pH Molecular Dynamics (CpHMD) and Replica Exchange Molecular Dynamics along the pH-dimension (pH-REMD). Due to the similarity between the Henderson-Hasselbalch equation (applied to acid-base reactions) and the Nernst equation (applied to electrochemistry, reduction reactions), the mathematical derivations used for CpHMD and pH-REMD can be extended to the redox potential. By making use of CUDA implementation, we obtain a high-performance code that can be used on simulations of large systems. REMD is an important technique that enhances the statistical ensemble of a simulation while takes advantage of parallelism. The REMD implementation along both the pH and redox potential dimensions (pH,E-REMD) is important because several experimental measures are done both at constant pH and at constant redox potential.
We also show how our results are in agreement with theoretical/experimental expectations, and how computational benchmarks show the high-performance of calculations using GPU in comparison with serial or MPI calculations for large systems.

Human GAR-Tfase is a regulatory enzyme in de-novo purine biosynthesis which has been proven to be an anti-cancer target. Drugs such as Lometrexol, AG-2034 and pemetrexed have been designed using GAR-Tfase as a target enzyme. Folylpolyglutamate synthetase (FPGS) adds charged glutamates to these drugs which improves their binding affinities. However, the added charge leads to excessive cell retentivity and toxicity. E.coli GAR-Tfase, a well studied system, is sequentially similar to Human GAR-Tfase and most of its functional residues remain conserved. In the present work, we use E.coli GAR-Tfase to study pH-dependent conformational changes, ligand-binding and catalysis. We use pH replica exchange molecular dynamics (pH-REMD) simulations implemented in the AMBER suite to run in GPUs. We find that conformational changes in GAR-Tfase take place at low pH, which are consistent with experimental and previous theoretical studies. We propose refinements to the existing GAR-Tfase catalytic mechanism which has the potential to design effective inhibitors. To explain the effect of pH on catalysis, we present a pH-activity curve combining the population of catalytically competent protonation states and pH-dependent ligand-binding; which is consistent with experiments. This work can be applied to identify and extend the pH-range of drug molecules’ inhibition.

Intracellular membraneless organelles, formed by liquid-liquid phase separation (LLPS) of mixtures of proteins and possibly RNA, mediate myriad cellular functions. Cells use a variety of biochemical signals such as expression level and posttranslational modification to regulate the formation and dissolution of membraneless organelles during functional processes. We have developed a powerful computational method called FMAP [1,2] for determining the thermodynamic conditions for LLPS, where proteins separate into a low-concentration (“dissolved”) phase and a high-concentration (“droplet”) phase. FMAP involves calculating the excess chemical potentials of protein molecules over a wide range of concentrations. By using fast Fourier transform to efficiently evaluate protein-protein interactions, FMAP enables an atomistic representation of the protein molecules. Here we applied FMAP to three homologues of g-crystalin to elucidate why minor changes in amino-acid sequence can lead to drastic differences in critical temperature (T_c, the highest temperature at which LLPS exists). Our calculations reproduce the experimental observation that gB and gD have similar T_c but gF has a much higher T_c. The calculations further reveal that weak intermolecular binding through the ridge between the two domains of g-crystalin drives droplet formation. A single residue at position 130, a Ser in both gB and gD but a Trp in gF, that lines the inter-domain ridge makes a prominent contribution to the disparity in T_c. Our study opens the door to quantitative modeling of the regulation of membraneless organelle formation by biochemical signals.

1. S. Qin and H.-X. Zhou (2014). Further development of the FFT-based method for atomistic modeling of protein folding and binding under crowding: optimization of accuracy and speed. J. Chem. Theory Comput. 10, 2824-2835.
2. S. Qin and H.-X. Zhou (2016). Fast method for computing chemical potentials and liquid-liquid phase equilibria of macromolecular solutions. J. Phys. Chem. B. 120, 8164-8174.