<fix>

	Now out_chg is the parameter deciding if to output charge and out_freq_elec and out_freq_ion decide the frequency.

Author: Qianruipku

source/input.cpp

@@ -252,6 +252,8 @@ void Input::Default(void)
     printe = 100; // must > 0
     init_chg = "atomic";
     chg_extrap = "atomic"; // xiaohui modify 2015-02-01
+    out_freq_elec = 0;
+    out_freq_ion = 0;
     out_chg = 0;
     out_dm = 0;
 
@@ -932,6 +934,14 @@ bool Input::Read(const std::string &fn)
         {
             read_value(ifs, chg_extrap); // xiaohui modify 2015-02-01
         }
+        else if (strcmp("out_freq_elec", word) == 0)
+        {
+            read_value(ifs, out_freq_elec);
+        }
+        else if (strcmp("out_freq_ion", word) == 0)
+        {
+            read_value(ifs, out_freq_ion);
+        }
         else if (strcmp("out_chg", word) == 0)
         {
             read_value(ifs, out_chg);
@@ -1914,6 +1924,8 @@ void Input::Bcast()
     Parallel_Common::bcast_int(printe);
     Parallel_Common::bcast_string(init_chg);
     Parallel_Common::bcast_string(chg_extrap); // xiaohui modify 2015-02-01
+    Parallel_Common::bcast_int(out_freq_elec);
+    Parallel_Common::bcast_int(out_freq_ion);
     Parallel_Common::bcast_int(out_chg);
     Parallel_Common::bcast_int(out_dm);
 

source/input.h

@@ -197,7 +197,9 @@ class Input
     int mem_saver; // 1: save psi when nscf calculation.
 
     int printe; // mohan add 2011-03-16
-    int out_chg; // output charge density.
+    int out_freq_elec;  // the frequency ( >= 0) of electronic iter to output charge density and wavefunction. 0: output only when converged
+    int out_freq_ion;  // the frequency ( >= 0 ) of ionic step to output charge density and wavefunction. 0: output only when ion steps are finished
+    int out_chg; // output charge density. 0: no; 1: yes
     int out_dm; // output density matrix.
     int out_pot; // yes or no
     int out_wfc_pw; // 0: no; 1: txt; 2: dat

source/input_conv.cpp

@@ -423,6 +423,8 @@ void Input_Conv::Convert(void)
     //----------------------------------------------------------
     // wavefunction / charge / potential / (2/4)
     //----------------------------------------------------------
+    GlobalV::OUT_FREQ_ELEC = INPUT.out_freq_elec;
+    GlobalV::OUT_FREQ_ION = INPUT.out_freq_ion;
     GlobalC::pot.init_chg = INPUT.init_chg;
     GlobalC::pot.chg_extrap = INPUT.chg_extrap; // xiaohui modify 2015-02-01
     GlobalC::CHR.out_chg = INPUT.out_chg;

source/module_base/global_variable.cpp

@@ -43,6 +43,8 @@ double PRESS3 = 0.0;
 double PRESSURE = 0.0;
 std::string RELAX_METHOD = "bfgs";
 std::string OUT_LEVEL = "ie";
+int OUT_FREQ_ELEC = 0;
+int OUT_FREQ_ION = 0;
 int RELAX_NMAX = 20;
 int MD_NSTEP = 20;
 int SCF_NMAX = 50;

source/module_base/global_variable.h

@@ -44,6 +44,8 @@ extern double PRESS3;
 extern double PRESSURE;
 extern std::string RELAX_METHOD;
 extern std::string OUT_LEVEL;
+extern int OUT_FREQ_ELEC;
+extern int OUT_FREQ_ION;
 
 extern int RELAX_NMAX; // 8.3
 extern int SCF_NMAX; // 8.4

source/module_esolver/esolver_ks.cpp

@@ -15,30 +15,31 @@ namespace ModuleESolver
 
 ESolver_KS::ESolver_KS()
 {
+    classname = "ESolver_KS";
     basisname = "PLEASE ADD BASISNAME FOR CURRENT ESOLVER.";
     diag_ethr = GlobalV::PW_DIAG_THR; 
     scf_thr = GlobalV::SCF_THR; 
     drho = 0.0;
     maxniter = GlobalV::SCF_NMAX;
     niter = maxniter;
-
-
+    out_freq_elec = GlobalV::OUT_FREQ_ELEC;
 }
 
 void ESolver_KS:: hamilt2density(int istep, int iter, double ethr)
 {
+    ModuleBase::timer:: tick(this->classname,"hamilt2density");
     //Temporarily, before HSolver is constructed, it should be overrided by
     //LCAO, PW, SDFT and TDDFT.
     //After HSolver is constructed, LCAO, PW, SDFT should delete their own
     //hamilt2density() and use:
     //this->phsol->solve(this->phamilt, this->pes, this->wf, ETHR);
-    
+    ModuleBase::timer:: tick(this->classname,"hamilt2density");
 }
 
 
 void ESolver_KS:: Run(int istep, UnitCell_pseudo& cell)
 {
-    ModuleBase::timer:: tick("ESolver_KS","run");
+    ModuleBase::timer:: tick(this->classname,"Run");
 
     this->printhead(); //print the headline on the screen.
     this->beforeiter(); //Something else to do before the iter loop
@@ -67,14 +68,7 @@ void ESolver_KS:: Run(int istep, UnitCell_pseudo& cell)
             // The error of HSolver is larger than drho, so a more precise HSolver should be excuted.
             if(hsover_error > drho)  
             {
-                GlobalV::ofs_running << " Notice: Threshold on eigenvalues was too large.\n";
-                ModuleBase::WARNING("scf","Threshold on eigenvalues was too large.");
-                GlobalV::ofs_running << " hsover_error=" << hsover_error << " > drho=" << drho << std::endl;
-
-                GlobalV::ofs_running << " Origin diag_ethr = " << this->diag_ethr << std::endl;
-                this->diag_ethr =0.1 *  drho / GlobalC::CHR.nelec; 
-                GlobalV::ofs_running << " New    diag_ethr = " << this->diag_ethr << std::endl;
-                
+                reset_diagethr(GlobalV::ofs_running, hsover_error);
                 this->hamilt2density(istep, iter, this->diag_ethr);
                 drho = GlobalC::CHR.get_drho();
             }   
@@ -92,7 +86,7 @@ void ESolver_KS:: Run(int istep, UnitCell_pseudo& cell)
         // Hamilt should be used after it is constructed.
         // this->phamilt->update(conv_elec);
         updatepot(conv_elec);
-        eachiterfinish(iter);  
+        eachiterfinish(iter,conv_elec);  
         iterend = std::clock();
         double duration = double(iterend-iterstart) / CLOCKS_PER_SEC;
         printiter(conv_elec, iter, drho, duration, diag_ethr);
@@ -104,11 +98,11 @@ void ESolver_KS:: Run(int istep, UnitCell_pseudo& cell)
     }
     afteriter(conv_elec); 
 
-    ModuleBase::timer:: tick("ESolver_KS","run");
+    ModuleBase::timer:: tick(this->classname,"Run");
     return;
 };
 
