• High performance computing (HPC) and large scale, data intensive scientific computing leveraging machine learning with focus on (bio)chemistry and nanotechnology.
• Highly parallel algorithms and software implementations for predictive computational science through quantum mechanical (QM) and classical atomistic numerical simulations on well established and emerging HPC platforms including graphics processing units (GPUs).
• Implementations mainly into the open-source GPU accelerated quantum chemistry software QUICK, the ADF software package for density functional theory (DFT), and the AMBER software package for classical MM and mixed QM/MM molecular dynamics (MD) simulations. Both ADF and AMBER are used by a large number of academic and industrial researchers world wide while QUICK is a more recent code.
• Education and training of the next generation of scientists in software engineering and numerical simulation methods, quantum chemistry, molecular dynamics, and data science.

We are grateful to our sponsors for financial support. Our research has been supported by federal and private organizations:

U.S. Department of Energy	National Institutes of Health	National Science Foundation
Intel	Nvidia	NSF Supercomputer access

Graphical processing units (GPUs) offer a tremendous amount of computing power in a compact package (see graphics to the right). Both the peak floating point operations per second and the memory bandwidth of NVIDIA GPUs compare favorably to Intel CPUs. However, this comes at the cost of reduced flexibility and increased programming complexity as compared to CPUs.

I am involved in an effort to port the all-atom classical MD engine PMEMD of the AMBER biomolecular simulation package to GPUs [1-3]. The AMBER implementation uses CUDA which is a relatively simple extension of the standard C programming language that allows one to code in an inherently parallel fashion and perform all necessary operations to access and manipulate data on a GPU device.

Peak floating point operations per second (Flop/s; left) and memory bandwidth (right) for Intel CPUs and NVIDIA GPUs. (Taken from our publication [1] on GPU accelerated MD) The data is from 2012 but the performance differential looks similar as of 2018.

We have developed a new precision model in which the contributions to the forces are calculated in single precision floating point arithmetics but accumulated in accumulated in 64-bit fixed-point integers (SPFP). The numerical noise due to rounding errors in a pure single precision implementation is sufficiently so large that it leads to accumulation of errors that can result in unphysical trajectories for long timescale simulations. The numerical results of the mixed-precision SPFP model are comparable to a full double precision model and the reference double precision CPU implementation but at significantly reduced computational cost.

The AMBER GPU implementatin provides performance for MD simulations on a single desktop machine that is on par with, and in some cases exceeds, that of traditional supercomputers (see table on the left). Routine microsecond MD simulations are now possible for systems with several hundred thousands of atoms under periodic boundary conditions. With GPUs becoming ubiquitous in workstations and also as accelerators in HPC platforms, the impact of this implementation on the field of molecular dynamics is broad and transformative.

QM/MM approaches in AMBER

Hybrid quantum mechanical and molecular mechanical (QM/MM) approaches are used extensively to study local electronic events in large molecular systems with a diverse area of applications ranging from enzymatic catalysis to properties of materials systems. We are continuously developing new QM/MM approaches. Among others we are working on a GPU accelerated software library for QM and QM/MM simulations.

We have developed an implementation that couples the MD software AMBER to electronic structure software packages including ADF, GAMESS-US, Gaussian, Orca, NWChem and TeraChem [1,2]. With this interface, ab initio wave function theory and DFT methods become available for QM/MM MD simulations with AMBER.

This represents a set of widely used programs, both commercial and freely available, each with its own strengths for different electronic structure methods and computing platforms ranging from desktop workstations to supercomputers. In the case of TeraChem this also includes accelerator hardware in form of graphics processing units (GPUs). Data exchange between AMBER and the external QM software is currently implemented by means of files and system calls or the message passing interface (MPI) standard.

This interface was used for example to study the absorption spectrum of the photoactive yellow protein (PYP) (see picture on right). It was shown that a large QM region (hundreds of atoms) is required to achieve a converged spectrum [2].

[1] A. W. Götz, M. A. Clark, R. C. Walker, J. Comput. Chem. 35, 95-108 (2014). DOI: 10.1002/jcc.23444

[2] C. M. Isborn , A. W. Götz , M. A. Clark , R. C. Walker, T. J. Martínez, J. Chem. Theory Comput. 8, 5092-5106 (2012). DOI: 10.1021/ct3006826

Adaptive QM/MM approaches

I have implemented a parallelized adaptive QM/MM approach into the AMBER MD software package [1]. It is based on a method developed by Rosa Bulo [2].

Adaptive QM/MM enables the QM treatment of an active region of a solvated molecular system including the solvent molecules in its vicinity (see left). The environment is handled with MM to reach high computational efficiency. During a simulation, solvent molecules can move via a transition region between the active and environment regions.

