Analysis of the application of a high order symplectic method in Shardlow's method for dissipative particle dynamics

TL;DR

This study addresses the efficiency and reliability of the M-Shardlow method for DPD by identifying improper parameter usage in the Yoshida-based fourth-order symplectic update and evaluating its practical performance against Shardlow and ABOBA. By both theoretical BCH analysis and numerical experiments, it shows that even with corrected parameters, the M-Shardlow-2 scheme incurs significant CPU overhead and does not outperform existing methods in preserving system structure or achieving favorable observables. The results demonstrate at least third-order convergence to the invariant measure, but the cost-benefit balance favors established methods like ABOBA, highlighting the challenges of high-order symplectic schemes in dissipative particle dynamics. The findings inform the design of efficient DPD integrators and caution against unrealized gains from high-order Hamiltonian solvers without commensurate reductions in computational effort.

Abstract

This study investigates the efficiency and reliability of the modified Shardlow's (M-Shardlow) method for dissipative particle dynamics (DPD). We show that the M-Shardlow method in which for its construction, the second order velocity Verlet method in the Shardlows method to integrate the Hamiltonian part has been replaced by a symplectic fourth order method, improperly uses some parameters. %In other words, in this paper, it is shown that the initial M-Shardlow method employed some parameters improperly in the fourth order symplectic method that was used for the M-Shardlow method. By numerical experiments and computing, some important configurational quantities such as configurational temperature and radial distribution function (RDF), the M-Shardlow's method is compared with the Shardlow and ABOBA methods. These results indicate that the new method obtained in this way, even with the proper parameters is too costly in the sense of the CPU-time that is required per each step which makes it an inefficient DPD integrator. Besides, by a comparison of the radial distribution function of this method with Shardlow and ABOBA for large time increments, we can observe no considerable improvement in preserving the structure of the system by this new DPD solver.

Figures (5)

• Figure 1: The relative error in the computed of the kinetic and configurational temperature for Shardlow, ABOBA, M-Shardlow1, M-Shardlow2 and M-ABOBA methods with $(a,\gamma)=(18.75,4.5)$. The system was simulated for 1000 reduced time units but only the last $80\%$ of the data was collected to compute the static quantity to be sure that the system was well equilibrated. The black and blue solid lines respectively represent the first order and the third order and the black dashed line shows the second order convergence to the invariant measure.
• Figure 2: Double logarithmic of the relative error in the computed configurational temperature in Shardlow, ABOBA, M-Shardlow2, and M-ABOBA methods for $\gamma=4.5,40.5,200,450$ ($a=25$). The format of the plots is the same as in figure \ref{['fig1']}, except that the stepsize starts at $\Delta t=0.01$.
• Figure 3: Comparison of the radial distribution function (RDF), $g(r)$, for different methods with a friction coefficient of $\gamma=4.5$ and constant repulsion $a=25$. In all the methods the stepsize start at $\Delta t=0.09$ and RDFs have been plotted until large artifacts observed. The solid black line is the refrence solution obtained by using the Shardlow's method with a very small stepsize of $\Delta t=0.001$. The dashed colour lines correspond to different stepsize as indicated.
• Figure 4: Comparison of the radial distribution function (RDF), $g(r)$, for different methods with a large coefficient friction of $\gamma=40.5$ and constant repulsion $a=25$. In all the methods the stepsize start at $\Delta t=0.09$ and RDFs have been plotted until large artifacts observed. The format of the plot is the same as in Figure 3.
• Figure 5: Comparison of the radial distribution function (RDF), $g(r)$, for different methods with a extremely large coefficient friction of $\gamma=450$ and $a=25$. In all the methods the stepsize start at $\Delta t=0.09$ and RDFs have been plotted until large artifacts observed. The format of the plots is the same as in Figure 3.