Residue-level coarse-grained (CG) molecular dynamics (MD) simulation is widely used to investigate slow biological processes that involve multiple proteins, nucleic acids, and their complexes. Biomolecules in a large simulation system are distributed non-uniformly, limiting computational efficiency with conventional methods. Here, we develop a hierarchical domain decomposition scheme with dynamic load balancing for heterogeneous biomolecular systems to keep computational efficiency even after drastic changes in particle distribution. The new schemes are applied to intrinsically disordered protein (IDP) droplet fusions. The droplet shape changes correlate with mixing IDP chains from two droplets. We also simulate formations of large IDP droplets, whose sizes are almost equivalent to those observed in microscopy. The methods have been implemented in CGDYN of the GENESIS software, which provides a new tool for investigating mesoscopic biological phenomena using the residue-level CG models.