-
+//<Temporary> It should be a function of Diag_H class in the future.
 void ESolver_KS:: set_ethr(int istep, int iter)
 {
 //It is too complex now and should be modified.
@@ -191,6 +185,16 @@ void ESolver_KS:: writehead(std::ofstream &ofs_running, int istep, int iter)
         << "--------------------------------\n";
 }
 
+void ESolver_KS:: reset_diagethr(std::ofstream &ofs_running, double hsover_error)
+{
+    ofs_running << " Notice: Threshold on eigenvalues was too large.\n";
+    ModuleBase::WARNING("scf","Threshold on eigenvalues was too large.");
+    ofs_running << " hsover_error=" << hsover_error << " > DRHO=" << drho << std::endl;
+    ofs_running << " Origin diag_ethr = " << this->diag_ethr << std::endl;
+    this->diag_ethr =0.1 *  drho / GlobalC::CHR.nelec; 
+    ofs_running << " New    diag_ethr = " << this->diag_ethr << std::endl;
+}
+
 int ESolver_KS:: getniter()
 {
     return this->niter;

source/module_esolver/esolver_ks.h

@@ -18,6 +18,7 @@ class ESolver_KS: public ESolver_FP
         double drho;      // the difference between rho_in (before HSolver) and rho_out (After HSolver)
         int maxniter;     // maximum iter steps for scf
         int niter;        // iter steps actually used in scf
+        int out_freq_elec;// frequency for output
 
         virtual void Run(int istep, UnitCell_pseudo& cell) override;
 
@@ -32,7 +33,7 @@ class ESolver_KS: public ESolver_FP
         // Something to do before hamilt2density function in each iter loop.
         virtual void eachiterinit(int iter){}; 
         // Something to do after hamilt2density function in each iter loop.
-        virtual void eachiterfinish(int iter){}; 
+        virtual void eachiterfinish(int iter, bool conv){}; 
         // Something to do after the iter loop when scf is converged or comes to the max iter step.
         virtual void afteriter(bool){};
         // <Temporary> It should be replaced by a function in Hamilt Class
@@ -44,14 +45,15 @@ class ESolver_KS: public ESolver_FP
         // Set ethr for hsolver
         void set_ethr(int istep, int iter);
         // Print the headline on the screen:
-        // ITER   ETOT(eV)       EDIFF(eV)      SCF_THR    TIME(s) 
+        // ITER   ETOT(eV)       EDIFF(eV)      DRHO    TIME(s) 
         void printhead();
         // Print inforamtion in each iter
         // G1    -3.435545e+03  0.000000e+00   3.607e-01  2.862e-01
         void printiter(bool conv, int iter, double drho, double duration, double ethr);
         // Write the headline in the running_log file
         // "PW/LCAO" ALGORITHM --------------- ION=   1  ELEC=   1--------------------------------
         void writehead(std::ofstream &ofs_running,int istep, int iter);
+        void reset_diagethr(std::ofstream &ofs_running, double hsover_error);
 
 
 

source/module_esolver/esolver_ks_pw.cpp

@@ -305,45 +305,55 @@ void ESolver_KS_PW:: updatepot(bool conv_elec)
     GlobalC::pot.set_vr_eff();
 }
 
