From NWChem to NWChemEx: Evolving with the Computational Chemistry Landscape.

Since the advent of the first computers, chemists have been at the forefront of using computers to understand and solve complex chemical problems. As the hardware and software have evolved, so have the theoretical and computational chemistry methods and algorithms. Parallel computers clearly changed the common computing paradigm in the late 1970s and 80s, and the field has again seen a paradigm shift with the advent of graphical processing units. This review explores the challenges and some of the solutions in transforming software from the terascale to the petascale and now to the upcoming exascale computers. While discussing the field in general, NWChem and its redesign, NWChemEx, will be highlighted as one of the early codesign projects to take advantage of massively parallel computers and emerging software standards to enable large scientific challenges to be tackled.

A Nonorthogonal Configuration Interaction Approach to Singlet Fission in Perylenediimide Compounds.

Perylenediimide molecules constitute a family of chromophores that undergo singlet fission, a process in which an excited singlet state converts into lower energy triplets on two neighboring molecules, potentially increasing the efficiency of organic solar cells. Here, the nonorthogonal configuration interaction method is applied to study the effect of the different crystal packing of various perylenediimide derivatives on the relative energies of the singlet and triplet states, the intermolecular electronic couplings, and the relative rates for singlet fission. The analysis of the wave functions and electronic couplings reveals that charge transfer states play an important role in the singlet fission mechanism. Dimer conformations where the PDI molecules are at large displacements along the long axis and short on the short axis are posed as the most favorable for singlet fission. The role of the substituent at the imide group has been inspected concluding that, although it has no effect in the energies, for some conformations it significantly influences the electronic couplings, and therefore, replacing this substituent with hydrogen may introduce artifacts in the computational modeling of the PDI molecules.

On the role of dynamic electron correlation in non-orthogonal configuration interaction with fragments.

Two different approaches have been implemented to include the effect of dynamic electron correlation in the Non-Orthogonal Configuration Interaction for Fragments (NOCI-F) method. The first is based on shifting the diagonal matrix elements of the NOCI matrix, while the second incorporates the dynamic correlation explicitly in the fragment wave functions used to construct the many-electron basis functions of the NOCI. The two approaches are illustrated for the calculation of the electronic coupling relevant in singlet fission and the coupling of spin moments in organic radicals. Comparison of the calculated diabatic couplings, the NOCI energies and wave functions shows that dynamic electron correlation is not only efficiently but also effectively incorporated by the shifting approach and can largely affect the coupling between electronic states. Also, it brings the NOCI coupling of the spin moments in close agreement with benchmark calculations.

GronOR: Massively parallel and GPU-accelerated non-orthogonal configuration interaction for large molecular systems.

GronOR is a program package for non-orthogonal configuration interaction calculations for an electronic wave function built in terms of anti-symmetrized products of multi-configuration molecular fragment wave functions. The two-electron integrals that have to be processed may be expressed in terms of atomic orbitals or in terms of an orbital basis determined from the molecular orbitals of the fragments. The code has been specifically designed for execution on distributed memory massively parallel and Graphics Processing Unit (GPU)-accelerated computer architectures, using an MPI+OpenACC/OpenMP programming approach. The task-based execution model used in the implementation allows for linear scaling with the number of nodes on the largest pre-exascale architectures available, provides hardware fault resiliency, and enables effective execution on systems with distinct central processing unit-only and GPU-accelerated partitions. The code interfaces with existing multi-configuration electronic structure codes that provide optimized molecular fragment orbitals, configuration interaction coefficients, and the required integrals. Algorithm and implementation details, parallel and accelerated performance benchmarks, and an analysis of the sensitivity of the accuracy of results and computational performance to thresholds used in the calculations are presented.

NWChem: Past, present, and future.

Specialized computational chemistry packages have permanently reshaped the landscape of chemical and materials science by providing tools to support and guide experimental efforts and for the prediction of atomistic and electronic properties. In this regard, electronic structure packages have played a special role by using first-principle-driven methodologies to model complex chemical and materials processes. Over the past few decades, the rapid development of computing technologies and the tremendous increase in computational power have offered a unique chance to study complex transformations using sophisticated and predictive many-body techniques that describe correlated behavior of electrons in molecular and condensed phase systems at different levels of theory. In enabling these simulations, novel parallel algorithms have been able to take advantage of computational resources to address the polynomial scaling of electronic structure methods. In this paper, we briefly review the NWChem computational chemistry suite, including its history, design principles, parallel tools, current capabilities, outreach, and outlook.

