The OpenMolcas Web: A Community-Driven Approach to Advancing Computational Chemistry.

The developments of the open-source OpenMolcas chemistry software environment since spring 2020 are described, with a focus on novel functionalities accessible in the stable branch of the package or via interfaces with other packages. These developments span a wide range of topics in computational chemistry and are presented in thematic sections: electronic structure theory, electronic spectroscopy simulations, analytic gradients and molecular structure optimizations, ab initio molecular dynamics, and other new features. This report offers an overview of the chemical phenomena and processes OpenMolcas can address, while showing that OpenMolcas is an attractive platform for state-of-the-art atomistic computer simulations.

Automated Active Space Selection with Dipole Moments.

Multireference calculations can provide accurate information of systems with strong correlation, which have increasing importance in the development of new molecules and materials. However, selecting a suitable active space for multireference calculations is nontrivial, and the selection of an unsuitable active space can sometimes lead to results that are not physically meaningful. Active space selection often requires significant human input, and the selection that leads to reasonable results often goes beyond chemical intuition. In this work, we have developed and evaluated two protocols for automated selection of the active space for multireference calculations based on a simple physical observable, the dipole moment, for molecules with nonzero ground-state dipole moments. One protocol is based on the ground-state dipole moment, and the other is based on the excited-state dipole moments. To evaluate the protocols, we constructed a dataset of 1275 active spaces from 25 molecules, each with 51 active space sizes considered, and have mapped out the relationship between the active space, dipole moments, and vertical excitation energies. We have demonstrated that, within this dataset, our protocols allow one to choose among a number of accessible active spaces one that is likely to give reasonable vertical excitation energies, especially for the first three excitations, with no parameters manually decided by the user. We show that, with large active spaces removed from consideration, the accuracy is similar and the time-to-solution can be reduced by more than 10 fold. We also show that the protocols can be applied to potential energy surface scans and determining the spin states of transition metal oxides.

A mechanistic study on the cellular uptake, intracellular trafficking, and antisense gene regulation of bottlebrush polymer-conjugated oligonucleotides.

We have developed a non-cationic transfection vector in the form of bottlebrush polymer-antisense oligonucleotide (ASO) conjugates. Termed pacDNA (polymer-assisted compaction of DNA), these agents show improved biopharmaceutical characteristics and antisense potency in vivo while suppressing non-antisense side effects. Nonetheless, there still is a lack of the mechanistic understanding of the cellular uptake, subcellular trafficking, and gene knockdown with pacDNA. Here, we show that the pacDNA enters human non-small cell lung cancer cells (NCI-H358) predominantly by scavenger receptor-mediated endocytosis and macropinocytosis and trafficks via the endolysosomal pathway within the cell. The pacDNA significantly reduces a target gene expression (KRAS) in the protein level but not in the mRNA level, despite that the transfection of certain free ASOs causes ribonuclease H1 (RNase H)-dependent degradation of KRAS mRNA. In addition, the antisense activity of pacDNA is independent of ASO chemical modification, suggesting that the pacDNA always functions as a steric blocker.

A Long-Circulating Vector for Aptamers Based upon Polyphosphodiester-Backboned Molecular Brushes.

Aptamers face challenges for use outside the ideal conditions in which they are developed. These difficulties are most palpable in vivo due to nuclease activities, rapid clearance, and off-target binding. Herein, we demonstrate that a polyphosphodiester-backboned molecular brush can suppress enzymatic digestion, reduce non-specific cell uptake, enable long blood circulation, and rescue the bioactivity of a conjugated aptamer in vivo. The backbone along with the aptamer is assembled via solid-phase synthesis, followed by installation of poly(ethylene glycol) (PEG) side chains using a two-step process with near-quantitative efficiency. The synthesis allows for precise control over polymer size and architecture. Consisting entirely of building blocks that are generally recognized as safe for therapeutics, this novel molecular brush is expected to provide a highly translatable route for aptamer-based therapeutics.

