Some small updates for the format

source/module_elecstate/cal_ux.cpp

@@ -4,62 +4,91 @@
 namespace elecstate {
 
 void cal_ux(UnitCell& ucell) {
+
     if (PARAM.inp.nspin != 4)
     {
         return;
     }
 
-    double amag, uxmod;
+    const double absolute_mag_thr = 1.0e-6;
+
+    double amag = 0.0;
+    double uxmod = 0.0;
+
     int starting_it = 0;
     int starting_ia = 0;
-    bool is_paraller;
+    bool is_paraller = false;
+
     // do not sign feature in teh general case
     ucell.magnet.lsign_ = false;
     ModuleBase::GlobalFunc::ZEROS(ucell.magnet.ux_, 3);
 
-    for (int it = 0; it < ucell.ntype; it++) {
-        for (int ia = 0; ia < ucell.atoms[it].na; ia++) {
+	for (int it = 0; it < ucell.ntype; it++) 
+	{
+		for (int ia = 0; ia < ucell.atoms[it].na; ia++) 
+		{
+            // m_loc_: local magnetization vector for each atom
             amag = pow(ucell.atoms[it].m_loc_[ia].x, 2)
                    + pow(ucell.atoms[it].m_loc_[ia].y, 2)
                    + pow(ucell.atoms[it].m_loc_[ia].z, 2);
-            if (amag > 1e-6) {
-                ucell.magnet.ux_[0] = ucell.atoms[it].m_loc_[ia].x;
-                ucell.magnet.ux_[1] = ucell.atoms[it].m_loc_[ia].y;
-                ucell.magnet.ux_[2] = ucell.atoms[it].m_loc_[ia].z;
-                starting_it = it;
-                starting_ia = ia;
-                ucell.magnet.lsign_ = true;
-                break;
-            }
-        }
-        if (ucell.magnet.lsign_) {
-            break;
-        }
-    }
+
+            // find the first atom (it,ia) whose magnetism is not zero
+            // compute ux
+			if (amag > absolute_mag_thr) 
+			{
+				ucell.magnet.ux_[0] = ucell.atoms[it].m_loc_[ia].x;
+				ucell.magnet.ux_[1] = ucell.atoms[it].m_loc_[ia].y;
+				ucell.magnet.ux_[2] = ucell.atoms[it].m_loc_[ia].z;
+
+				starting_it = it;
+				starting_ia = ia;
+
+				ucell.magnet.lsign_ = true;
+				break;
+			}
+		}
+
+        // if any atom has magnetism, then break the for iteration
+		if (ucell.magnet.lsign_) 
+		{
+			break;
+		}
+	}
+
     // whether the initial magnetizations is parallel
-    for (int it = starting_it; it < ucell.ntype; it++) {
-        for (int ia = 0; ia < ucell.atoms[it].na; ia++) {
-            if (it > starting_it || ia > starting_ia) {
-                ucell.magnet.lsign_
-                    = ucell.magnet.lsign_
-                      && judge_parallel(ucell.magnet.ux_, ucell.atoms[it].m_loc_[ia]);
-            }
-        }
-    }
-    if (ucell.magnet.lsign_) {
-        uxmod = pow(ucell.magnet.ux_[0], 2) + pow(ucell.magnet.ux_[1], 2)
-                + pow(ucell.magnet.ux_[2], 2);
-        if (uxmod < 1e-6) {
-            ModuleBase::WARNING_QUIT("cal_ux", "wrong uxmod");
-        }
-        for (int i = 0; i < 3; i++) {
-            ucell.magnet.ux_[i] *= 1 / sqrt(uxmod);
-        }
-        //       std::cout<<"    Fixed quantization axis for GGA: "
-        //<<std::setw(10)<<ux[0]<<"  "<<std::setw(10)<<ux[1]<<"
-        //"<<std::setw(10)<<ux[2]<<std::endl;
-    }
-    return;
+	for (int it = starting_it; it < ucell.ntype; it++) 
+	{
+		for (int ia = 0; ia < ucell.atoms[it].na; ia++) 
+		{
+			if (it > starting_it || ia > starting_ia) 
+			{
+				ucell.magnet.lsign_
+					= ucell.magnet.lsign_
+					&& judge_parallel(ucell.magnet.ux_, ucell.atoms[it].m_loc_[ia]);
+			}
+		}
+	}
+
+    // if all of the atoms have the same parallel magnetism direction,
+    // then set the direction to a unit vector
+	if (ucell.magnet.lsign_) 
+	{
+		uxmod = pow(ucell.magnet.ux_[0], 2) 
+			+ pow(ucell.magnet.ux_[1], 2)
+			+ pow(ucell.magnet.ux_[2], 2);
+
+		if (uxmod < absolute_mag_thr) 
+		{
+			ModuleBase::WARNING_QUIT("elecstate::cal_ux", "wrong uxmod");
+		}
+
+        // reset the magnetism for each direction
+		for (int i = 0; i < 3; i++) 
+		{
+			ucell.magnet.ux_[i] *= 1 / sqrt(uxmod);
+		}
+	}
+	return;
 }
 
 bool judge_parallel(double a[3], ModuleBase::Vector3<double> b) {
@@ -71,4 +100,4 @@ bool judge_parallel(double a[3], ModuleBase::Vector3<double> b) {
     jp = (fabs(cross) < 1e-6);
     return jp;
 }
-}
+}

