General force fields such as General Amber Force Field (GAFF) have been designed for broad applicability and are widely used in protein-ligand binding simulations in structure-based drug discovery. However, the force field parameters are not always transferable across ligand molecules, and custom reparameterization is sometimes necessary for accurate binding free energy simulations. This is especially true for torsion parameters, which are highly dependent on stereoelectronic and steric effects. Here, we report a novel, flexible, and user-friendly computational tool called the Automated Force Field Developer and Optimizer (AFFDO) platform that allows generating accurate, tailored GAFF2 torsion parameters for drug-like molecules. For a given ligand, AFFDO selects the most important torsions, carries out GPU-accelerated density functional theory calculations to collect reference data and fits torsion terms using a fast gradient-based optimizer that leverages automated differentiation. We benchmark AFFDO by parametrizing a series of drug-like molecules and carrying out protein-ligand relative binding free energy (RBFE) simulations. The results show that AFFDO can significantly improve GAFF2 torsion parameters against QM reference data, which in some cases translates into better agreement with experimental RBFE values within a reasonable computational time.

The acoustic inverse scattering problem is of critical importance in a number of fields, including medical imaging, sonar, and non-destructive evaluation. The problem of interest can vary from the detection of the shape to the properties of an obstacle. The challenge is that this problem is severely ill-posed and highly nonlinear. Significant efforts have been expended over the years to develop solutions to this problem. However, existing fast data-driven methods primarily focus on the two-dimensional scattering case. This paper explores the potential of using machine learning to accelerate the solution to the three-dimensional (3D) version of the problem. To this end, we develop inverse scattering shape reconstruction network (ISSRNet), a deep learning framework for 3D shape reconstruction using phaseless far-field data. The framework is implemented by (a) using a compact probabilistic shape latent space learned by a 3D variational auto-encoder, and (b) a convolutional neural network trained to extract far-field features due to multiple incident waves and map the acoustic scattering information to this shape representation. We demonstrate ISSRNet's 3D shape reconstruction capabilities on random rock-like particles, and airplane objects from the popular ShapeNet data set. We also evaluate the framework's performance when trained on lower-resolution scattering data and when receiver locations include uncertainty. Our experiments show that the proposed framework is able to capture both global and local details, differentiate between different types of shapes and performs several orders of magnitude faster than its numerical iterative counterparts.

This paper reformulates complementarity-based time-stepping for frictionless nonsmooth contact between smooth rigid bodies as a recursively generated linear complementarity problem (ReLCP), involving a sequence of LCPs of increasing dimension. Starting from a classical single-constraint shared-normal signed-distance (SNSD) LCP, the method adds unilateral constraints only when the discrete-time update predicted by the current contact set would violate nonpenetration of the underlying smooth surfaces. The resulting procedure acts directly on smooth geometry, enforces nonpenetration to a prescribed tolerance, and avoids the oversampling inherent to proxy-surface contact models such as tessellations or multi-sphere decompositions, for which improved geometric fidelity can drive rapid growth in constraint count and cost. For strictly convex bodies, we prove that an initially overlap free configuration with sufficiently small timestep sizes, imply finite termination of the adaptive augmentation, and yield a unique discrete-time velocity update. In the small timestep limit and for any fixed overlap-free discrete state with a fixed geometric overlap tolerance, we prove that the recursion terminates after the initial solve, reducing the method to the classical single-constraint SNSD LCP and retaining the usual consistency of complementarity time-stepping with the underlying differential variational inequality. Numerical tests on colliding ellipsoids, compacting ellipsoid suspensions, growing bacterial colonies, and taut chainmail networks demonstrate stable large-timestep behavior, bounded interpenetration without discretization-induced surface roughness, and substantial reductions in both active constraint counts and runtime relative to representative discrete-surface complementarity formulations.