Targeting oncogenic KRAS with molecular brush-conjugated antisense oligonucleotides.

The mutant form of the guanosine triphosphatase (GTPase) KRAS is a key driver in human tumors but remains a challenging therapeutic target, making KRASMUT cancers a highly unmet clinical need. Here, we report a class of bottlebrush polyethylene glycol (PEG)-conjugated antisense oligonucleotides (ASOs) for potent in vivo KRAS depletion. Owing to their highly branched architecture, these molecular nanoconstructs suppress nearly all side effects associated with DNA-protein interactions and substantially enhance the pharmacological properties of the ASO, such as plasma pharmacokinetics and tumor uptake. Systemic delivery to mice bearing human non-small-cell lung carcinoma xenografts results in a significant reduction in both KRAS levels and tumor growth, and the antitumor performance well exceeds that of current popular ASO paradigms, such as chemically modified oligonucleotides and PEGylation using linear or slightly branched PEG. Importantly, these conjugates relax the requirement on the ASO chemistry, allowing unmodified, natural phosphodiester ASOs to achieve efficacy comparable to that of chemically modified ones. Both the bottlebrush polymer and its ASO conjugates appear to be safe and well tolerated in mice. Together, these data indicate that the molecular brush-ASO conjugate is a promising therapeutic platform for the treatment of KRAS-driven human cancers and warrant further preclinical and clinical development.

Machine learning dielectric screening for the simulation of excited state properties of molecules and materials.

Accurate and efficient calculations of absorption spectra of molecules and materials are essential for the understanding and rational design of broad classes of systems. Solving the Bethe-Salpeter equation (BSE) for electron-hole pairs usually yields accurate predictions of absorption spectra, but it is computationally expensive, especially if thermal averages of spectra computed for multiple configurations are required. We present a method based on machine learning to evaluate a key quantity entering the definition of absorption spectra: the dielectric screening. We show that our approach yields a model for the screening that is transferable between multiple configurations sampled during first principles molecular dynamics simulations; hence it leads to a substantial improvement in the efficiency of calculations of finite temperature spectra. We obtained computational gains of one to two orders of magnitude for systems with 50 to 500 atoms, including liquids, solids, nanostructures, and solid/liquid interfaces. Importantly, the models of dielectric screening derived here may be used not only in the solution of the BSE but also in developing functionals for time-dependent density functional theory (TDDFT) calculations of homogeneous and heterogeneous systems. Overall, our work provides a strategy to combine machine learning with electronic structure calculations to accelerate first principles simulations of excited-state properties.

OpenMolcas: From Source Code to Insight.

In this Article we describe the OpenMolcas environment and invite the computational chemistry community to collaborate. The open-source project already includes a large number of new developments realized during the transition from the commercial MOLCAS product to the open-source platform. The paper initially describes the technical details of the new software development platform. This is followed by brief presentations of many new methods, implementations, and features of the OpenMolcas program suite. These developments include novel wave function methods such as stochastic complete active space self-consistent field, density matrix renormalization group (DMRG) methods, and hybrid multiconfigurational wave function and density functional theory models. Some of these implementations include an array of additional options and functionalities. The paper proceeds and describes developments related to explorations of potential energy surfaces. Here we present methods for the optimization of conical intersections, the simulation of adiabatic and nonadiabatic molecular dynamics, and interfaces to tools for semiclassical and quantum mechanical nuclear dynamics. Furthermore, the Article describes features unique to simulations of spectroscopic and magnetic phenomena such as the exact semiclassical description of the interaction between light and matter, various X-ray processes, magnetic circular dichroism, and properties. Finally, the paper describes a number of built-in and add-on features to support the OpenMolcas platform with postcalculation analysis and visualization, a multiscale simulation option using frozen-density embedding theory, and new electronic and muonic basis sets.

Nature of the 11Bu and 21Ag Excited States of Butadiene and the Goldilocks Principle of Basis Set Diffuseness.