source/module_elecstate/magnetism.cpp

@@ -22,6 +22,8 @@ void Magnetism::compute_magnetization(const double& omega,
                                       const double* const * rho, 
                                       double* nelec_spin)
 {
+    assert(nxyz>0);
+
     if (PARAM.inp.nspin==2)
     {
         this->tot_magnetization = 0.00;
@@ -57,36 +59,53 @@ void Magnetism::compute_magnetization(const double& omega,
 	// noncolliear :
 	else if(PARAM.inp.nspin==4)
 	{
-		for(int i=0;i<3;i++) {this->tot_magnetization_nc[i] = 0.00;
-}
+		for(int i=0;i<3;i++) 
+		{
+			this->tot_magnetization_nc[i] = 0.00;
+		}
+
 		this->abs_magnetization = 0.00;
 		for (int ir=0; ir<nrxx; ir++)
 		{
 			double diff = sqrt(pow(rho[1][ir], 2) + pow(rho[2][ir], 2) +pow(rho[3][ir], 2));
 
-			for(int i=0;i<3;i++) {this->tot_magnetization_nc[i] += rho[i+1][ir];
-}
+			for(int i=0;i<3;i++) 
+			{
+				this->tot_magnetization_nc[i] += rho[i+1][ir];
+			}
 			this->abs_magnetization += std::abs(diff);
 		}
 #ifdef __MPI
         Parallel_Reduce::reduce_pool(this->tot_magnetization_nc, 3);
         Parallel_Reduce::reduce_pool(this->abs_magnetization);
 #endif
-		for(int i=0;i<3;i++) {this->tot_magnetization_nc[i] *= omega/ nxyz;
-}
+		for(int i=0;i<3;i++) 
+		{
+			this->tot_magnetization_nc[i] *= omega/ nxyz;
+		}
+
 		this->abs_magnetization *= omega/ nxyz;
-		GlobalV::ofs_running<<"total magnetism (Bohr mag/cell)"<<'\t'<<this->tot_magnetization_nc[0]<<'\t'<<this->tot_magnetization_nc[1]<<'\t'<<this->tot_magnetization_nc[2]<<'\n';
+
+		GlobalV::ofs_running<<"total magnetism (Bohr mag/cell)"
+			<<'\t'<<this->tot_magnetization_nc[0]
+			<<'\t'<<this->tot_magnetization_nc[1]
+			<<'\t'<<this->tot_magnetization_nc[2]<<'\n';
+
 		ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"absolute magnetism (Bohr mag/cell)",this->abs_magnetization);
 	}
 
     return;
 }
 
+
 bool Magnetism::judge_parallel(double a[3], ModuleBase::Vector3<double> b)
 {
    bool jp=false;
-   double cross;
-   cross = pow((a[1]*b.z-a[2]*b.y),2) +  pow((a[2]*b.x-a[0]*b.z),2) + pow((a[0]*b.y-a[1]*b.x),2);
+   double cross=0.0;
+   cross = pow((a[1]*b.z-a[2]*b.y),2) 
+	   + pow((a[2]*b.x-a[0]*b.z),2) 
+	   + pow((a[0]*b.y-a[1]*b.x),2);
+
    jp = (fabs(cross)<1e-6);
    return jp;
 }

source/module_elecstate/magnetism.h

@@ -13,7 +13,7 @@ class Magnetism
     Magnetism();
     ~Magnetism();
 
-    // notice : becast is done in unitcell
+    // notice : bcast (MPI operation) is done in unitcell
     double *start_magnetization=nullptr;
 
     // tot_magnetization : majority spin - minority spin (nelup - neldw).
@@ -28,10 +28,9 @@ class Magnetism
 
                                double* nelec_spin = nullptr);
 
