Fix: DeePKS MD/Relax/Cell-relax calculation convergence

[dyzheng]

Reminder

• Have you linked an issue with this pull request?
• Have you added adequate unit tests and/or case tests for your pull request?
• Have you noticed possible changes of behavior below or in the linked issue?
• Have you explained the changes of codes in core modules of ESolver, HSolver, ElecState, Hamilt, Operator or Psi? (ignore if not applicable)

Linked Issue

Fix #4525

Unit Tests and/or Case Tests for my changes

• A unit test is added for each new feature or bug fix.

What's changed?

• Example: My changes might affect the performance of the application under certain conditions, and I have tested the impact on various scenarios...

Any changes of core modules? (ignore if not applicable)

• Example: I have added a new virtual function in the esolver base class in order to ...

[xuan112358]

@dyzheng I note that this pr is mainly about reusing pdm, which is beneficial to reduce the number of MD/Relax/Cell-relax steps.
However, pr #3256 seems to be mainly related to the way dm is referenced and the use of BLAS library functions for computation (Am I right here?). And I have no idea why it affects the number of steps. I am worring that #3256 introduces some problems that are not addressed here.

[mohanchen]

We will keep improving DeePKS in the near future.