Butadiene, being the simplest conjugated organic molecule, has been studied extensively by experiments and various high-level theoretical methods. Previous studies conclude that the complete active space (CAS) self-consistent field (SCF) method was unable to obtain the correct 11Bu and 21Ag state ordering and that it introduces artificial valence-Rydberg mixing into the 11Bu state because of the lack of external correlation. Basis sets and initial guesses specifically constructed for this problem were able to improve the vertical excitation energy of the 11Bu state but did not resolve the controversy of the nature of the 11Bu state. In the present work, we demonstrate that, using standard intermediately diffuse basis sets such as jul-cc-pVTZ and ma-TZVP, CASSCF is able to obtain the 11Bu and 21Ag states of trans-butadiene with much improved vertical excitation energy and reasonable wave function characteristics, and it provides a good reference wave function (capable of yielding quantitative excitation energies for excited states) for three methods that treat electron correlation in different ways, namely, multiconfiguration pair-density functional theory (MC-PDFT), CAS second-order perturbation theory (PT2), and multistate (MS) CASPT2. We demonstrate that a combined analysis of the orbital second moments, state second moments, and MC-PDFT energy components is a very useful approach in determining excited-state characteristics and assigning states, and we show that basis sets without diffuse functions or with very diffuse basis functions are not stable or accurate in predicting the excited states of butadiene. We show that the 21Ag state is valence-like and has an atypical double/single excitation character and the 11Bu state has a certain degree of Rydberg character that is not artificial, settling the decades of controversy of the characters of these states.

State-Interaction Pair-Density Functional Theory Can Accurately Describe a Spiro Mixed Valence Compound.

Mixed-valence compounds with strong couplings between electronic states constitute one of the most challenging types of multireference systems for electronic structure theory. Previous work on a model mixed-valence compound, the 2,2',6,6'-tetrahydro-4 H,4' H-5,5'-spirobi[cyclopenta[ c]pyrrole] cation, showed that multireference perturbation theory (MRPT) can give a physical energy surface for the mixed-valence compound only by going to the third order or by using a scheme involving averaging orbital energies in a way specific to mixed-valence systems. In this study, we show that second-order MRPT methods (CASPT2, MS-CASPT2, and XMS-CASPT2) can give good results by calculating the Fock operator for the zeroth-order Hamiltonian using the state-averaged density matrix. We also show that state-interaction pair-density functional theory (SI-PDFT) is free from the unphysical behavior of previously tested second-order MRPT methods for this prototype mixed-valence compound near the avoided crossing. This is very encouraging because of the much lower cost in applying SI-PDFT to large or complex systems.

Assessment of MC-PDFT Excitation Energies for a Set of QM/MM Models of Rhodopsins.

A methodology for the automatic production of quantum mechanical/molecular mechanical (QM/MM) models of retinal-binding rhodopsin proteins and subsequent prediction of their spectroscopic properties has been proposed recently by some of the authors. The technology employed for the evaluation of the excitation energies is called Automatic Rhodopsin Modeling (ARM), and it involves the use of the complete active space self-consistent field (CASSCF) method followed by a multiconfiguration second-order perturbation theory (in particular, CASPT2) calculation of external correlation energies. Although it was shown that ARM is capable of successfully reproducing and predicting spectroscopic property trends in chromophore-embedding protein sets, practical applications of such technology are limited by the high computational costs of the multiconfiguration perturbation theory calculations. In the present work we benchmark the more affordable multiconfiguration pair-density functional theory (MC-PDFT) method whose accuracy has been recently validated for retinal chromophores in the gas phase, indicating that MC-PDFT could potentially be used to analyze large (e.g., few hundreds) sets of rhodopsin proteins. Here, we test this theory for a set of rhodopsin QM/MM models whose experimental absorption maxima (λ a max) have been measured. The results indicate that MC-PDFT may be employed to calculate λ a max values for this important class of photoresponsive proteins.

Extended Hamiltonian molecular dynamics: semiclassical trajectories with improved maintenance of zero point energy.