MFDn (Many-body Fermion Dynamics for nucleons) is a cutting-edge Configuration Interaction (CI) code designed to tackle nuclear quantum many-body problems. The core computational task in MFDn involves solving a large sparse eigenvalue problem using iterative methods, where the most costly computational step is the sparse matrix-vector multiplication (SpMV). Recently, an MPI/OpenACC based SpMV was developed to enable MFDn to run efficiently on NVIDIA GPUs. In this work, we explore various strategies to further enhance the performance of MFDn by exploiting architecture features of GPUs and making use of alternative programming models such as CUDA kernels as well as communication protocols such as the asynchronous point-to-point (P2P) message passing and NVIDIA Collective Communication Library (NCCL). We demonstrate the performance gain achieved by using these techniques. In particular, we show that, on problem sizes with up to 1.3 trillion nonzeros, we can obtain up to $2.0\times$ improvement in the overall solver time by switching from OpenACC to a hand-optimized CUDA implementation across 1540 GPUs. Additionally, when combined with asynchronous P2P message passing and NCCL, we observe performance boosts of $2.9\times$ (P2P) and $4.9\times$ (NCCL) compared to the baseline optimized version using the MPI/OpenACC programming model. Our study has uncovered a few limitations of the existing directive based programming models such as OpenACC and highlights the challenges of using CUDA-aware MPI for certain collective communications. While the primary focus of this work is on optimizing the performance of MFDn, the insights and solutions we provide are likely to be relevant to a broad range of applications.

We introduce OpenRAND, a C++17 library aimed at facilitating reproducible scientific research by generating statistically robust yet replicable random numbers in as little as two lines of code, overcoming some of the unnecessary complexities of existing RNG libraries. OpenRAND accommodates single and multi-threaded applications on CPUs and GPUs and offers a simplified, user-friendly API that complies with the C++ standard’s random number engine interface. It is lightweight; provided as a portable, header-only library. It is statistically robust: a suite of built-in tests ensures no pattern exists within single or multiple streams. Despite its simplicity and portability, it remains performant—matching and sometimes outperforming native libraries. Our tests, including a Brownian walk simulation, affirm its reproducibility and ease-of-use while highlight its computational efficiency, outperforming CUDA’s cuRAND by up to 1.8 times.

The inverse scattering problem is of critical importance in a number of fields, including medical imaging, sonar, sensing, non-destructive evaluation, and several others. The problem of interest can vary from detecting the shape to the constitutive properties of the obstacle. The challenge in both is that this problem is ill-posed, more so when there is limited information. That said, significant effort has been expended over the years in developing solutions to this problem. Here, we use a different approach, one that is founded on data. Specifically, we develop a deep learning framework for shape reconstruction using limited information with single incident wave, single frequency, and phase-less far-field data. This is done by (a) using a compact probabilistic shape latent space, learned by a 3D variational auto-encoder, and (b) a convolutional neural network trained to map the acoustic scattering information to this shape representation. The proposed framework is evaluated on a synthetic 3D particle dataset, as well as ShapeNet, a popular 3D shape recognition dataset. As demonstrated via a number of results, the proposed method is able to produce accurate reconstructions for large batches of complex scatterer shapes (such as airplanes and automobiles), despite the significant variation present within the data.

We report the development and testing of new integrated cyberinfrastructure for performing free energy simulations with generalized hybrid quantum mechanical/molecular mechanical (QM/MM) and machine learning potentials (MLPs) in Amber. The Sander molecular dynamics program has been extended to leverage fast, density-functional tight-binding models implemented in the DFTB+ and xTB packages, and an interface to the DeePMD-kit software enables the use of MLPs. The software is integrated through application program interfaces that circumvent the need to perform "system calls" and enable the incorporation of long-range Ewald electrostatics into the external software's self-consistent field procedure. The infrastructure provides access to QM/MM models that may serve as the foundation for QM/MM-Delta MLP potentials, which supplement the semiempirical QM/MM model with a MLP correction trained to reproduce ab initio QM/MM energies and forces. Efficient optimization of minimum free energy pathways is enabled through a new surface-accelerated finite-temperature string method implemented in the FE-ToolKit package. Furthermore, we interfaced Sander with the i-PI software by implementing the socket communication protocol used in the i-PI client-server model. The new interface with i-PI allows for the treatment of nuclear quantum effects with semiempirical QM/MM-Delta MLP models. The modular interoperable software is demonstrated on proton transfer reactions in guanine-thymine mispairs in a B-form deoxyribonucleic acid helix. The current work represents a considerable advance in the development of modular software for performing free energy simulations of chemical reactions that are important in a wide range of applications.

AbstractThis chapter describes recent advances in the use of machine learning techniques in reactive atomistic simulations. In particular, it provides an overview of techniques used in training force fields with closed form potentials, developing machine-learning-based potentials, use of machine learning in accelerating the simulation process, and analytics techniques for drawing insights from simulation results. The chapter covers basic machine learning techniques, training procedures and loss functions, issues of off-line and in-lined training, and associated numerical and algorithmic issues. The chapter highlights key outstanding challenges, promising approaches, and potential future developments. While the chapter relies on reactive atomistic simulations to motivate models and methods, these are more generally applicable to other modeling paradigms for reactive flows.

