An O(n) Framework for Internal Coordinate Molecular Dynamics Applicable to Molecules with Arbitrary Constraints and Geometries

Abstract

Molecular Dynamics (MD) is a numerical simulation technique which is used to obtain the time evolution trajectory of a system of interacting particles. Consideration of the molecular movement in the space of generalized internal coordinate rather than in Cartesian coordinate is an efficient way to deal with the constraints such as fixed angles and lengths. This type of molecular dynamics is called Internal Coordinate Molecular Dynamics (ICMD). To integrate the equation we need to invert the dense mass matrix. Since the direct calculation of mass matrix inversion is cubically scaled with the number of atoms in a molecule, a more efficient method is required for macro-molecules simulation. Fixman's work in 1974 and the follow-up studies have developed a method that can factorize the inverse of mass matrix into an arithmetic combination of three sparse matrices—one of them is positive definite and need to be further factorized by using the Cholesky decomposition or similar methods. When the molecule subjected to study is of serial chain structure, this method can achieve $\mathcal{O}(n)$ computational scaling, where $n$ is the number of atoms within a specific molecule. However, for molecules with long branches and loops, the nonzero structure of this positive definite matrix makes its decomposition in scaling of $\mathcal{O}(n^3)$. We have presented a new method which can guarantee for no fill-in in doing the Cholesky decomposition. As a result, the inverting of mass matrix will remain the $\mathcal{O}(n)$ scaling, no matter the molecule structure has long branches or not. Based on the above $\mathcal{O}(n)$ method in inverting mass matrix, an $\mathcal{O}(n)$ framework for internal coordinated molecular dynamics has been built. It has the following properties: has $\mathcal{O}(n)$ time complexity; constraints such as constant bond lengths and angles can be arbitrarily applied\added{ on tree structures}; suitable to complex geometries such as branches, \added{flexible }loops, dummy atoms/sites; and is singularity free in inverting mass matrix. When compared to traditional MD methods, the new method is especially powerful in macromolecule simulations. The ICMD utilizes the bond angle $\theta_i$, torsional angle $\phi_i$, and bond length $b_i$ as the basic types of internal coordinate parameters. However, there is a serious problem of using $\theta_i$ and $\phi_i$ as the internal coordinates. When $\theta_i$ is very close to $0$ or $\pi$, the mass matrix becomes singular and its inverse does not exist, so that the numerical integration of the equations of motion is unstable in this situation. In this work, we will introduce a convention switch method which can prevent the singularity by adopting two alternative rotation conventions. According to our test, this method is suitable to any size of molecule. In addition, it only requires a very small amount of additional computational cost. The conservation of total energy in microcanonical (NVE) ensemble simulation is extremely important because the NVE ensemble simulation is a statistical sampling over the constant energy surface in the hyperdimensional phase space. A symplectic numerical scheme can make the error of the total energy be bounded even for exponentially large simulation time. On the other hand, non-symplect time integrators such as 4th order Runge-Kutta (RK4) would give monotonic increasing or decreasing total energy. In internal coordinate molecular dynamics, the Hamiltonian is a strongly coupled equation which depends both on the momentum and position vectors, and very few explicit symplectic numerical methods existing for inseparable Hamiltonian. In this work, the influence of symplectic and non-symplectic time integrators on the energy conservation of ICMD have been studied.