I was awarded a TRO (Triton Research Opportunity) grant by SDSC to collaborate with Francesco Paesani and use my implementation to study processes of relevance for atmospheric chemistry. Together with Kyoyeon Park we have shown the potential of adaptive QM/MM for MD simulations of aqueous systems [3].

Adaptive QM/MM can be applied to lots of exciting problems, for example solvent effects on absorption spectra, binding free energies of ions or drug-like molecules to proteins or enzymes.

[1] A. W. Götz, K. Park, R. E. Bulo, F. Paesani, R. C. Walker, publication in preparation (2015).

[2] R. E. Bulo, B. Ensing, J. Sikkema, L. Visscher, J. Chem. Theory Comput. 5, 2212-2221 (2009).

[3] K. Park, A. W. Götz, R. C. Walker, F. Paesani,
J. Chem. Theory Comput. 8, 2868-2877 (2012). DOI: 10.1021/ct300331f

The Goetz group is part of the NSF Center for Aerosol Impacts on Chemistry of the Environment (CAICE), a multi-institution, interdisciplinary research center with headquarters at UC San Diego. The goal of CAICE is to understand fundamental processes regarding the chemistry of aerosol particles and their impact on the environment.

We develop highly accurate many-body potentials based on correlated electronic structure reference data and machine learning techniques for atmospherically relevant molecules and their interactions with aqueous aerosol particles containing ions and surfactants. We use these potentials for computer simulations of gas scattering and reactive uptake of atmospheric gases to sea spray mimics.

Molecular dynamics simulations with rare events sampling methods give access to structural, thermodynamic, and spectroscopic properties from small model clusters to condensed phase systems with explicit liquid-gas interfaces. We also employ QM/MM methods to study chemical reaction dynamics.

The figure to the right shows an early application with a simple Amber GAFF force field for the important pollutant dinitrogenpentoxide (N₂O₅) and the TIP4P water potential. We have determined the free energy change when N₂O₅ transfers from liquid water to the gas phase. We are currently developing a highly accurate potential that will allow a quantitative prediction of the solubility, a quantity that is hard to obtain experimentally.

[1] B. Hirshberg, E. Rossich molina, A. W. Götz, A. Hammerich, G. Nathanson, T. H. Bertram, M. A. Johnson, R. B. Gerber, Phys. Chem. Chem. Phys. 20, 17961-17976 (2018). DOI: 10.1039/c8cp03022g

Cytochrome c oxidases (CcOs) are redox driven proton pumps in the membranes of mitochondria and many aerobic bacteria. CcOs are the terminal oxidases (Complex IV) of the electron transport chain and drive ATP synthesis by producing the required electrochemical potential gradient across the membrane, in conjunction with Complexes I and III.

We employ density functional theory (DFT) calculations and classical and QM/MM molecular dynamics simulations to provide a detailed mechanistic understanding of the catalytic reaction pathways and mechanisms of proton pumping in CcO.

We work closely with the group of Prof. Louis Noodleman at Scripps Research. Based on DFT calculations we have developed a detailed reaction mechanism for the oxygen reduction to water in the dinucler reaction center (DNC) of CcO. [1] This is a chemically difficult but biologically important reaction.

Using molecular dynamics simulations, we have identified water exit pathways from the water pool above the DNC to the P-side of the membrane, which could also act as proton transfer pathways. [2] We are currently investigating plausible proton transport mechanisms.

[1] L. Noodleman, W. Han Du, J. Fee, A. W. Götz, R. C. Walker Inorg. Chem. 53, 6458-6472 (2014). DOI: 10.1021/ic500363h

[2] L. Yang, A. A. Skjevik, W.-G. Han Du, L. Noodleman, R. C. Walker, A. W. Götz, BBA, Bioenergetics. 1857, 1594-1606 (2016). DOI: 10.1016/j.bbabio.2016.06.005

Computer simulations are a powerful tool for the modeling of molecules and materials but are still limited by a tradeoff between the accuracy of the molecular models and the associated computational cost. The many-body expansion of the total energy can be effectively employed in combination with machine learning techniques to develop highly accurate models that are useful for practical applications.

Long-range and higher-order terms are represented by classical dispersion, permanent electrostatics, and induction. Machine learning techniques are then used to encode the highly complex quantum mechanical many-body interactions that arise when molecules interact at short range. A successful example of this approach is the MB-pol water model and its extension to atomic ions in water, see for example [1,2].

[1] S. K. Reddy, S. C. Straight, P. Bajaj, C. H. Pham, M. Riera, D. R. Moberg, M. A. Morales, C. Knight, A. W. Götz, F. Paesani J. Chem. Phys. 145, 194594 (2016). DOI: 10.1063/1.4967719