We examine and compare several iterative methods for solving large-scale eigenvalue problems arising from nuclear structure calculations. In particular, we discuss the possibility of using block Lanczos method, a Chebyshev filtering based subspace iterations and the residual minimization method accelerated by direct inversion of iterative subspace (RMM-DIIS) and describe how these algorithms compare with the standard Lanczos algorithm and the locally optimal block preconditioned conjugate gradient (LOBPCG) algorithm. Although the RMM-DIIS method does not exhibit rapid convergence when the initial approximations to the desired eigenvectors are not sufficiently accurate, it can be effectively combined with either the block Lanczos or the LOBPCG method to yield a hybrid eigensolver that has several desirable properties. We will describe a few practical issues that need to be addressed to make the hybrid solver efficient and robust.

We have ported and optimized the graphics processing unit (GPU)-accelerated QUICK and AMBER-based ab initio quantum mechanics/molecular mechanics (QM/MM) implementation on AMD GPUs. This encompasses the entire Fock matrix build and force calculation in QUICK including one-electron integrals, two-electron repulsion integrals, exchange-correlation quadrature, and linear algebra operations. General performance improvements to the QUICK GPU code are also presented. Benchmarks carried out on NVIDIA V100 and AMD MI100 cards display similar performance on both hardware for standalone HF/DFT calculations with QUICK and QM/MM molecular dynamics simulations with QUICK/AMBER. Furthermore, with respect to the QUICK/AMBER release version 21, significant speedups are observed for QM/MM molecular dynamics simulations. This significantly increases the range of scientific problems that can be addressed with open-source QM/MM software on state-of-the-art computer hardware.

AmberTools is a free and open-source collection of programs used to set up, run, and analyze molecular simulations. The newer features contained within AmberTools23 are briefly described in this Application note.

The reactive force field (ReaxFF) interatomic potential is a powerful tool for simulating the behavior of molecules in a wide range of chemical and physical systems at the atomic level. Unlike traditional classical force fields, ReaxFF employs dynamic bonding and polarizability to enable the study of reactive systems. Over the past couple decades, highly optimized parallel implementations have been developed for ReaxFF to efficiently utilize modern hardware such as multi-core processors and graphics processing units (GPUs). However, the complexity of the ReaxFF potential poses challenges in terms of portability to new architectures (AMD and Intel GPUs, RISC-V processors, etc.), and limits the ability of computational scientists to tailor its functional form to their target systems. In this regard, the convergence of cyber-infrastructure for high performance computing (HPC) and machine learning (ML) presents new opportunities for customization, programmer productivity and performance portability. In this paper, we explore the benefits and limitations of JAX, a modern ML library in Python representing a prime example of the convergence of HPC and ML software, for implementing ReaxFF. We demonstrate that by leveraging auto-differentiation, just-in-time compilation, and vectorization capabilities of JAX, one can attain a portable, performant, and easy to maintain ReaxFF software. Beyond enabling MD simulations, end-to-end differentiability of trajectories produced by ReaxFF implemented with JAX makes it possible to perform related tasks such as force field parameter optimization and meta-analysis without requiring any significant software developments. We also discuss scalability limitations using the current version of JAX for ReaxFF simulations.

Many integral equations used to analyze scattering, such as the standard combined field integral equation (CFIE), are not well-conditioned for a wide range of frequencies and multi-scale geometries. There has been significant effort to alleviate this problem. A more recent one is using a set of decoupled potential integral equations (DPIE). These equations have been shown to be robust at low frequencies and immune to topology breakdown. But they mimic the ill-conditioning behavior of CFIE at high frequencies. This paper addresses this deficiency through new Calderón-type identities derived from the Vector Potential Integral Equation (VPIE). We construct novel analytic preconditioners for the vector potential integral equation (VPIE) and scalar potential integral equation (SPIE) constrained to perfect electric conductors (PEC). These new formulations are wide-band well-conditioned and converge rapidly for multi-scale geometries. This is demonstrated though a number of examples that use analytic and piecewise basis sets.