It is well known that classical trajectories, even if they are initiated with zero point energy (ZPE) in each mode (trajectories initiated this way are commonly called quasiclassical trajectories), do not maintain ZPE in the final states. The energy of high-frequency modes will typically leak into low-frequency modes or relative translation of subsystems during the time evolution. This can lead to severe problems such as unphysical dissociation of a molecule, production of energetically disallowed reaction products, and unphysical product energy distributions. Here a new molecular dynamics method called extended Hamiltonian molecular dynamics (EHMD) is developed to improve the ZPE problem in classical molecular dynamics. In EHMD, two images of a trajectory are connected by one or more springs. The EHMD method is tested with the Henon-Heiles Hamiltonian in reduced and real units and with a Hamiltonian with quartic anharmonicity in real units, and the method is found to improve zero-point maintenance as intended.

Excitation spectra of retinal by multiconfiguration pair-density functional theory.

Retinal is the chromophore in proteins responsible for vision. The absorption maximum of retinal is sensitive to mutations of the protein. However, it is not easy to predict the absorption spectrum of retinal accurately, and questions remain even after intensive investigation. Retinal poses a challenge for Kohn-Sham density functional theory (KS-DFT) because of the charge transfer character in its excitations, and it poses a challenge for wave function theory because the large size of the molecule makes multiconfigurational perturbation theory methods expensive. In this study, we demonstrate that multiconfiguration pair-density functional theory (MC-PDFT) provides an efficient way to predict the vertical excitation energies of 11-Z retinal, and it reproduces the experimentally determined absorption band widths and peak positions better than complete active space second-order perturbation theory (CASPT2). The consistency between complete active space self-consistent field (CASSCF) and KS-DFT dipole moments is demonstrated to be a useful criterion in selecting the active space. We also found that the nature of the terminal groups and the conformations of retinal play a significant role in the absorption spectrum. By considering a thermal distribution of conformations, we predict an absorption spectrum of retinal that is consistent with the experimental gas-phase spectrum. The location of the absorption peak and the spectral broadening based on MC-PDFT calculations agree better with experiments than those of CASPT2.

Automatic Selection of an Active Space for Calculating Electronic Excitation Spectra by MS-CASPT2 or MC-PDFT.

Multireference methods such as multistate complete active space second-order perturbation theory (MS-CASPT2) and multiconfiguration pair-density functional theory (MC-PDFT) can be very accurate, but they have the disadvantage that they are not black-box methods, and finding a good active space for the reference wave function often requires nonsystematic procedures based on intimate knowledge of the system under study. In this paper, we propose a scheme called the ABC scheme, which is a three-step calculation using three parameters, for automatic selection (without looking at the orbitals and without using any knowledge of the specific system at hand) of the active space (including both the size of the active space and the orbitals in the active space) for MS-CASPT2 or MC-PDFT calculations, and we report tests of the scheme on the first five excitation energies for a set of ten doublet radicals. The results show that the ABC scheme is very successful for both MS-CASPT2 and MC-PDFT for all ten systems considered here.

Identifying multiple active conformations in the G protein-coupled receptor activation landscape using computational methods.

G protein-coupled receptors (GPCRs) are membrane proteins critical in cellular signaling, making them important targets for therapeutics. The activation of GPCRs is central to their function, requiring multiple conformations of the GPCRs in their activation landscape. To enable rational design of GPCR-targeting drugs, it is essential to obtain the ensemble of atomistic structures of GPCRs along their activation pathways. This is most challenging for structure determination experiments, making it valuable to develop reliable computational structure prediction methods. In particular, since the active-state conformations are higher in energy (less stable) than inactive-state conformations, they are difficult to stabilize. In addition, the computational methods are generally biased toward lowest energy structures by design and miss these high energy but functionally important conformations. To address this problem, we have developed a computationally efficient ActiveGEnSeMBLE method that systematically predicts multiple conformations that are likely in the GPCR activation landscape, including multiple active- and inactive-state conformations. ActiveGEnSeMBLE starts with a systematic coarse grid sampling of helix tilts/rotations (~13 trillion transmembrane domain conformations) and identifies multiple potential active-state energy wells, using the TM3-TM6 intracellular distance as a surrogate activation coordinate. These energy wells are then sampled locally using a finer grid in conformational space to find a locally minimized conformation in each energy well, which can be further relaxed using molecular dynamics (MD) simulations. This method, combining homology modeling, hierarchical complete conformational sampling, and nanosecond scale MD, provides one of the very few computational methods that predict multiple candidates for active-state conformations and is one of the most computationally affordable.