Retention of conformational entropy upon calmodulin binding to target peptides is driven by transient salt bridges.

Calmodulin (CaM) is a highly flexible calcium-binding protein that mediates signal transduction through an ability to differentially bind to highly variable binding sequences in target proteins. To identify how binding affects CaM motions, and its relationship to conformational entropy and target peptide sequence, we have employed fully atomistic, explicit solvent molecular dynamics simulations of unbound CaM and CaM bound to five different target peptides. The calculated CaM conformational binding entropies correlate with experimentally derived conformational entropies with a correlation coefficient R(2) of 0.95. Selected side-chain interactions with target peptides restrain interhelical loop motions, acting to tune the conformational entropy of the bound complex via widely distributed CaM motions. In the complex with the most conformational entropy retention (CaM in complex with the neuronal nitric oxide synthase binding sequence), Lys-148 at the C-terminus of CaM forms transient salt bridges alternating between Glu side chains in the N-domain, the central linker, and the binding target. Additional analyses of CaM structures, fluctuations, and CaM-target interactions illuminate the interplay between electrostatic, side chain, and backbone properties in the ability of CaM to recognize and discriminate against targets by tuning its conformational entropy, and suggest a need to consider conformational dynamics in optimizing binding affinities.

GronOR: Scalable and Accelerated Nonorthogonal Configuration Interaction for Molecular Fragment Wave Functions.

GronOR is a program package for nonorthogonal configuration interaction calculations. Electronic wave functions are constructed in terms of antisymmetrized products of multiconfiguration molecular fragment wave functions. The computational complexity of the nonorthogonal methodologies implemented in GronOR applied to large molecular assemblies requires a design that takes full advantage of massively parallel supercomputer architectures and accelerator technologies. This work describes the implementation strategy and resulting performance characteristics. In addition to parallelization and acceleration, the software development strategy includes aspects of fault resiliency and heterogeneous computing. The program was designed for large-scale supercomputers but also runs effectively on small clusters and workstations for small molecular systems. GronOR is available as open source to the scientific community.

Characterization of the outer membrane protein OprF of Pseudomonas aeruginosa in a lipopolysaccharide membrane by computer simulation.

The N-terminal domain of outer membrane protein OprF of Pseudomonas aeruginosa forms a membrane spanning eight-stranded antiparallel beta-barrel domain that folds into a membrane channel with low conductance. The structure of this protein has been modeled after the crystal structure of the homologous protein OmpA of Escherichia coli. A number of molecular dynamics simulations have been carried out for the homology modeled structure of OprF in an explicit molecular model for the rough lipopolysaccharide (LPS) outer membrane of P. aeruginosa. The structural stability of the outer membrane model as a result of the strong electrostatic interactions compared with simple lipid bilayers is restricting both the conformational flexibility and the lateral diffusion of the porin in the membrane. Constricting side-chain interactions within the pore are similar to those found in reported simulations of the protein in a solvated lipid bilayer membrane. Because of the strong interactions between the loop regions of OprF and functional groups in the saccharide core of the LPS, the entrance to the channel from the extracellular space is widened compared with the lipid bilayer simulations in which the loops are extruding in the solvent. The specific electrostatic signature of the LPS membrane, which results in a net intrinsic dipole across the membrane, is found to be altered by the presence of OprF, resulting in a small electrically positive patch at the position of the channel.

Open "back door" in a molecular dynamics simulation of acetylcholinesterase.

The enzyme acetylcholinesterase generates a strong electrostatic field that can attract the cationic substrate acetylcholine to the active site. However, the long and narrow active site gorge seems inconsistent with the enzyme's high catalytic rate. A molecular dynamics simulation of acetylcholinesterase in water reveals the transient opening of a short channel, large enough to pass a water molecule, through a thin wall of the active site near tryptophan-84. This simulation suggests that substrate, products, or solvent could move through this "back door," in addition to the entrance revealed by the crystallographic structure. Electrostatic calculations show a strong field at the back door, oriented to attract the substrate and the reaction product choline and to repel the other reaction product, acetate. Analysis of the open back door conformation suggests a mutation that could seal the back door and thus test the hypothesis that thermal motion of this enzyme may open multiple routes of access to its active site.

