Is it possible to couple MD and DPD in lammps?

Thank you so much for your detailed explanation, Steve.
I still have 3 questions:

1. In my opinion, using fixes like "fix deposit" and "fix evaporate", I can easily add or remove atoms and rebuild the neighbor list at the same time, right?

Not for your purpose.

2. In my simulation, I would like to use different time steps for MD region and DPD region, is it possible?

No.

3. The purpose of using multiscale methods is to save computational time, suppose I use 128 processes in my simulation, and the MD region is half of the whole computational region, then there are 64 processes used for MD simulation; if MD region is a quater of the whole region, then there are 32 processes used for MD simulation... It seems the computational time is always the same no matter how many percentage MD region account for (the time consumed by MD region is dominant). Thus it's not effecient at all. Therefore I wonder if it is possible to use two serial Lammps programs, one for MD simulation, the other for DPD simulation. For example, after 128 processes are used for MD region for 1000 time steps, they are then used for DPD region for 100 time steps (suppose 10 MD time steps equals 1 DPD time step).

It is not that simple. You need to write a new integrator for that,
since you need a consistent way to couple the two regions. This has
been done before. I suggest to read up on QM/MM coupling, coarse grain
with all-atom coupling and integrators like r-RESPA before trying to
reinvent the wheel and repeating errors that others have done.

Axel