<Feature> In order to check the validity of MPN KEDF, calculate and output TF KEDF when running MPN KEDF.

Author: sunliang98

source/module_esolver/esolver_fp.cpp

@@ -52,7 +52,7 @@ void ESolver_FP::init(Input& inp, UnitCell& cell)
 #ifdef __MPI
 	this->pw_rho->initmpi(GlobalV::NPROC_IN_POOL, GlobalV::RANK_IN_POOL, POOL_WORLD);
 #endif
-	if (this->classname == "ESolver_OF")
+	if (this->classname == "ESolver_OF" || GlobalV::of_ml_gene_data == 1)
     {
 		this->pw_rho->setfullpw(inp.of_full_pw, inp.of_full_pw_dim);
     }

source/module_esolver/esolver_of.cpp

@@ -584,6 +584,35 @@ void ESolver_OF::after_opt(const int istep, UnitCell& ucell)
             &(ucell),
             this->pelec->pot->get_fixed_v());
     }
+    if (this->of_kinetic_ == "mpn" || this->of_kinetic_ == "ml")
+    {
+        this->tf_->get_energy(this->pelec->charge->rho);
+        std::cout << "ML Term = " << this->ml_->ml_energy << " Ry, TF Term = " << this->tf_->tf_energy << " Ry." << std::endl;
+        if (this->ml_->ml_energy >= this->tf_->tf_energy)
+        {
+            std::cout << "WARNING: ML >= TF" << std::endl;
+        }
+    }
+
+    if (GlobalV::of_ml_gene_data)
+    {
+        this->pelec->pot->update_from_charge(pelec->charge, &GlobalC::ucell); // Hartree + XC + external
+        this->kinetic_potential(pelec->charge->rho, this->pphi_, this->pelec->pot->get_effective_v()); // (kinetic + Hartree + XC + external) * 2 * phi
+        
+        const double* vr_eff = this->pelec->pot->get_effective_v(0);
+        for (int ir = 0; ir < this->pw_rho->nrxx; ++ir)
+        {
+            this->pdEdphi_[0][ir] = vr_eff[ir];
+        }
+        this->mu_[0] = this->cal_mu(this->pphi_[0], this->pdEdphi_[0], this->nelec_[0]);
+
+        // === temporary ===
+        // assert(GlobalV::of_kinetic == "wt" || GlobalV::of_kinetic == "ml");
+        // =================
+        std::cout << "Generating Training data..." << std::endl;
+        std::cout << "mu = " << this->mu_[0] << std::endl;
+        this->ml_->generateTrainData(pelec->charge->rho, *(this->wt_), *(this->tf_), this->pw_rho, vr_eff);
+    }
 }
 
 /**

source/module_esolver/esolver_of_interface.cpp

@@ -12,7 +12,11 @@ namespace ModuleESolver
 void ESolver_OF::init_kedf(Input& inp)
 {
     //! Thomas-Fermi (TF) KEDF, TF+ KEDF, and Want-Teter (WT) KEDF
-    if (this->of_kinetic_ == "tf" || this->of_kinetic_ == "tf+" || this->of_kinetic_ == "wt")
+    if (this->of_kinetic_ == "tf"
+     || this->of_kinetic_ == "tf+"
+     || this->of_kinetic_ == "wt"
+     || this->of_kinetic_ == "ml"
+     || this->of_kinetic_ == "mpn")
     {
 		if (this->tf_ == nullptr)
 		{
@@ -103,6 +107,7 @@ void ESolver_OF::kinetic_potential(double** prho, double** pphi, ModuleBase::mat
     if (this->of_kinetic_ == "ml" || this->of_kinetic_ == "mpn")
     {
         this->ml_->ml_potential(prho, this->pw_rho, rpot);
+        this->tf_->get_energy(prho); // temp
     }
 
     // Before call vw_potential, change rpot to rpot * 2 * pphi
@@ -164,6 +169,11 @@ double ESolver_OF::kinetic_energy()
     if (this->of_kinetic_ == "ml" || this->of_kinetic_ == "mpn")
     {
         kinetic_energy += this->ml_->ml_energy;
+        if (this->ml_->ml_energy >= this->tf_->tf_energy)
+        {
+            std::cout << "WARNING: ML >= TF" << std::endl;
+            std::cout << "ML Term = " << this->ml_->ml_energy << " Ry, TF Term = " << this->tf_->tf_energy << " Ry." << std::endl;
+        }
     }
 
     return kinetic_energy;

source/module_hamilt_pw/hamilt_ofdft/kedf_ml_pot.cpp

@@ -43,6 +43,7 @@ void KEDF_ML::get_potential_(const double * const * prho, ModulePW::PW_Basis *pw
                     //   + xinlterm[ir] + tanhxinlterm[ir] + tanhxi_nlterm[ir]
                     //   + tanhptanh_pnlterm[ir] + tanhqtanh_qnlterm[ir]
                     //   + tanhptanhp_nlterm[ir] + tanhqtanhq_nlterm[ir];
+        // if (pauli_potential[ir] < 0) pauli_potential[ir] = 0;
         rpotential(0, ir) += pauli_potential[ir];
 
         // if (this->enhancement_cpu_ptr[ir] < 0)