Refactor RDMFT_LCAO and RDMFT classes: Clean up constructors, enhance initialization, and update energy calculation methods

Author: Kai Luo

source/source_esolver/directmin_problems/rdmft_lcao.cpp

@@ -40,20 +40,67 @@ RDMFT_LCAO<TK, TR>::RDMFT_LCAO(UnitCell& ucell, const Input_para& inp)
     }
 #endif
 
-    // ESolver_KS_LCAO<TK, TR>::before_all_runners(ucell, inp);
+    ESolver_KS_LCAO<TK, TR>::before_all_runners(ucell, inp);
 
     // initialize rdmft 
+    rdmft_solver.init(this->GG,
+                          this->GK,
+                          this->pv,
+                          ucell,
+                          this->gd,
+                          this->kv,
+                          *(this->pelec),
+                          this->orb_,
+                          this->two_center_bundle_,
+                          inp.dft_functional,
+                          inp.rdmft_power_alpha);
+    
 
 
+    // const integer nbands = static_cast<integer>(inp.nbands > 0 ? inp.nbands : 1);
+    // const integer stiefel_n = nbands;
+    // const integer stiefel_p = nbands;
+
+    nbands_total = inp.nbands;
+    nspin = inp.nspin;
+    nk_total = this->kv.get_nks();
+    nk_total *= nspin;
+
+    nbasis_total = this->psi ->get_nbasis();
+
+    // next set up the manifold
+    integer stiefel_n = nbasis_total;
+    integer stiefel_p = nbands_total;
+
+    integer multitude = nk_total;
+
+    integer numoftypes = 2; // wavefunctions and occupation numbers
+
+
+    this->mani_euc = new Euclidean(nbands_total, "real");
+
+    // create the product manifold
+
+    if(inp.gamma_only)
+    {
+        this->mani_st = new Stiefel(stiefel_n, stiefel_p);
+        this->mani_st->ChooseParamsSet1(); // change later
+        this->mani = new ProductManifold(numoftypes, this->mani_st, multitude, this->mani_euc, multitude);
+    }
+    else
+    {
+        this->mani_cst = new CStiefel(stiefel_n, stiefel_p);
+        this->mani_cst->ChooseParamsSet1(); // change later
+        this->mani = new ProductManifold(numoftypes, this->mani_cst, multitude, this->mani_euc, multitude);
+    }
+
+    // X = this->mani->RandominManifold(); // need to update it, either from KS calculation or file
 
-    const integer nbands = static_cast<integer>(inp.nbands > 0 ? inp.nbands : 1);
-    const integer stiefel_n = nbands;
-    const integer stiefel_p = nbands;
 
-    this->mani = new CStiefel(stiefel_n, stiefel_p);
     this->SetDomain(this->mani);
-    this->mani->ChooseParamsSet4();
-    this->mani->CheckParams();
+    // this->mani->CheckParams();
+
+    // this->mani->ChooseParamsSet4();
 }
 
 // template <typename TK, typename TR>
@@ -75,21 +122,20 @@ RDMFT_LCAO<TK, TR>::RDMFT_LCAO(UnitCell& ucell, const Input_para& inp)
 template <typename TK, typename TR>
 RDMFT_LCAO<TK, TR>::~RDMFT_LCAO()
 {
-    if (this->mani != nullptr)
-    {
-        delete this->mani;
-        this->mani = nullptr;
-    }
+//     if (this->mani != nullptr)
+//     {
+//         delete this->mani;
+//         this->mani = nullptr;
+//     }
 }
 
 
 
 template <typename TK, typename TR>
 realdp RDMFT_LCAO<TK, TR>::f(const Variable& x) const
 {
-    // Energy evaluation for the RDMFT problem will be implemented later.
-    (void)x;
-    return 0.0;
+    this->rdmft_solver.cal_Energy(1); // 1 means not to calculate forces here
+    return this->rdmft_solver.Etotal;
 }
 
 template <typename TK, typename TR>
@@ -101,8 +147,8 @@ Vector& RDMFT_LCAO<TK, TR>::EucGrad(const Variable& x, Vector* result) const
         throw std::runtime_error("RDMFT_LCAO::EucGrad requires a valid result storage");
     }
 
-    *result = x;
-    result->ScalarTimesThis(0.0);
+    *result = x; // copy structure
+    result->ScalarTimesThis(0.0); // zero until real gradient assembly implemented
     return *result;
 }
 