The growing interest in the effects of external electric fields on reactive processes requires predictive methods that can reach longer length and time scales than quantum mechanical simulations. Recently, many studies have included electric fields in ReaxFF, a widely used reactive molecular dynamics method. In the case of modeling an external electric field, the charge distribution method used in ReaxFF is critical. The most common charge distribution method used in previous studies of electric fields is the charge equilibration (QEq) method, which assumes that the system is a contiguous conductor and that charge transfer can occur across any distance. In contrast, many systems of interest are insulators or semiconductors, and long-distance charge transfer should not occur in response to a small difference in potential. This study focuses on the limitations of the QEq method in the context of water in an external electric field. We demonstrate that QEq can predict unphysical charge distributions and exhibits properties that do not converge as a function of system size. Furthermore, we show that electric fields within the recently developed atom-condensed Kohn-Sham density functional theory (DFT) approximated to the second-order (ACKS2) approach address the major limitations of electric fields in QEq. With ACKS2, we observe more physical charge distributions and properties that converge as a function of system size. We do not suggest that ACKS2 is perfect in all circumstances but rather show specific cases where it addresses the major shortcomings of QEq in the context of an external electric field.

Evaluation of pair potentials is critical in a number of areas of physics. The classical N-body problem has its root in evaluating the Laplace potential, and has spawned tree-algorithms, the fast multipole method (FMM), as well as kernel independent approaches. Over the years, FMM for Laplace potential has had a profound impact on a number of disciplines as it has been possible to develop highly scalable parallel versions of these algorithms. This is in stark contrast to parallel algorithms for oscillatory potentials such as the Helmholtz potential. The principal bottlenecks to scalable parallelism are the computation and communication costs of operations necessary to traverse up, across, and down the tree. In this article, we analyze asymptotic costs for both computation and communication in a parallel implementation, and describe techniques to overcome bottlenecks and achieve high performance evaluation of the Helmholtz potential for different distributions of particles. We demonstrate that the resulting implementation has a load balancing effect that significantly reduces the time-to-solution and enhances the scale of problems that can be treated using full wave physics.

The reactive force field (ReaxFF) model bridges the gap between traditional classical models and quantum mechanical (QM) models by incorporating dynamic bonding and polarizability. To achieve realistic simulations using ReaxFF, model parameters must be optimized against high fidelity training data which typically come from QM calculations. Existing parameter optimization methods for ReaxFF consist of black box techniques using genetic algorithms or Monte Carlo methods. Due to the stochastic behavior of these methods, the optimization process oftentimes requires millions of error evaluations for complex parameter fitting tasks, thereby significantly hampering the rapid development of high quality parameter sets. Rapid optimization of the parameters is essential for developing and refining Reax force fields because producing a force field which exhibits empirical accuracy in terms of dynamics typically requires multiple refinements to the training data as well as to the parameters under optimization. In this work, we present JAX-ReaxFF, a novel software tool that leverages modern machine learning infrastructure to enable fast optimization of ReaxFF parameters. By calculating gradients of the loss function using the JAX library, JAX-ReaxFF utilizes highly effective local optimization methods that are initiated from multiple guesses in the high dimensional optimization space to obtain high quality results. Leveraging the architectural portability of the JAX framework, JAX-ReaxFF can execute efficiently on multicore CPUs, graphics processing units (GPUs), or even tensor processing units (TPUs). As a result of using the gradient information and modern hardware accelerators, we are able to decrease ReaxFF parameter optimization time from days to mere minutes. Furthermore, the JAX-ReaxFF framework can also serve as a sandbox environment for domain scientists to explore customizing the ReaxFF functional form for more accurate modeling.

A novel locally polarizable multisite model based on the original cation dummy atom (CDA) model is described for molecular dynamics simulations of ions in condensed phases. Polarization effects are introduced by the electronegativity equalization model (EEM) method where charges on the metal ion and its dummy atoms can fluctuate to respond to the environment. This model includes explicit polarization and ion-induced interactions and can be coupled with nonpolarizable or polarizable water models, making it more transferable to simpler force fields. This approach allows us to enhance the original fixed charge CDA model where the charge distribution cannot adapt to the local solvent structure. To illustrate the new CDApol model, we examined properties of the Zn2+, Al3+, and Zr4+ ions in aqueous solution. The polarizable model and Lennard-Jones parameters were refined for octahedrally coordinated Zn2+, Al3+, and Zr4+ CDAs to reproduce thermodynamic and geometrical properties. Using this locally polarizable model, we were able to obtain the experimental hydration free energy, ion-oxygen distance, and coordination number coupled with the standard 12-6 Lennard-Jones model. This model can be used in myriad additional applications where local polarization and charge transfer is important.

