Performance Tuning of Polarizable Gaussian Multipole Model in Molecular Dynamics Simulations

Abstract

Molecular dynamics (MD) simulations are essential for understanding molecular phenomena at the atomic level, with their accuracy largely dependent on both the employed force field and sampling. Polarizable force fields, which incorporate atomic polarization effects, represent a significant advancement in simulation technology. The polarizable Gaussian multipole (pGM) model has been noted for its accurate reproduction of ab initio electrostatic interactions. In this study, we document our effort to enhance the computational efficiency and scalability of the pGM simulations within the AMBER framework using MPI (message passing interface). Performance evaluations reveal that our MPI-based pGM model significantly reduces runtime and scales effectively while maintaining computational accuracy. Additionally, we investigated the stability and reliability of the MPI implementation under the NVE simulation ensemble. Optimal Ewald and induction parameters for the pGM model are also explored, and its statistical properties are assessed under various simulation ensembles. Our findings demonstrate that the MPI-implementation maintains enhanced stability and robustness during extended simulation times. We further evaluated the model performance under both NVT (constant number, volume, and temperature) and NPT (constant number, pressure, and temperature) ensembles and assessed the effects of varying timesteps and convergence tolerance on induced dipole calculations. The lessons learned from these exercises are expected to help the users to make informed decisions on simulation setup. The improved performance under these ensembles enables the study of larger molecular systems, thereby expanding the applicability of the pGM model in detailed MD simulations.

Figure 1.

Workflowchart for the pGM model and the MPI communications applied within the pGM model. The yellow blocks indicate the operation completed within a single thread, while the green blocks represent operations performed locally in each thread. The induced dipole evaluation involves proposing an initial guess, followed by self-iteration, which includes calculating the dipole–dipole interactions and updating the induced dipole.

Figure 2.

In the pGM model, the field calculation can be divided into three parts: the direct part, the reciprocal part and the self-correction part. The direct and self-correction parts use atom-based decomposition in MPI, while the reciprocal part is based on the spatial decomposition where in SANDER, we used slab decomposition.

Figure 3.

Scalability of the pGM MPI model for the water model across different system sizes with SANDER. In the figure, all simulations are conducted under the NVE ensemble, with each simulated for 10 ps. The efficiency is calculated as the time taken divided by the time taken for the serial code simulation.

Figure 4.

Efficiency of the pGM water model utilizing SANDER compared to the point-charge nonpolarizable TIP3P water model.

Figure 5.

pGM simulation of 512 water model of 1 ns under NVE ensemble, the MPI and serial results are compared. Identical parameters were used in both simulations. Although there is small difference, both codes achieve energy conservation and thus proves the pGM MPI is stable and robust. The energy drift values for the are 1.13 × 10⁻⁴ kcal/mol ps mol in the MPI simulation and 9.33 × 10⁻⁵ kcal/mol ps in the serial simulation.

Figure 6.

Energy drift with grid densities in pGM simulations of 512-water box of 1 ns under NVE ensemble. The label shows the density of the grid along one direction. The energy drift values are −4.65 × 10⁻³, −5.05 × 10⁻⁴, 8.00 × 10⁻⁶, −9.00 × 10⁻⁶ kcal/mol ps for grid densities of 0.74, 1.11, 1.48, and 1.85 (per Å), respectively.

Figure 7.

Energy drift with different B-spline orders in pGM simulations of 512 water model of 1 ns under NVE ensemble. The energy drift values are −6.03 × 10⁻³, 3.39 × 10⁻⁴, −9.61 × 10⁻⁵,6.93 × 10⁻⁴, 4.23 × 10⁻⁵ kcal/mol ps for B-spline orders of 4, 5, 6, 7, and 8, respectively.

Figure 8.

Histograms of temperature distributions of pGM and TIP3P models at 300 K with the Langevin thermostat under NVT. The mean and standard deviation are 300.11 and 7.66 for the pGM model, and 299.49 and 7.66 for the TIP3P model, respectively.

Figure 9.

Histograms of relative potential energy distributions of pGM and TIP3P models at 300 K with the Langevin thermostat under NVT. The potential energy is relative to the mean value so that both distributions are centered at 0 kcal/mol. The standard deviations are 38.96 for the pGM water and 32.34 for the TIP3P water, respectively.

Figure 10.

Comparison of the velocity correlation functions of pGM and TIP3P models.

Figure 11.

Comparison of the density distributions of pGM and TIP3P models at 300 K with Langevin thermostat and Monte Carlo barostat under NPT. Mean values are 1.006 g/cm³ and 0.977 g/cm³ for pGM and TIP3P waters, respectively. The variances are 0.012 for both water models.