Biomolecular Simulations in a Continuum Ionic Solvent with Polarizable Force Fields, Using Python and GPUs

Abstract

In molecular simulations, implicit-solvent models use continuum electrostatic theory to average the solvent degrees of freedom, drastically decreasing the amount of computations compared with explicit molecular dynamics. Mathematically, it can be formulated as a coupled system of the Poisson and Poisson-Boltzmann equations, interfaced by the biomolecule's solvent-excluded surface. Usually, these models use point-charge-based (classical) force fields to parameterize the biomolecule of interest, yet, more elaborate force fields that consider higher order multipoles and polarizability have emerged lately. Polarizable force fields generate more realistic charge distributions, but involve larger computations than a classical force field. One popular polarizable force field is AMOEBA, which uses point multipoles up to the quadrupole order to describe the charge distribution, and allows for the dipole component to polarize. In this work, we present an extension of the Poisson-Boltzmann solver PyGBe, which couples with the AMOEBA force field by solving iteratively for a self-consistent reaction potential. PyGBe is a boundary element method code, accelerated with a treecode algorithm and GPUs, that works from a Python interface. The boundary integral formulation treats the point multipoles analytically, which is an advantage over volumetric-based methods where very fine meshes are required to effectively place the point multipoles onto the numerical mesh. We will show results verifying our implementation in comparison with analytical solutions, available for spherical molecules, and computing the solvation energy of larger test proteins.