fix final mag in STRU_ION_D (#4044)

deepmodeling · Apr 23, 2024 · 2524477 · 2524477
1 parent f709ba0
commit 2524477
Show file tree

Hide file tree

Showing 3 changed files with 35 additions and 13 deletions.
diff --git a/source/module_cell/read_atoms.cpp b/source/module_cell/read_atoms.cpp
@@ -1003,6 +1003,7 @@ void UnitCell::print_stru_file(const std::string &fn, const int &type, const int
 
 	if(type==1)
 	{
+		int nat_tmp = 0;
 		ofs << "Cartesian" << std::endl;
 		for(int it=0; it<ntype; it++)
 		{
@@ -1028,28 +1029,31 @@ void UnitCell::print_stru_file(const std::string &fn, const int &type, const int
                     ofs << context.str();
                 }
 
-                if (GlobalV::NSPIN == 2)
+                if (GlobalV::NSPIN == 2 && GlobalV::out_mul)
                 {
                     // output magnetic information
                     ofs << " mag ";
                     context.set_context("double_w6_f2");
-                    context << atoms[it].mag[ia];
+                    context << atom_mulliken[nat_tmp][1];
                     ofs << context.str();
                 }
-                else if (GlobalV::NSPIN == 4)
+                else if (GlobalV::NSPIN == 4 && GlobalV::out_mul)
                 {
                     // output magnetic information
                     ofs << " mag ";
-                    context.set_context("vector3d");
-                    context << atoms[it].m_loc_[ia].x << " " << atoms[it].m_loc_[ia].y << " " << atoms[it].m_loc_[ia].z;
-                    ofs << context.str();
+					ofs << std::fixed << std::setprecision(5);
+                    ofs << std::setw(8) << atom_mulliken[nat_tmp][1] 
+						<< std::setw(8) << atom_mulliken[nat_tmp][2] 
+						<< std::setw(8) << atom_mulliken[nat_tmp][3];
                 }
                 ofs << std::endl;
+				nat_tmp++;
 			}
 		}
 	}
 	else if(type==2)
 	{
+		int nat_tmp = 0;
 		ofs << "Direct" << std::endl;
 		for(int it=0; it<ntype; it++)
 		{
@@ -1075,23 +1079,25 @@ void UnitCell::print_stru_file(const std::string &fn, const int &type, const int
                     ofs << context.str();
                 }
 
-                if (GlobalV::NSPIN == 2)
+                if (GlobalV::NSPIN == 2 && GlobalV::out_mul)
                 {
                     // output magnetic information
                     ofs << " mag ";
                     context.set_context("double_w6_f2");
-                    context << atoms[it].mag[ia];
+                    context << atom_mulliken[nat_tmp][1];
                     ofs << context.str();
                 }
-                else if (GlobalV::NSPIN == 4)
+                else if (GlobalV::NSPIN == 4 && GlobalV::out_mul)
                 {
                     // output magnetic information
                     ofs << " mag ";
-                    context.set_context("vector3d");
-                    context << atoms[it].m_loc_[ia].x << " " << atoms[it].m_loc_[ia].y << " " << atoms[it].m_loc_[ia].z;
-                    ofs << context.str();
+					ofs << std::fixed << std::setprecision(5);
+                    ofs << std::setw(8) << atom_mulliken[nat_tmp][1] 
+						<< std::setw(8) << atom_mulliken[nat_tmp][2] 
+						<< std::setw(8) << atom_mulliken[nat_tmp][3];
                 }
                 ofs << std::endl;
+				nat_tmp++;
 			}
 		}
 	}

diff --git a/source/module_cell/unitcell.h b/source/module_cell/unitcell.h
@@ -23,7 +23,7 @@ class UnitCell
     Magnetism magnet;  // magnetism Yu Liu 2021-07-03
     void cal_ux();
     bool judge_parallel(double a[3],ModuleBase::Vector3<double> b);
-	double *atom_mag;
+    std::vector<std::vector<double>> atom_mulliken;  //[nat][nspin]
 	int n_mag_at;
 
     std::string& Coordinate = lat.Coordinate;

diff --git a/source/module_io/mulliken_charge.cpp b/source/module_io/mulliken_charge.cpp
@@ -312,6 +312,16 @@ void ModuleIO::out_mulliken(const int& step, LCAO_Matrix* LM, const elecstate::E
 		os << "Decomposed Mulliken populations" << std::endl;
         GlobalV::ofs_running << std::endl <<std::endl;
 
+        // alloacte space for ucell.atom_mulliken
+        if (GlobalC::ucell.atom_mulliken.empty())
+        {
+            GlobalC::ucell.atom_mulliken.resize(GlobalC::ucell.nat);
+            for (int iat = 0; iat < GlobalC::ucell.nat; iat++)
+            {
+                GlobalC::ucell.atom_mulliken[iat].resize(GlobalV::NSPIN);
+            }
+        }
+
         for (size_t i = 0; i != GlobalC::ucell.nat; ++i)
         {
             double total_charge = 0.0, atom_mag = 0.0;
@@ -428,6 +438,8 @@ void ModuleIO::out_mulliken(const int& step, LCAO_Matrix* LM, const elecstate::E
                 os << "Total Charge on atom:  " << GlobalC::ucell.atoms[t].label <<  std::setw(20) << total_charge <<std::endl;
                 os << "Total Magnetism on atom:  " << GlobalC::ucell.atoms[t].label <<  std::setw(20) << ModuleIO::output_cut(atom_mag) <<std::endl;
                 GlobalV::ofs_running << "Total Magnetism on atom:  " << GlobalC::ucell.atoms[t].label <<  std::setw(20) << std::setprecision(10) << atom_mag <<std::endl;
+                GlobalC::ucell.atom_mulliken[i][0] = total_charge;
+                GlobalC::ucell.atom_mulliken[i][1] = atom_mag;
             }
             else if (GlobalV::NSPIN==4)
             {
@@ -443,6 +455,10 @@ void ModuleIO::out_mulliken(const int& step, LCAO_Matrix* LM, const elecstate::E
                 GlobalV::ofs_running << "Total Magnetism on atom:  " << GlobalC::ucell.atoms[t].label <<  std::setw(10)
                 << "("  << std::setprecision(10) << spin2 << ", " << std::setprecision(10) << spin3 << ", " << std::setprecision(10) << spin4 << ")"
                 <<std::endl;
+                GlobalC::ucell.atom_mulliken[i][0] = spin1;
+                GlobalC::ucell.atom_mulliken[i][1] = spin2;
+                GlobalC::ucell.atom_mulliken[i][2] = spin3;
+                GlobalC::ucell.atom_mulliken[i][3] = spin4;
             }
             os << std::endl <<std::endl;
         }