-    ModuleBase::Vector3<double> *m_loc_=nullptr;//magnetization for each element along c-axis
-	double *angle1_=nullptr;                    //angle between c-axis and real spin std::vector
-	double *angle2_=nullptr;                    //angle between a-axis and real spin std::vector projection in ab-plane
-    //void cal_ux(const int ntype);
+    ModuleBase::Vector3<double> *m_loc_=nullptr; //magnetization for each element along c-axis
+	double *angle1_=nullptr;                     //angle between c-axis and real spin std::vector
+	double *angle2_=nullptr;                     //angle between a-axis and real spin std::vector projection in ab-plane
     double ux_[3]={0.0};
 	bool lsign_=false;
 

source/module_esolver/esolver_fp.cpp

@@ -258,6 +258,7 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep)
 {
     ModuleBase::TITLE("ESolver_FP", "before_scf");
 
+    // if the cell has changed
     if (ucell.cell_parameter_updated)
     {
         // only G-vector and K-vector are changed due to the change of lattice
@@ -266,62 +267,71 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep)
         this->pw_rho->collect_local_pw();
         this->pw_rho->collect_uniqgg();
 
+        // if double grid used in USPP, update related quantities in dense grid
         if (PARAM.globalv.double_grid)
         {
             this->pw_rhod->initgrids(ucell.lat0, ucell.latvec, pw_rhod->nx, pw_rhod->ny, pw_rhod->nz);
             this->pw_rhod->collect_local_pw();
             this->pw_rhod->collect_uniqgg();
         }
 
+        // reset local pseudopotentials
         this->locpp.init_vloc(ucell, this->pw_rhod);
         ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "LOCAL POTENTIAL");
 
         this->pelec->omega = ucell.omega;
 
+        // perform symmetry analysis
         if (ModuleSymmetry::Symmetry::symm_flag == 1)
         {
             ucell.symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, GlobalV::ofs_running);
             ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SYMMETRY");
         }
 
+        // reset k-points
         kv.set_after_vc(PARAM.inp.nspin, ucell.G, ucell.latvec);
         ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT K-POINTS");
     }
 
+    //----------------------------------------------------------
     // charge extrapolation
+    //----------------------------------------------------------
     if (ucell.ionic_position_updated)
     {
         this->CE.update_all_dis(ucell);
-        this->CE.extrapolate_charge(&(this->Pgrid),
+        this->CE.extrapolate_charge(&this->Pgrid,
                                     ucell,
                                     &this->chr,
-                                    &(this->sf),
+                                    &this->sf,
                                     GlobalV::ofs_running,
                                     GlobalV::ofs_warning);
     }
 
     //----------------------------------------------------------
-    // about vdw, jiyy add vdwd3 and linpz add vdwd2
+    //! calculate D2 or D3 vdW
     //----------------------------------------------------------
     auto vdw_solver = vdw::make_vdw(ucell, PARAM.inp, &(GlobalV::ofs_running));
     if (vdw_solver != nullptr)
     {
         this->pelec->f_en.evdw = vdw_solver->get_energy();
     }
 
-    // calculate ewald energy
+    //----------------------------------------------------------
+    //! calculate ewald energy
+    //----------------------------------------------------------
     if (!PARAM.inp.test_skip_ewald)
     {
         this->pelec->f_en.ewald_energy = H_Ewald_pw::compute_ewald(ucell, this->pw_rhod, this->sf.strucFac);
     }
 
     //----------------------------------------------------------
-    //! cal_ux should be called before init_scf because
-    //! the direction of ux is used in noncoline_rho
+    //! set direction of magnetism, used in non-collinear case 
     //----------------------------------------------------------
     elecstate::cal_ux(ucell);
 
+    //----------------------------------------------------------
     //! output the initial charge density
+    //----------------------------------------------------------
     if (PARAM.inp.out_chg[0] == 2)
     {
         for (int is = 0; is < PARAM.inp.nspin; is++)
@@ -339,7 +349,9 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep)
         }
     }
 
+    //----------------------------------------------------------
     //! output total local potential of the initial charge density
+    //----------------------------------------------------------
     if (PARAM.inp.out_pot == 3)
     {
         for (int is = 0; is < PARAM.inp.nspin; is++)

source/module_esolver/esolver_ks_lcao.cpp

@@ -160,7 +160,9 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_pa
 #ifdef __EXX
     // 5) initialize exx
     // PLEASE simplify the Exx_Global interface
-    if (PARAM.inp.calculation == "scf" || PARAM.inp.calculation == "relax" || PARAM.inp.calculation == "cell-relax"
+    if (PARAM.inp.calculation == "scf" 
+        || PARAM.inp.calculation == "relax" 
+        || PARAM.inp.calculation == "cell-relax"
         || PARAM.inp.calculation == "md")
     {
         if (GlobalC::exx_info.info_global.cal_exx)