Conformational and Thermodynamic Landscape of GPCR Activation from Theory and Computation.

We present a hybrid computational methodology to predict multiple energetically accessible conformations for G protein-coupled receptors (GPCRs) that might play a role in binding to ligands and different signaling partners. To our knowledge, this method, termed ActiveGEnSeMBLE, enables the first quantitative energy profile for GPCR activation that is consistent with the qualitative profile deduced from experiments. ActiveGEnSeMBLE starts with a systematic coarse grid sampling of helix tilts/rotations (∼13 trillion transmembrane-domain conformations) and selects the conformational landscape based on energy. This profile identifies multiple potential active-state energy wells, with the TM3-TM6 intracellular distance as an approximate activation coordinate. These energy wells are then sampled locally using a finer grid to find locally minimized conformation in each energy well. We validate this strategy using the inactive and active experimental structures of ß2 adrenergic receptor (hß2AR) and M2 muscarinic acetylcholine receptor. Structures of membrane-embedded hß2AR along its activation coordinate are subjected to molecular-dynamics simulations for relaxation and interaction energy analysis to generate a quantitative energy landscape for hß2AR activation. This landscape reveals several metastable states along this coordinate, indicating that for hß2AR, the agonist alone is not enough to stabilize the active state and that the G protein is necessary, consistent with experimental observations. The method's application to somatostatin receptor SSTR5 (no experimental structure available) shows that to predict an active conformation it is better to start from an inactive structure template based on a close homolog than to start from an active template based on a distant homolog. The energy landscape for hSSTR5 activation is consistent with hß2AR in the role of the G protein. These results demonstrate the utility of the ActiveGEnSeMBLE method for predicting multiple conformations along the pathways for activating GPCRs and the corresponding energy landscapes, thereby providing detailed structural insights into the initial molecular events of GPCR function that are not easily accessible by experiments.

Computational predictions of corroles as a class of Hsp90 inhibitors.

Corroles have been shown experimentally to cause cell cycle arrest, and there is some evidence that this might be attributed to an inhibitory effect of corroles on Heat shock protein 90 (Hsp90), which is known to play a vital role in cancer cell proliferation. In this study, we used molecular dynamics to examine the interaction of gallium corroles with Hsp90, and found that they can bind preferentially to the ATP-binding N-terminal site. We also found that structural variations of the corrole ring can influence the binding energies and affinities of the corrole to Hsp90. We predict that both the bis-carboxylated corrole (4-Ga) and a proposed 3,17-bis-sulfonated corrole (7-Ga) are promising alternatives to Ga(III) 5,10,15-tris(pentafluorophenyl)-2,17-bis(sulfonic acid)-corrole (1-Ga) as anti-cancer agents.

The predicted ensemble of low-energy conformations of human somatostatin receptor subtype 5 and the binding of antagonists.

