Scalable evaluation of polarization energy and associated forces in polarizable molecular dynamics: I. Toward massively parallel direct space computations

In this paper, we investigate various numerical strategies to compute the direct space polarization energy and associated forces in the context of the point dipole approximation (including damping) used in polarizable molecular dynamics. We present a careful mathematical analysis of the algorithms that have been implemented in popular production packages and applied to large test systems. We show that the classical Jacobi Over-Relaxation method (JOR) should not be used as its convergence requires a proper value of the relaxation parameter, whereas other strategies should be preferred. On a single node, Preconditioned Conjugate Gradient methods (PCG) and Jacobi algorithm coupled with the Direct Inversion in the Iterative Subspace (JI/DIIS) provide reliable stability/convergence and are roughly twice as fast as JOR. Moreover, both algorithms are suitable for massively parallel implementations. The lower requirements in terms of processes communications make JI/DIIS the method of choice for MPI and hybrid OpenMP/MPI paradigms for real life tests. Furthermore, using a predictor step as a guess along a molecular dynamics simulation provides another inexpensive, yet very effective, form of convergence acceleration. Overall, two to three orders of magnitude in time can be gained compared to the initial JOR single node approach to the final PGC or JI/DIIS parallel one combined with the predictors MD refinements. Such a speedup traces a new route for the high performance implementation of polarizable molecular dynamics and therefore extends the applicability of the technique as it will facilitate future multiscale QM/MM/continuum computations.