Feature: add Hessian operator <\phi|\nabla_x\nabla_y|\phi>

[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 #6889

Unit Tests and/or Case Tests for my changes

================================================================================
TEST RESULTS

Total Tests: 15
Passed: 15 (100%)
Failed: 0 (0%)
Execution Time: < 0.01 seconds

Test Breakdown:

1. Zeros ✓ PASS
2. HessianSymmetryL5 ✓ PASS
3. HessianSymmetryL6 ✓ PASS
4. HessianFiniteDifferenceL5 ✓ PASS
5. HessianFiniteDifferenceL6 ✓ PASS
6. HessianL7NotImplemented ✓ PASS
7. HessianAllComponentsL2 ✓ PASS
8. HessianM0DifferentL ✓ PASS
9. HessianSpecialPointsL4 ✓ PASS
10. HessianTraceL3 ✓ PASS
11. HessianRotationalInvariance ✓ PASS
12. HessianL0Constant ✓ PASS
13. HessianL1Linear ✓ PASS
14. HessianNumericalStability ✓ PASS
15. HessianLargeCoordinates ✓ PASS

================================================================================
COVERAGE SUMMARY

Angular Momentum (l):
l=0: ✓ Tested (analytical - constant function)
l=1: ✓ Tested (analytical - linear functions)
l=2: ✓ Tested (all 6 Hessian components)
l=3: ✓ Tested (trace property, stability)
l=4: ✓ Tested (special points on axes)
l=5: ✓ Tested (finite difference validation)
l=6: ✓ Tested (finite difference validation)
l>6: ✓ Tested (boundary condition - not implemented)

Hessian Components (for l=2):
H_xx: ✓ Tested
H_xy: ✓ Tested
H_xz: ✓ Tested
H_yy: ✓ Tested
H_yz: ✓ Tested
H_zz: ✓ Tested

Test Points:
General coordinates: ✓ Tested
Coordinate axes (x,y,z): ✓ Tested
Small values (~1e-3): ✓ Tested
Large values (~100-300): ✓ Tested

================================================================================
VALIDATION METHODS

1. Analytical Verification (l=0, l=1)

   • Exact zero derivatives expected
   • Tolerance: 1e-10
   • Result: ✓ All zeros confirmed

2. Central Finite Difference (l=2, l=5, l=6)

   • Formula: H_ij ≈ [∂f/∂x_i(x_j+h) - ∂f/∂x_i(x_j-h)] / (2h)
   • Step size: h = 1e-5
   • Tolerance: 1e-3
   • Result: ✓ Within tolerance

3. Property Verification

   • Trace (Laplacian) properties
   • Rotational invariance
   • Result: ✓ Properties satisfied

4. Stability Testing

   • No NaN or Inf values
   • Consistent across coordinate ranges
   • Result: ✓ Numerically stable

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