Human somatostatin receptor subtype 5 (hSSTR5) regulates cell proliferation and hormone secretion. However, the identification of effective therapeutic small-molecule ligands is impeded because experimental structures are not available for any SSTR subtypes. Here, we predict the ensemble of low-energy 3D structures of hSSTR5 using a modified GPCR Ensemble of Structures in Membrane BiLayer Environment (GEnSeMBLE) complete sampling computational method. We find that this conformational ensemble displays most interhelical interactions conserved in class A G protein-coupled receptors (GPCRs) plus seven additional interactions (e.g., Y2.43-D3.49, T3.38-S4.53, K5.64-Y3.51) likely conserved among SSTRs. We then predicted the binding sites for a series of five known antagonists, leading to predicted binding energies consistent with experimental results reported in the literature. Molecular dynamics (MD) simulation of 50 ns in explicit water and lipid retained the predicted ligand-bound structure and formed new interaction patterns (e.g. R3.50-T6.34) consistent with the inactive µ-opioid receptor X-ray structure. We suggest more than six mutations for experimental validation of our prediction. The final predicted receptor conformations and antagonist binding sites provide valuable insights for designing new small-molecule drugs targeting SSTRs.

Highly efficient and diastereoselective gold(I)-catalyzed synthesis of tertiary amines from secondary amines and alkynes: substrate scope and mechanistic insights.

An efficient method for the synthesis of tertiary amines through a gold(I)-catalyzed tandem reaction of alkynes with secondary amines has been developed. In the presence of ethyl Hantzsch ester and [{(tBu)(2)(o-biphenyl)P}AuCl]/AgBF(4) (2 mol %), a variety of secondary amines bearing electron-deficient and electron-rich substituents and a wide range of alkynes, including terminal and internal aryl alkynes, aliphatic alkynes, and electron-deficient alkynes, underwent a tandem reaction to afford the corresponding tertiary amines in up to 99 % yield. For indolines bearing a preexisting chiral center, their reactions with alkynes in the presence of ethyl Hantzsch ester catalyzed by [{(tBu)(2)(o-biphenyl)P}AuCl]/AgBF(4) (2 mol %) afforded tertiary amines in excellent yields and with good to excellent diastereoselectivity. All of these organic transformations can be conducted as a one-pot reaction from simple and readily available starting materials without the need of isolation of air/moisture-sensitive enamine intermediates, and under mild reaction conditions (mostly room temperature and mild reducing agents). Mechanistic studies by NMR spectroscopy, ESI-MS, isotope labeling studies, and DFT calculations on this gold(I)-catalyzed tandem reaction reveal that the first step involving a monomeric cationic gold(I)-alkyne intermediate is more likely than a gold(I)-amine intermediate, a three-coordinate gold(I) intermediate, or a dinuclear gold(I)-alkyne intermediate. These studies also support the proposed reaction pathway, which involves a gold(I)-coordinated enamine complex as a key intermediate for the subsequent transfer hydrogenation with a hydride source, and reveal the intrinsic stereospecific nature of these transformations observed in the experiments.

Electronic structures of Group 9 metallocorroles with axial ammines.

The electronic structures of metallocorroles (tpfc)M(NH(3))(2) and (tfc)M(NH(3))(2) (tpfc is the trianion of 5,10,15-(tris)pentafluorophenylcorrole, tfc is the trianion of 5,10,15-trifluorocorrole, and M = Co, Rh, Ir) have been computed using first principles quantum mechanics [B3LYP flavor of Density Functional Theory (DFT) with Poisson-Boltzmann continuum solvation]. The geometry was optimized for both the neutral systems (formal M(III) oxidation state) and the one-electron oxidized systems (formally M(IV)). As expected, the M(III) systems have a closed shell d(6) configuration; for all three metals, the one-electron oxidation was calculated to occur from a ligand-based orbital (highest occupied molecular orbital (HOMO) of B(1) symmetry). The ground state of the formal M(IV) system has M(III)-Cπ character, indicating that the metal remains d(6), with the hole in the corrole π system. As a result the calculated M(IV/III) reduction potentials are quite similar (0.64, 0.67, and 0.56 V vs SCE for M = Ir, Rh and Co, respectively), whereas the differences would have been large for purely metal-based oxidations. Vertically excited states with substantial metal character are well separated from the ground state in one-electron-oxidized cobalt (0.27 eV) and rhodium (0.24 eV) corroles, but become closer in energy in the iridium (0.15 eV) analogues. The exact splittings depend on the chosen functional and basis set combination and vary by ~0.1 eV.