@@ -127,6 +173,7 @@ Vector& RDMFT_LCAO<TK, TR>::EucHessianEta(const Variable& x, const Vector& etax,
 
 
 
+
 // Explicit template instantiations for the typical scalar combinations used in ABACUS.
 template class RDMFT_LCAO<double, double>;
 template class RDMFT_LCAO<std::complex<double>, double>;

source/source_esolver/directmin_problems/rdmft_lcao.h

@@ -2,14 +2,20 @@
 #define RDMFT_LCAO_H
 
 // #include "Problems/CStieBrockett.h"
+#include "Manifolds/Stiefel.h"
 #include "Manifolds/CStiefel.h"
+#include "Manifolds/Euclidean.h"
+#include "Manifolds/MultiManifolds.h"
+
 #include "Problems/Problem.h"
 
 #include <memory>
 
+/*
 #include "source_cell/unitcell.h"
 #include "source_io/module_parameter/input_parameter.h"
 
+
 // for k-points in Brillouin zone
 #include "source_cell/klist.h"
 // for plane-wave basis
@@ -33,6 +39,7 @@
 #include "source_psi/psi.h"
 // for Hamiltonian
 #include "source_hamilt/hamilt.h"
+*/ // Since esolver_ks_lcao.h already includes these headers
 
 #include "source_esolver/esolver_ks_lcao.h"
 
@@ -42,20 +49,40 @@ namespace ModuleESolver
 {
 
 template <typename TK, typename TR>
-class RDMFT_LCAO : public Problem //, public ESolver_KS_LCAO<TK, TR>
+class RDMFT_LCAO : public Problem, public ESolver_KS_LCAO<TK, TR>
 {
 public:
     // RDMFT_LCAO();
     RDMFT_LCAO(UnitCell& ucell, const Input_para& inp);
     ~RDMFT_LCAO();
 
-    virtual realdp f(const Variable &x) const;
+    // Override pure virtuals from ROPTLITE::Problem (must be const signatures)
+    realdp f(const Variable &x) const override;
+    Vector &EucGrad(const Variable &x, Vector *result) const override;
+    Vector &EucHessianEta(const Variable &x, const Vector &etax, Vector *result) const override;
+
+
+    // some variables  for the dimension of the problem, they should be set in the constructor
+    int nbands_total = 0;
+    int nspin = 1;
+    int nk_total = 1;
+    int nbasis_total= 0;
+
+    // int nbasis_local = 0;
 
-    virtual Vector &EucGrad(const Variable &x, Vector *result) const;
-    virtual Vector &EucHessianEta(const Variable &x, const Vector &etax, Vector *result) const;
+    // In this code, we are dealing with joint optimization of both occupation numbers and wavefunctions
+    // So we use multimanifolds, composed of CStiefel or Stiefel manifolds for the orbitals
+    // and Euclidean manifolds for the occupation numbers
+    CStiefel* mani_cst = nullptr; 
+    Stiefel* mani_st = nullptr;
+    Euclidean* mani_euc = nullptr;
+    
 
-    CStiefel* mani = nullptr;
 
+    rdmft::RDMFT<TK, TR> rdmft_solver;
+    
+    ProductManifold* mani = nullptr;
+    // const TwoCenterBundle* two_center_bundle = nullptr;
 
 };
 

source/source_lcao/module_rdmft/rdmft.cpp

@@ -292,7 +292,7 @@ double RDMFT<TK, TR>::cal_E_grad_wfc_occ_num()
 
 // cal_type = 2 just support XC-functional without exx
 template <typename TK, typename TR>
-void RDMFT<TK, TR>::cal_Energy(const int cal_type)
+void RDMFT<TK, TR>::cal_Energy(const int cal_type) const
 {
     double E_Ewald = pelec->f_en.ewald_energy;
     double E_entropy = pelec->f_en.demet;

source/source_lcao/module_rdmft/rdmft.h

@@ -75,8 +75,8 @@ class RDMFT
     psi::Psi<TK> occNum_HamiltWfc;
 
     //! E_RDMFT[4] stores ETV, Ehartree, Exc, Etotal
-    double E_RDMFT[4] = {0.0};
-    double Etotal = 0.0;
+    mutable double E_RDMFT[4] = {0.0};
+    mutable double Etotal = 0.0;
     // std::vector<double> E_RDMFT(4);
 
     //! initialization of rdmft calculation
@@ -103,7 +103,7 @@ class RDMFT
     //! obtain the gradient of total energy with respect to occupation number and wfc
     double cal_E_grad_wfc_occ_num();
 
-    void cal_Energy(const int cal_type = 1);
+    void cal_Energy(const int cal_type = 1) const;
 
     //! update occ_number for optimization algorithms that depend on Hamilton
     void update_occNumber(const ModuleBase::matrix& occ_number_in);
@@ -184,8 +184,8 @@ class RDMFT
     bool exx_spacegroup_symmetry = false;
 #endif
 
-    double etxc = 0.0;
-    double vtxc = 0.0;
+    mutable double etxc = 0.0;
+    mutable double vtxc = 0.0;
     bool only_exx_type = false;
     const int cal_E_type = 1;   // cal_type = 2 just support XC-functional without exx