-void ESolver_KS_PW:: eachiterfinish(int iter)
+void ESolver_KS_PW:: eachiterfinish(int iter, bool conv_elec)
 {
     //print_eigenvalue(GlobalV::ofs_running);
     GlobalC::en.calculate_etot();
-    if(GlobalC::CHR.out_chg != 0 && iter % GlobalC::CHR.out_chg == 0)
+    //We output it for restarting the scf.
+    bool print = false;
+    if(this->out_freq_elec == 0 )
     {
-        for(int is=0; is<GlobalV::NSPIN; is++)
+        if(conv_elec) print = true;
+    }
+    else
+    {
+        if(iter % this->out_freq_elec == 0 || conv_elec) print = true;
+    }
+    
+    if(print)
+    {
+        if(GlobalC::CHR.out_chg > 0)
+        {
+            for(int is=0; is<GlobalV::NSPIN; is++)
+            {
+                std::stringstream ssc;
+                std::stringstream ss1;
+                ssc << GlobalV::global_out_dir << "tmp" << "_SPIN" << is + 1 << "_CHG";
+                GlobalC::CHR.write_rho(GlobalC::CHR.rho_save[is], is, iter, ssc.str(), 3);//mohan add 2007-10-17
+	        	ss1 << GlobalV::global_out_dir << "tmp" << "_SPIN" << is + 1 << "_CHG.cube";
+	        	GlobalC::CHR.write_rho_cube(GlobalC::CHR.rho_save[is], is, ssc.str(), 3);
+            }
+        }
+        //output wavefunctions
+        if(GlobalC::wf.out_wfc_pw == 1 || GlobalC::wf.out_wfc_pw == 2)
         {
-            std::stringstream ssc;
-            std::stringstream ss1;
-            ssc << GlobalV::global_out_dir << "tmp" << "_SPIN" << is + 1 << "_CHG";
-            GlobalC::CHR.write_rho(GlobalC::CHR.rho_save[is], is, iter, ssc.str(), 3);//mohan add 2007-10-17
-	    	ss1 << GlobalV::global_out_dir << "tmp" << "_SPIN" << is + 1 << "_CHG.cube";
-	    	GlobalC::CHR.write_rho_cube(GlobalC::CHR.rho_save[is], is, ssc.str(), 3);
+            std::stringstream ssw;
+            ssw << GlobalV::global_out_dir << "WAVEFUNC";
+            //WF_io::write_wfc( ssw.str(), GlobalC::wf.evc );
+            // mohan update 2011-02-21
+	    	//qianrui update 2020-10-17
+            WF_io::write_wfc2( ssw.str(), GlobalC::wf.evc, GlobalC::pw.gcar);
+            //ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running,"write wave functions into file WAVEFUNC.dat");
         }
+        
     }
 
+    
 }
 
 
 void ESolver_KS_PW::afteriter(bool conv_elec)
 {
-    //output wavefunctions
-    if(GlobalC::wf.out_wfc_pw == 1 || GlobalC::wf.out_wfc_pw == 2)
-    {
-        std::stringstream ssw;
-        ssw << GlobalV::global_out_dir << "WAVEFUNC";
-        //WF_io::write_wfc( ssw.str(), GlobalC::wf.evc );
-        // mohan update 2011-02-21
-		//qianrui update 2020-10-17
-        WF_io::write_wfc2( ssw.str(), GlobalC::wf.evc, GlobalC::pw.gcar);
-        //ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running,"write wave functions into file WAVEFUNC.dat");
-    }
-
-
-    //--------------------------------------
-    // output charge density for converged,
-    // 0 means don't need to consider iter,
-    //--------------------------------------
 #ifdef __LCAO
     if(GlobalC::chi0_hilbert.epsilon)                 // pengfei 2016-11-23
     {

source/module_esolver/esolver_ks_pw.h

@@ -24,7 +24,7 @@ class ESolver_KS_PW: public ESolver_KS
     virtual void eachiterinit(int iter) override; 
     virtual void hamilt2density(int istep, int iter, double ethr) override;
     virtual void updatepot(bool conv) override;
-    virtual void eachiterfinish(int iter) override; 
+    virtual void eachiterfinish(int iter, bool conv) override; 
     virtual void afteriter(bool) override;
 
     // <Temporary> Get wavefunctions and eigen energies.