Computer simulation of uranyl uptake by the rough lipopolysaccharide membrane of Pseudomonas aeruginosa.

Heavy metal environmental contaminants cannot be destroyed but require containment, preferably in concentrated form, in a solid or immobile form for recycling or final disposal. Microorganisms are able to take up and deposit high levels of contaminant metals, including radioactive metals such as uranium and plutonium, into their cell wall. Consequently, these microbial systems are of great interest as the basis for potential environmental bioremediation technologies. The outer membranes of Gram-negative microbes are highly nonsymmetric and exhibit a significant electrostatic potential gradient across the membrane. This gradient has a significant effect on the uptake and transport of charged and dipolar compounds. However, the effectiveness of microbial systems for environmental remediation will depend strongly on specific properties that determine the uptake of targeted contaminants by a particular cell wall. To aid in the design of microbial remediation technologies, knowledge of the factors that determine the affinity of a particular bacterial outer membrane for the most common ionic species found in contaminated soils and groundwater is of great importance. Using our previously developed model for the lipopolysaccharide (LPS) membrane of Pseudomonas aeruginosa, this work presents the potentials of mean force as the estimate of the free energy profile for uptake of sodium, calcium, chloride, uranyl ions, and a water molecule by the bacterial LPS membrane. A compatible classical parameter set for uranyl has been developed and validated. Results show that the uptake of uranyl is energetically a favorable process relative to the other ions studied. At neutral pH, this nuclide is shown to be retained on the surface of the LPS membrane through chelation with the carboxyl and hydroxyl groups located in the outer core.

Electronic coupling between heme electron-transfer centers and its decay with distance depends strongly on relative orientation.

A method for calculating the electron-transfer matrix element V(RP) using density functional theory Kohn-Sham orbitals is presented and applied to heme dimers of varying relative orientation. The electronic coupling decays with increased iron separation according to V(RP) = V(0)(RP)exp(-beta r/2) with a distance dependence parameter beta approximately 2 A(-1) for hemes with parallel porphyrins and either 1.1 or 4.0 A(-1) when the porphyrin planes are perpendicular, depending on the alignment of the iron d(pi) orbital. These findings are used to interpret the observed orientation of the hemes in tetraheme redox proteins such as Flavocytochrome c(3) fumarate reductase (Ifc(3), PDB code 1QJD) of Shewanella frigidimarina, another flavocytochrome from the same bacterium (Fcc(3), 1E39) and a small tetraheme cytochrome of Shewanella oneidensis strain MR1 (1M1P). Our results show that shifting and rotating the hemes controls the adiabaticity of the three electron hopping steps.

Molecular docking of balanol to dynamics snapshots of protein kinase A.

Even if the structure of a receptor has been determined experimentally, it may not be a conformation to which a ligand would bind when induced fit effects are significant. Molecular docking using such a receptor structure may thus fail to recognize a ligand to which the receptor can bind with reasonable affinity. Here, we examine one way to alleviate this problem by using an ensemble of receptor conformations generated from a molecular dynamics simulation for molecular docking. Two molecular dynamics simulations were conducted to generate snapshots for protein kinase A: one with the ligand bound, the other without. The ligand, balanol, was then docked to conformations of the receptors presented by these trajectories. The Lamarckian genetic algorithm in Autodock [Goodsell et al. J Mol Recognit 1996;9(1):1-5; Morris et al. J Comput Chem 1998;19(14):1639-1662] was used in the docking. Three ligand models were used: rigid, flexible, and flexible with torsional potentials. When the snapshots were taken from the molecular dynamics simulation of the protein-ligand complex, the correct docking structure could be recovered easily by the docking algorithm in all cases. This was an easier case for challenging the docking algorithm because, by using the structure of the protein in a protein-ligand complex, one essentially assumed that the protein already had a pocket to which the ligand can fit well. However, when the snapshots were taken from the ligand-free protein simulation, which is more useful for a practical application when the structure of the protein-ligand complex is not known, several clusters of structures were found. Of the 10 docking runs for each snapshot, at least one structure was close to the correctly docked structure when the flexible-ligand models were used. We found that a useful way to identify the correctly docked structure was to locate the structure that appeared most frequently as the lowest energy structure in the docking experiments to different snapshots.

