Refactor: replace raw pointers to vector in atom_spec.h (#5521)

* Refactor: new real space projection DeltaSpin method

* Refactor: replace raw pointers to vector in atom_spec.h

* Fix: add namespace

* Fix: save mass

* [pre-commit.ci lite] apply automatic fixes

* Fix: error in CI test

* Fix: error in CI

* [pre-commit.ci lite] apply automatic fixes

---------

Co-authored-by: dyzheng <[email protected]>
Co-authored-by: pre-commit-ci-lite[bot] <117423508+pre-commit-ci-lite[bot]@users.noreply.github.com>

Author: 3 people

docs/advanced/input_files/input-main.md

@@ -2925,7 +2925,7 @@ These variables are used to control DFT+U correlated parameters
 
   - where $\gamma$ is a parameter that adjusts the relative weight of the error function to the derivative error function.
 - **Unit**: Bohr
-- **Default**: 5.0
+- **Default**: 3.0
 
 [back to top](#full-list-of-input-keywords)
 

source/Makefile.Objects

@@ -317,7 +317,7 @@ OBJS_HAMILT_LCAO=hamilt_lcao.o\
     op_dftu_lcao.o\
     deepks_lcao.o\
     op_exx_lcao.o\
-    sc_lambda_lcao.o\
+    dspin_lcao.o\
     dftu_lcao.o\
 
 OBJS_HCONTAINER=base_matrix.o\
@@ -696,14 +696,11 @@ OBJS_DFTU=dftu.o\
       dftu_hamilt.o
 
 OBJS_DELTASPIN=basic_funcs.o\
-      cal_h_lambda.o\
       cal_mw_from_lambda.o\
-      cal_mw_helper.o\
       cal_mw.o\
       init_sc.o\
       lambda_loop_helper.o\
       lambda_loop.o\
-      sc_parse_json.o\
       spin_constrain.o\
       template_helpers.o\
 

source/module_cell/atom_spec.cpp

@@ -13,55 +13,20 @@ Atom::Atom()
     type = 0;
     stapos_wf = 0;
     mass = 0.0;
-    tau = new ModuleBase::Vector3<double>[1];
-    dis = new ModuleBase::Vector3<double>[1];
-    taud = new ModuleBase::Vector3<double>[1];
-    vel = new ModuleBase::Vector3<double>[1];
-    mag = new double[1];
-    angle1 = new double[1];
-    angle2 = new double[1];
-    m_loc_ = new ModuleBase::Vector3<double>[1];
-    l_nchi = new int[1];
-    iw2l = new int[1];
-    iw2n = new int[1];
-    iw2m = new int[1];
-    iw2_ylm = new int[1];
-    iw2_new = new bool[1];
-    mbl = new ModuleBase::Vector3<int>[1];
 }
 
 Atom::~Atom()
 {
-    delete[] tau;
-    delete[] dis;
-    delete[] taud;
-    delete[] vel;
-    delete[] mag;
-    delete[] angle1;
-    delete[] angle2;
-    delete[] m_loc_;
-    delete[] l_nchi;
-    delete[] iw2l;
-    delete[] iw2n;
-    delete[] iw2m;
-    delete[] iw2_ylm;
-    delete[] iw2_new;
-    delete[] mbl;
 }
 
 void Atom::set_index(void)
 {
     assert(nw != 0);
-    delete[] iw2l;
-    delete[] iw2n;
-    delete[] iw2m;
-    delete[] iw2_ylm;
-    delete[] iw2_new;
-    iw2l = new int[nw];
-    iw2n = new int[nw];
-    iw2m = new int[nw];
-    iw2_ylm = new int[nw];
-    iw2_new = new bool[nw];
+    this->iw2l.resize(nw, 0);
+    this->iw2n.resize(nw, 0);
+    this->iw2m.resize(nw, 0);
+    this->iw2_ylm.resize(nw, 0);
+    this->iw2_new.resize(nw, false); // bool array to check if the local orbital is new
 
     int index = 0;
     for (int L = 0; L <= nwl; L++)
@@ -103,7 +68,7 @@ void Atom::print_Atom(std::ofstream& ofs)
     ModuleBase::GlobalFunc::OUT(ofs, "mass", mass);
     ofs << std::endl;
 
-    output::printv31_d(ofs, "atom_position(cartesian)", tau, na);
+    output::printv31_d(ofs, "atom_position(cartesian)", tau.data(), na);
     /*
     for (int i = 0;i < na;i++)
     {
@@ -131,34 +96,26 @@ void Atom::bcast_atom(void)
     Parallel_Common::bcast_bool(coulomb_potential);
     if (GlobalV::MY_RANK != 0)
     {
-        delete[] l_nchi;
-        l_nchi = new int[nwl + 1];
+        this->l_nchi.resize(nwl + 1, 0);
     }
-    Parallel_Common::bcast_int(l_nchi, nwl + 1);
+    Parallel_Common::bcast_int(l_nchi.data(), nwl + 1);
     Parallel_Common::bcast_bool(this->flag_empty_element);
     Parallel_Common::bcast_double(mass);
 
     if (GlobalV::MY_RANK != 0)
     {
         assert(na != 0);
-        delete[] tau;
-        delete[] dis;
-        delete[] taud;
-        delete[] vel;
-        delete[] mag;
-        delete[] angle1;
-        delete[] angle2;
-        delete[] m_loc_;
-        delete[] mbl;
-        tau = new ModuleBase::Vector3<double>[na];
-        dis = new ModuleBase::Vector3<double>[na];
-        taud = new ModuleBase::Vector3<double>[na];
-        vel = new ModuleBase::Vector3<double>[na];
-        mag = new double[na];
-        angle1 = new double[na];
-        angle2 = new double[na];
-        m_loc_ = new ModuleBase::Vector3<double>[na];
-        mbl = new ModuleBase::Vector3<int>[na];
+        this->tau.resize(na, ModuleBase::Vector3<double>(0, 0, 0));
+        this->dis.resize(na, ModuleBase::Vector3<double>(0, 0, 0));
+        this->taud.resize(na, ModuleBase::Vector3<double>(0, 0, 0));
+        this->vel.resize(na, ModuleBase::Vector3<double>(0, 0, 0));
+        this->mag.resize(na, 0);
+        this->angle1.resize(na, 0);
+        this->angle2.resize(na, 0);
+        this->m_loc_.resize(na, ModuleBase::Vector3<double>(0, 0, 0));
+        this->mbl.resize(na, ModuleBase::Vector3<int>(0, 0, 0));
+        this->lambda.resize(na, ModuleBase::Vector3<double>(0, 0, 0));
+        this->constrain.resize(na, ModuleBase::Vector3<int>(0, 0, 0));
     }
 
     for (int i = 0; i < na; i++)
@@ -184,6 +141,12 @@ void Atom::bcast_atom(void)
         Parallel_Common::bcast_int(mbl[i].x);
         Parallel_Common::bcast_int(mbl[i].y);
         Parallel_Common::bcast_int(mbl[i].z);
+        Parallel_Common::bcast_double(lambda[i].x);
+        Parallel_Common::bcast_double(lambda[i].y);
+        Parallel_Common::bcast_double(lambda[i].z);
+        Parallel_Common::bcast_int(constrain[i].x);
+        Parallel_Common::bcast_int(constrain[i].y);
+        Parallel_Common::bcast_int(constrain[i].z);
     }
 
     return;

source/module_cell/atom_spec.h

@@ -13,14 +13,14 @@ class Atom
 
     Atom_pseudo ncpp;
     double mass; // the mass of atom
-    ModuleBase::Vector3<int> *mbl; // whether the atoms can move or not
+    std::vector<ModuleBase::Vector3<int>> mbl; // whether the atoms can move or not
 	bool flag_empty_element = false; // whether is the empty element for bsse.	Peize Lin add 2021.04.07
 
-    int *iw2m; // use iw to find m
-    int *iw2n; // use iw to find n
-    int *iw2l; // use iw to find L
-	int *iw2_ylm;
-	bool *iw2_new;
+    std::vector<int> iw2m; // use iw to find m
+    std::vector<int> iw2n; // use iw to find n
+    std::vector<int> iw2l; // use iw to find L
+	std::vector<int> iw2_ylm;
+	std::vector<bool> iw2_new;
     int nw; // number of local orbitals (l,n,m) of this type
 
     void set_index(void);
@@ -30,21 +30,23 @@ class Atom
 
     int nwl; // max L(Angular momentum) (for local basis)
     double Rcut; //pengfei Li 16-2-29
-    int *l_nchi; // number of chi for each L
+    std::vector<int> l_nchi; // number of chi for each L
     int stapos_wf; // start position of wave functions
 
     std::string label; // atomic symbol
-    ModuleBase::Vector3<double> *tau;// Cartesian coordinates of each atom in this type.
-    ModuleBase::Vector3<double> *dis;// direct displacements of each atom in this type in current step  liuyu modift 2023-03-22
-	ModuleBase::Vector3<double> *taud;// Direct coordinates of each atom in this type.
-    ModuleBase::Vector3<double> *vel;// velocities of each atom in this type.
-    ModuleBase::Vector3<double> *force; // force acting on each atom in this type.
+    std::vector<ModuleBase::Vector3<double>> tau;// Cartesian coordinates of each atom in this type.
+    std::vector<ModuleBase::Vector3<double>> dis;// direct displacements of each atom in this type in current step  liuyu modift 2023-03-22
+	std::vector<ModuleBase::Vector3<double>> taud;// Direct coordinates of each atom in this type.
+    std::vector<ModuleBase::Vector3<double>> vel;// velocities of each atom in this type.
+    std::vector<ModuleBase::Vector3<double>> force; // force acting on each atom in this type.
+    std::vector<ModuleBase::Vector3<double>> lambda; // Lagrange multiplier for each atom in this type. used in deltaspin
+    std::vector<ModuleBase::Vector3<int>> constrain; // constrain for each atom in this type. used in deltaspin
     std::string label_orb; // atomic Element symbol in the orbital file of lcao
 
-	double* mag;
-	double* angle1;//spin angle, added by zhengdy-soc
-	double* angle2;
-    ModuleBase::Vector3<double> *m_loc_;
+	std::vector<double> mag;
+	std::vector<double> angle1;//spin angle, added by zhengdy-soc
+	std::vector<double> angle2;
+    std::vector<ModuleBase::Vector3<double>> m_loc_;
     // Coulomb potential v(r) = z/r
     // It is a local potentail, and has no non-local potential parts.
     bool coulomb_potential = false;

source/module_cell/module_neighbor/test/prepare_unitcell.h

@@ -138,34 +138,25 @@ class UcellTestPrepare
 			ucell->atoms[it].label = this->elements[it];
 			ucell->atoms[it].nw = 0;
 			ucell->atoms[it].nwl = 2;
-			delete[] ucell->atoms[it].l_nchi;
-			ucell->atoms[it].l_nchi = new int[ ucell->atoms[it].nwl+1];
+			ucell->atoms[it].l_nchi.resize(ucell->atoms[it].nwl+1);
 			for(int L=0; L<ucell->atoms[it].nwl+1; L++)
 			{
 				ucell->atoms[it].l_nchi[L] = 1;
 				ucell->atoms[it].nw += (2*L + 1) * ucell->atoms[it].l_nchi[L];
 			}
 			ucell->atoms[it].na = this->natom[it];
 			//coordinates and related physical quantities
-			delete[] ucell->atoms[it].tau;
-			delete[] ucell->atoms[it].dis;
-			delete[] ucell->atoms[it].taud;
-			delete[] ucell->atoms[it].vel;
-			delete[] ucell->atoms[it].mag;
-			delete[] ucell->atoms[it].angle1;
-			delete[] ucell->atoms[it].angle2;
-			delete[] ucell->atoms[it].m_loc_;
-			delete[] ucell->atoms[it].mbl;
-			ucell->atoms[it].tau = new ModuleBase::Vector3<double>[ucell->atoms[it].na];
-			ucell->atoms[it].dis = new ModuleBase::Vector3<double>[ucell->atoms[it].na];
-			ucell->atoms[it].taud = new ModuleBase::Vector3<double>[ucell->atoms[it].na];
-			ucell->atoms[it].vel = new ModuleBase::Vector3<double>[ucell->atoms[it].na];
-			ucell->atoms[it].mag = new double[ucell->atoms[it].na];
-			ucell->atoms[it].angle1 = new double[ucell->atoms[it].na];
-			ucell->atoms[it].angle2 = new double[ucell->atoms[it].na];
-			ucell->atoms[it].m_loc_ = new ModuleBase::Vector3<double>[ucell->atoms[it].na];
-			ucell->atoms[it].mbl = new ModuleBase::Vector3<int>[ucell->atoms[it].na];
+			ucell->atoms[it].tau.resize(ucell->atoms[it].na);
+			ucell->atoms[it].dis.resize(ucell->atoms[it].na);
+			ucell->atoms[it].taud.resize(ucell->atoms[it].na);
+			ucell->atoms[it].vel.resize(ucell->atoms[it].na);
+			ucell->atoms[it].mag.resize(ucell->atoms[it].na);
+			ucell->atoms[it].angle1.resize(ucell->atoms[it].na);
+			ucell->atoms[it].angle2.resize(ucell->atoms[it].na);
+			ucell->atoms[it].m_loc_.resize(ucell->atoms[it].na);
+			ucell->atoms[it].mbl.resize(ucell->atoms[it].na);
 			ucell->atoms[it].mass = ucell->atom_mass[it]; // mass set here
+
 			for(int ia=0; ia<ucell->atoms[it].na; ++ia)
 			{
 				if (ucell->Coordinate == "Direct")

source/module_cell/module_symmetry/test/symmetry_test.cpp

@@ -25,8 +25,8 @@ void SymmetryTest::construct_ucell(stru_ &stru)
     {
         ucell.atoms[i].label = coord[i].atomname;
         ucell.atoms[i].na = coord[i].coordinate.size();
-        ucell.atoms[i].tau = new ModuleBase::Vector3<double>[ucell.atoms[i].na];
-        ucell.atoms[i].taud = new ModuleBase::Vector3<double>[ucell.atoms[i].na];
+        ucell.atoms[i].tau.resize(ucell.atoms[i].na);
+        ucell.atoms[i].taud.resize(ucell.atoms[i].na);
         for (int j = 0; j < ucell.atoms[i].na; j++)
         {
             std::vector<double> this_atom = coord[i].coordinate[j];
@@ -76,10 +76,5 @@ void SymmetryTest::construct_ucell(stru_ &stru)
 
 void SymmetryTest::ClearUcell()
 {
-    for (int i = 0; i < ucell.ntype; i++)
-    {
-        delete[] ucell.atoms[i].tau;
-        delete[] ucell.atoms[i].taud;
-    }
     delete[] ucell.atoms;
 }