Fix: calculate energy density matrix using 0.5(S^-1HD + DHS^-1)

[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

This PR should be tested and check if the current code is wrong.
Any one can repeat the error with multi-k calculation in rt-TDDFT MD simulation.

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 ...

[dyzheng]

Force of rt-TDDFT is still different with AIMD, both before and after this PR.

[mohanchen]

We will close this PR, but new PRs related to this topic are welcome