K-means clustering, as a classic unsupervised machine learning algorithm, is the key step to select the interpolation sampling points in interpolative separable density fitting (ISDF) decomposition. Real-valued K-means clustering for accelerating the ISDF decomposition has been demonstrated for large-scale hybrid functional enabled \textit{ab initio} molecular dynamics (hybrid AIMD) simulations within plane-wave basis sets where the Kohn-Sham orbitals are real-valued. However, it is unclear whether such K-means clustering works for complex-valued Kohn-Sham orbitals. Here, we apply the K-means clustering into hybrid AIMD simulations for complex-valued Kohn-Sham orbitals and use an improved weight function defined as the sum of the square modulus of complex-valued Kohn-Sham orbitals in K-means clustering. Numerical results demonstrate that this improved weight function in K-means clustering algorithm yields smoother and more delocalized interpolation sampling points, resulting in smoother energy potential, smaller energy drift and longer time steps for hybrid AIMD simulations compared to the previous weight function used in the real-valued K-means algorithm. In particular, we find that this improved algorithm can obtain more accurate oxygen-oxygen radial distribution functions in liquid water molecules and more accurate power spectrum in crystal silicon dioxide compared to the previous K-means algorithm. Finally, we describe a massively parallel implementation of this ISDF decomposition to accelerate large-scale complex-valued hybrid AIMD simulations containing thousands of atoms (2,744 atoms), which can scale up to 5,504 CPU cores on modern supercomputers.