Since the classical molecular dynamics simulator LAMMPS was released as an open source code in 2004, it has become a widely-used tool for particle-based modeling of materials at length scales ranging from atomic to mesoscale to continuum. Reasons for its popularity are that it provides a wide variety of particle interaction models for different materials, that it runs on any platform from a single CPU core to the largest supercomputers with accelerators, and that it gives users control over simulation details, either via the input script or by adding code for new interatomic potentials, constraints, diagnostics, or other features needed for their models. As a result, hundreds of people have contributed new capabilities to LAMMPS and it has grown from fifty thousand lines of code in 2004 to a million lines today. In this paper several of the fundamental algorithms used in LAMMPS are described along with the design strategies which have made it flexible for both users and developers. We also highlight some capabilities recently added to the code which were enabled by this flexibility, including dynamic load balancing, on-the-fly visualization, magnetic spin dynamics models, and quantum-accuracy machine learning interatomic potentials. Program Summary Program Title: Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) CPC Library link to program files: https://doi .org /10 .17632 /cxbxs9btsv.1 Developer's repository link: https://github .com /lammps /lammps Licensing provisions: GPLv2 Programming language: C++, Python, C, Fortran Supplementary material: https://www.lammps .org Nature of problem: Many science applications in physics, chemistry, materials science, and related fields require parallel, scalable, and efficient generation of long, stable classical particle dynamics trajectories. Within this common problem definition, there lies a great diversity of use cases, distinguished by different particle interaction models, external constraints, as well as timescales and lengthscales ranging from atomic to mesoscale to macroscopic. Solution method: The LAMMPS code uses parallel spatial decomposition, distributed neighbor lists, and parallel FFTs for long-range Coulombic interactions [1]. The time integration algorithm is based on the Stormer-Verlet symplectic integrator [2], which provides better stability than higher-order non-symplectic methods. In addition, LAMMPS supports a wide range of interatomic potentials, constraints, diagnostics, software interfaces, and pre- and post-processing features. Additional comments including restrictions and unusual features: This paper serves as the definitive reference for the LAMMPS code. References [1] S. Plimpton, Fast parallel algorithms for short-range molecular dynamics. J. Comp. Phys. 117 (1995) 1-19. [2] L. Verlet, Computer experiments on classical fluids: I. Thermodynamical properties of Lennard-Jones molecules, Phys. Rev. 159 (1967) 98-103. (c) 2021 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

Programming applications on heterogeneous systems with hardware accelerators is challenging due to the disjoint address spaces between the host (CPU) and the device (GPU). The limited device memory further exacerbates the challenges as most data-intensive applications will not fit in the limited device memory. CUDA Unified Memory (UM) was introduced to mitigate such challenges. UM improves GPU programmability by supporting oversubscription, on-demand paging, and migration. However, when the working set of an application exceeds the device memory capacity, the resulting data movement can cause significant performance losses. We propose a tiling-based task-parallel framework, named DeepSparseGPU, to accelerate sparse eigensolvers on GPUs by minimizing data movement between the host and device. To this end, we tile all operations in a sparse solver and express the entire computation as a directed acyclic graph (DAG). We design and develop a memory manager (MM) to execute larger inputs that do not fit into GPU memory. MM keeps track of the data on CPU and GPU, and automatically moves data between them as needed. We use OpenMP target offload in our implementation to achieve portability beyond NVIDIA hardware. Performance evaluations show that DeepSparseGPU transfers 1.39x-2.18x less host to device (H2D) and device to host (D2H) data, while executing up to 2.93x faster than the UM-based baseline version.

Recently, integral equation formulations that use potentials as opposed to fields as unknown quantities have been developed for scattering from dielectric objects. It has been shown that these formulations can be construed so that they are well-conditioned across a broad frequency spectrum, a result that has been theoretically proven for spherical systems. Unfortunately, to date, this formulation has not been implemented on practical discretizations of objects. This is the goal of this article. Specifically, we present a well-conditioned and well-tested decoupled potential integral equation (DPIE) formulation and all the necessary implementation details for electromagnetic scattering from homogeneous, dielectric, arbitrarily shaped objects. The resulting decoupled systems do not suffer from low-frequency breakdown. Results that demonstrate these properties are presented for a number of different dielectric targets. Furthermore, in order to fully validate each of the two integral equations for the potentials, we develop analytical solutions for spherical systems.