Characterization of the active site of DNA polymerase beta by molecular dynamics and quantum chemical calculation.

It is well established that the fully formed polymerase active site of the DNA repair enzyme, polymerase beta (pol beta), including two bound Mg2+ cations and the nucleoside triphosphate (dNTP) substrate, exists at only one point in the catalytic cycle just prior to the chemical nucleotidyl transfer step. The structure of the active conformation has been the subject of much interest as it relates to the mechanism of the chemical step and also to the question of fidelity assurance. Although crystal structures of ternary pol beta-(primer-template) DNA-dNTP complexes have provided the main structural features of the active site, they are necessarily incomplete due to intentional alterations (e.g., removal of the 3'OH groups from primer and substrate) needed to obtain a structure from midcycle. Working from the crystal structure closest to the fully formed active site [Protein Data Bank (PDB) code: 1bpy], two molecular dynamics (MD) simulations of the solvated ternary complex were performed: one with the missing 3'OHs restored, via modeling, to the primer and substrate, and the other without restoration of the 3'OHs. The results of the simulations, together with ab initio optimizations on simplified active-site models, indicate that the missing primer 3'OH in the crystal structure is responsible for a significant perturbation in the coordination sphere of the catalytic cation and allow us to suggest several corrections and additions to the active-site structure as observed by crystallography. In addition, the calculations help to resolve questions raised regarding the protonation states of coordinating ligands.

Theoretical study of inhibition of adenosine deaminase by (8R)-coformycin and (8R)-deoxycoformycin.

Molecular dynamics and free energy simulations were performed to examine the binding of (8R)-deoxycoformycin and (8R)-coformycin to adenosine deaminase. The two inhibitors differ only at the 2' position of the sugar ring; the sugar moiety of conformycin is ribose, while it is deoxyribose for deoxycoformycin. The 100 ps molecular dynamics trajectories reveal that Asp 19 and His 17 interact strongly with the 5' hydroxyl group of the sugar moiety of both inhibitors and appear to play an important role in binding the sugar. The 2' and 3' groups of the sugars are near the protein-water interface and can be stabilized by either protein residues or water. The flexibility of the residues at the opening of the active site helps to explain the modest difference in binding of the two inhibitors and how substrates/inhibitors can enter an otherwise inaccessible binding site.

Force-Field Development and Molecular Dynamics of [NiFe] Hydrogenase.

Classical molecular force-field parameters describing the structure and motion of metal clusters in [NiFe] hydrogenase enzymes can be used to compare the dynamics and thermodynamics of [NiFe] under different oxidation, protonation, and ligation circumstances. Using density functional theory (DFT) calculations of small model clusters representative of the active site and the proximal, medial, and distal Fe/S metal centers and their attached protein side chains, we have calculated classical force-field parameters for [NiFe] in reduced and oxidized states, including internal coordinates, force constants, and atom-centered charges. Derived force constants revealed that cysteinate ligands bound to the metal ions are more flexible in the Ni-B active site, which has a bridging hydroxide ligand, than in the Ni-C active site, which has a bridging hydride. Ten nanosecond all-atom, explicit-solvent MD simulations of [NiFe] hydrogenase in oxidized and reduced catalytic states established the stability of the derived force-field parameters in terms of Cα and metal cluster fluctuations. Average active site structures from the protein MD simulations are consistent with [NiFe] structures from the Protein Data Bank, suggesting that the derived force-field parameters are transferrable to other hydrogenases beyond the structure used for testing. A comparison of experimental H2-production rates demonstrated a relationship between cysteinate side chain rotation and activity, justifying the use of a fully dynamic model of [NiFe] metal cluster motion.

Molecular basis of the structural stability of a Top7-based scaffold at extreme pH and temperature conditions.

The development of stable biomolecular scaffolds that can tolerate environmental extremes has considerable potential for industrial and defense-related applications. However, most natural proteins are not sufficiently stable to withstand non-physiological conditions. We have recently engineered the de novo designed Top7 protein to specifically recognize the glycoprotein CD4 by insertion of an eight-residue loop. The engineered variant exhibited remarkable stability under chemical and thermal denaturation conditions. In the present study, far-UV CD spectroscopy and explicit-solvent MD simulations are used to investigate the structural stability of Top7 and the engineered variant under extreme conditions of temperature and pH. Circular dichroism measurements suggest that the engineered variant Top7(CB1), like Top7, retains its structure at high temperatures. Changes in CD spectra suggest that there are minor structural rearrangements between neutral and acidic environments for both proteins but that these do not make the proteins less stable at high temperatures. The anti-parallel beta-sheet is well conserved within the timescale simulated whereas there is a decrease of helical content when low pH and high-temperature conditions are combined. Concerted alanine mutations along the alpha-helices of the engineered Top7 variant did not revert this trend when at pH 2 and 400K. The structural resilience of the anti-parallel beta-sheet suggests that the protein scaffold can accommodate varying sequences. The robustness of the Top7 scaffold under extreme conditions of pH and temperature and its amenability to production in inexpensive bacterial expression systems reveal great potential for novel biotechnological applications.

Interaction between the CBM of Cel9A from Thermobifida fusca and cellulose fibers.

Molecular docking and molecular dynamics (MD) simulations were used to investigate the binding of a cellodextrin chain in a crystal-like conformation to the carbohydrate-binding module (CBM) of Cel9A from Thermobifida fusca. The fiber was found to bind to the CBM in a single and well-defined configuration in-line with the catalytic cleft, supporting the hypothesis that this CBM plays a role in the catalysis by feeding the catalytic domain (CD) with a polysaccharide chain. The results also expand the current known list of residues involved in the binding. The polysaccharide-protein attachment is shown to be mediated by five amine/amide-containing residues. E478 and E559 were found not to interact directly with the sugar chain; instead they seem to be responsible to stabilize the binding motif via hydrogen bonds.

Characterization of Chitin and Chitosan Molecular Structure in Aqueous Solution.

Molecular dynamics simulations have been used to characterize the structure of single chitin and chitosan chains in aqueous solutions. Chitin chains, whether isolated or in the form of a ß-chitin nanoparticle, adopt the 2-fold helix with ϕ and φ values similar to its crystalline state. In solution, the intramolecular hydrogen bond HO3(n)···O5(n+1) responsible for the 2-fold helical motif in these polysaccharides is stabilized by hydrogen bonds with water molecules in a well-defined orientation. On the other hand, chitosan can adopt five distinct helical motifs, and its conformational equilibrium is highly dependent on pH. The hydrogen bond pattern and solvation around the O3 atom of insoluble chitosan (basic pH) are nearly identical to these quantities in chitin. Our findings suggest that the solubility and conformation of these polysaccharides are related to the stability of the intrachain HO3(n)···O5(n+1) hydrogen bond, which is affected by the water exchange around the O3-HO3 hydroxyl group.

Molecular Dynamics of Organophosphorous Hydrolases Bound to the Nerve Agent Soman.

The organophosphorous hydrolase (OPH) from Pseudomonas diminuta is capable of degrading extremely toxic organophosphorous compounds with a high catalytic turnover and broad substrate specificity. Although the natural substrate for OPH is unknown, its triple-mutant H254G/H257W/L303T exhibits a 3 order of magnitude increase in catalytic efficiency and modified stereospecificity toward the most toxic SpSc enantiomer of soman. Molecular dynamics simulations and binding free-energy calculations have been undertaken for the wild-type and triple-mutant H254G/H257W/L303T enzymes bound to the SpSc-soman enantiomer. Comparison of the simulations indicates that substrate binding induces conformational changes of the loops near the active site. The coordination of the zinc cations in the active site of OPH differs between the free enzyme and the complexes. This suggests that the active site of OPH can accommodate several catalytically active coordination geometries, consistent with the fact that the enzymatic activity of the wild-type OPH can be enhanced by alterations to the metal content of the enzyme. It is also argued that the enhanced efficiency of the triple mutant is determined by enzyme-transition-state complementarity. These results provide a qualitative, molecular-level explanation for the 3 order of magnitude increase in catalytic efficiency of the triple-mutant toward SpSc-soman.