[2] M. Riera, N. Mardirossian, P. Bajaj, A. W. Götz, F. Paesani J. Chem. Phys. 147, 161715 (2017). DOI: 10.1063/1.4993213

We are developing an open-source software infrastructure for the automated generation of many-body potentials for generic molecules. The central building block is a database that stores all relevant properties and molecular configurations that are required to train machine learning models of 1-body monomer distortions and 2-body and 3-body interaction potentials. We exploit distributed and volunteer computing for massively parallel electronic structure calculations.

We recently have shown that different machine learning techniques can be used to represent the short-range interactions of the many-body expansion. [3] We currently work with permutationally invariant polynomials (PIPs) and Behler-Parinello type neural networks (BPNN), see Figure on the left. We work closely with Professors Francesco Paesani from the Chemistry Department and Sean Gao from Computer Science to improve these machine learning algorithms.

These alternative machine learning models are both highly accurate, which means that researchers can choose the algorithms that best map to the available hardware. Modern many-core processors, for instance, are well-suited to evaluate the complex expressions of the PIPs, while massively parallel graphics processing units (GPUs) perform exceptionally well for neural networks.

We also work on efficient implementations of these machine learning based many-body potentials. We employ C++ with OpenMP for many-core CPUs and C++ with CUDA and Nvidia's deep learning library for massively parallel GPUs.

[3] T. T. Nguyen, E. Székely, G. Imbalzano, J. Behler, G. Csányi, M. Ceriotti, A. W. Götz, F. Paesani, J. Chem. Phys. 148, 241725 (2018). DOI: 10.1063/1.5024577

We develop models and software to better understand at a molecular level how catalysts function to ultimately improve processes for the production of fuels and raw materials for the chemical industry from renewable biomass. Detailed computer simulations can play a key role in understanding how these catalysts function at a molecular level and guide the development of improved catalysts for industrially viable processes in biorefineries.

In a collaboration with Philippe Sautet (UCLA), Paul Fleurat-Lessard (Université de Bourgogne) and Carine Michel and Stephan Steinmann (ENS Lyon) we have developed a new force field for water over platinum that is more accurate than any existing model. We implemented the force field in the Amber software.

GAL17 can be combined with different water models and reparameterized for metals other than platinum. We are currently extending GAL17 to chemical compounds with alcohol groups at metal surfaces in aqueous solution. This will allow efficient sampling of the structure of water and biomass processing reactants/products at the metal/liquid interface and determination of solvation free energies along catalytic reaction pathways.

[1] S. Steinmann, R. Ferreira de Morais, A. W. Götz, P. Fleurat-Lessard, M. Iannuzzi, P. Sautet, C. Michel, J. Chem. Theory Comput. 14, 3238-3251 (2018). DOI: 10.1021/acs.jctc.7b01177

In recent years, a new generation of DFT methods emerged that uses functionals that not only depend on the elctron density and its gradients but also on Kohn-Sham orbitals. These DFT methods require the solution of the optimized effective potential (OEP) equations which are plagued by numerical stability problems if localized basis sets are employed. Together with Dr. A. Heßelmann I have developed a numerically stable approach which can be used for ground state calculations with orbital-dependent functionals [1,2]. Our method has been implemented and tested for exact-exchange (EXX) DFT which employs the exact exchange functional instead of approximations to it.

Time-dependent density functional theory (TD-DFT) is a very efficient and in many cases encouragingly accurate method to describe electronically excited states of molecules. A fundamental problem which has not yet been solved concerns excited states with charge-transfer (CT) character. Such excited states, which are of importance for a great deal of applications, cannot be described by GGA XC functionals. A very promising avenue to solve this problem is the extension of the EXX approach to excited states [3]. Although the theory has been worked out already in 1998, no practical implementation for molecules exists so far.

LEDO is an acronym for limited expansion of diatomic overlap densities and is a novel density fitting approach to speed up molecular DFT calculations. Conventional density fitting methods (also called RI methods) as generally employed in DFT expand the complete electron density into a fit basis which is distributed over the whole system (molecule). As a result one- to three-center two-electron repulsion integrals (ERIs) have to be evaluated and the computational expense is O(N³) with respect to the size N of the system under investigation as compared to O(N⁴) for conventional DFT.

In the LEDO-DFT formalism each individual diatomic overlap density is expanded into a fit basis which is restricted to the atoms of that overlap density. Thus, two-center matrix elements of the Coulomb- and Exchange-Correlation contribution to the Kohn-Sham matrix can be expressed in terms of the one-center elements using the LEDO expansion coefficients. Consequently only one- and two-center ERIs have to be evaluated and the computational expense of the Kohn-Sham matrix formation in a LEDO-DFT calculation is reduced to O(N²).