Merge pull request #608 from deepmodeling/ABACUS_2.2.0_beta

merge linpz new code to develop branch

Author: dyzheng

source/Makefile.Objects

@@ -249,6 +249,16 @@ read_rho.o\
 read_atoms.o\
 read_cell_pseudopots.o\
 read_dm.o\
+read_txt_tools.o\
+read_txt_stru.o\
+read_txt_input_value.o\
+read_txt_input_item.o\
+read_txt_input_list.o\
+read_txt_input_process.o\
+read_txt_input_process_global.o\
+read_txt_input-general.o\
+read_txt_input-pw.o\
+read_txt_input-spectrum.o\
 write_pot.o\
 write_rho.o\
 write_rho_cube.o\

source/module_base/timer.cpp

@@ -8,7 +8,11 @@
 #include<vector>
 
 #ifdef __MPI
-#include "mpi.h"
+#include <mpi.h>
+#endif
+
+#ifdef _OPENMP
+#include <omp.h>
 #endif
 
 namespace ModuleBase
@@ -67,7 +71,11 @@ void timer::tick(const std::string &class_name,const std::string &name)
 	if (disabled)
 		return;
 
-	Timer_One &timer_one = timer_pool[class_name][name];
+#ifdef _OPENMP
+	if(!omp_get_thread_num())
+#endif
+	{
+		Timer_One &timer_one = timer_pool[class_name][name];
 
 //----------------------------------------------------------
 // CALL MEMBER FUNCTION :
@@ -79,32 +87,33 @@ void timer::tick(const std::string &class_name,const std::string &name)
 // if start_flag == false, means it's the end of this counting,
 // so we add the time during this two 'time point'  to the clock time storage.
 //----------------------------------------------------------
-	if(timer_one.start_flag)
-	{
+		if(timer_one.start_flag)
+		{
 #ifdef __MPI
-		timer_one.cpu_start = MPI_Wtime();
+			timer_one.cpu_start = MPI_Wtime();
 #else
-		timer_one.cpu_start = cpu_time();
+			timer_one.cpu_start = cpu_time();
 #endif
-		++timer_one.calls;
-		timer_one.start_flag = false;
-	}
-	else
-	{
+			++timer_one.calls;
+			timer_one.start_flag = false;
+		}
+		else
+		{
 #ifdef __MPI
-		timer_one.cpu_second += MPI_Wtime() - timer_one.cpu_start;
+			timer_one.cpu_second += MPI_Wtime() - timer_one.cpu_start;
 #else
-		// if(class_name=="electrons"&&name=="c_bands")
-		// {
-		// 	cout<<"call times"<<timer_one.calls<<endl;
-		// 	cout<<"electrons c_bands cost time:"<<endl;
-		// 	cout<<cpu_time()<<"-"<<timer_one.cpu_start<<endl;
-		// }
-
-		timer_one.cpu_second += (cpu_time() - timer_one.cpu_start);
+			// if(class_name=="electrons"&&name=="c_bands")
+			// {
+			// 	cout<<"call times"<<timer_one.calls<<endl;
+			// 	cout<<"electrons c_bands cost time:"<<endl;
+			// 	cout<<cpu_time()<<"-"<<timer_one.cpu_start<<endl;
+			// }
+
+			timer_one.cpu_second += (cpu_time() - timer_one.cpu_start);
 #endif
-		timer_one.start_flag = true;
-	}
+			timer_one.start_flag = true;
+		}
+	} // end if(!omp_get_thread_num())
 }
 
 long double timer::print_until_now(void)

source/src_external/src_test/print_tmp.h

@@ -1,4 +1,5 @@
-#pragma once
+#ifndef PRINT_TMP_H
+#define PRINT_TMP_H
 
 #include<iostream>
 #include<fstream>
@@ -42,4 +43,6 @@ void print_matrixess(const std::map<size_t,std::map<size_t,std::map<Abfs::Vector
 				ofs<<c1.first<<"\t"<<c2.first<<"\t"<<c3.first<<std::endl;
 				c3.second->print(ofs, 1E-10)<<std::endl;
 			}
-}
+}
+
+#endif

source/src_external/src_test/src_global/element_basis_index-test.h

@@ -3,8 +3,8 @@
 // DATE : 2016-06-02
 //==========================================================
 
-#pragma once
-
+#ifndef ELEMENT_BASIS_INDEX_TEST_H
+#define ELEMENT_BASIS_INDEX_TEST_H
 
 #include <iostream>
 using namespace std;
@@ -41,4 +41,6 @@ static std::ostream & operator << (std::ostream & os, const ModuleBase::Element_
 		}
 	}
 	return os;
-}
+}
+
+#endif

source/src_io/CMakeLists.txt

@@ -23,6 +23,16 @@ add_library(
     print_info.cpp
     read_dm.cpp
     read_rho.cpp
+    read_txt_tools.cpp
+    read_txt_input_value.cpp
+    read_txt_input_item.cpp
+	read_txt_input_list.cpp
+	read_txt_input_process.cpp
+	read_txt_input_process_global.cpp
+	read_txt_input-general.cpp
+	read_txt_input-pw.cpp
+	read_txt_input-spectrum.cpp
+    read_txt_stru.cpp
     restart.cpp
     rwstream.cpp
     to_wannier90.cpp

source/src_io/dos.cpp

@@ -347,8 +347,16 @@ void Dos::nscf_band(
 	}
 	MPI_Barrier(MPI_COMM_WORLD);
 
+	std::vector<double> klength;
+	klength.resize(nks);
+	klength[0] = 0.0;
 	for(int ik=0; ik<nks; ik++)
 	{
+		if (ik>0)
+		{
+			auto delta=GlobalC::kv.kvec_c[ik]-GlobalC::kv.kvec_c[ik-1];
+			klength[ik] = klength[ik-1] + delta.norm();
+		}
 		if ( GlobalV::MY_POOL == GlobalC::Pkpoints.whichpool[ik] )
 		{
 			const int ik_now = ik - GlobalC::Pkpoints.startk_pool[GlobalV::MY_POOL];
@@ -360,6 +368,7 @@ void Dos::nscf_band(
 					ofs << std::setprecision(8);
 					//start from 1
 					ofs << ik+1;
+					ofs << " " << klength[ik] << " ";
 					for(int ib = 0; ib < nband; ib++)
 					{
 						ofs << " " << (ekb[ik_now+is*nks][ib]-fermie) * ModuleBase::Ry_to_eV;

source/src_io/numerical_basis.cpp

@@ -95,18 +95,18 @@ void Numerical_Basis::output_overlap( const ModuleBase::ComplexMatrix *psi)
             GlobalV::ofs_running << " --------------------------------------------------------" << std::endl;
 
             // search for all k-points.
-            overlap_Q[ik] = this->cal_overlap_Q(ik, npw, psi[ik], derivative_order);
+            overlap_Q[ik] = this->cal_overlap_Q(ik, npw, psi[ik], static_cast<double>(derivative_order));
             ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running,"cal_overlap_Q");
 
             // (2) generate Sq matrix if necessary.
             if (winput::out_spillage == 2)
             {
-                overlap_Sq[ik] = this->cal_overlap_Sq( ik, npw, derivative_order );
+                overlap_Sq[ik] = this->cal_overlap_Sq( ik, npw, static_cast<double>(derivative_order) );
                 ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running,"cal_overlap_Sq");
             }
         }
 
-        const ModuleBase::matrix overlap_V = this->cal_overlap_V(psi, derivative_order);		// Peize Lin add 2020.04.23
+        const ModuleBase::matrix overlap_V = this->cal_overlap_V(psi, static_cast<double>(derivative_order));		// Peize Lin add 2020.04.23
 
     #ifdef __MPI
         for (int ik=0; ik<GlobalC::kv.nks; ik++)
@@ -139,7 +139,7 @@ ModuleBase::ComplexArray Numerical_Basis::cal_overlap_Q(
     const int &ik,
     const int &np,
     const ModuleBase::ComplexMatrix &psi,
-	const int derivative_order) const
+	const double derivative_order) const
 {
     ModuleBase::TITLE("Numerical_Basis","cal_overlap_Q");
     ModuleBase::timer::tick("Numerical_Basis","cal_overlap_Q");
@@ -156,6 +156,8 @@ ModuleBase::ComplexArray Numerical_Basis::cal_overlap_Q(
     for (int ig=0; ig<np; ig++)
         gk[ig] = GlobalC::wf.get_1qvec_cartesian(ik, ig);
 
+    const std::vector<double> gpow = Numerical_Basis::cal_gpow(gk, derivative_order);
+
 	const ModuleBase::realArray flq = this->cal_flq(ik, gk);
 
     const ModuleBase::matrix ylm = Numerical_Basis::cal_ylm(gk);
@@ -195,7 +197,7 @@ ModuleBase::ComplexArray Numerical_Basis::cal_overlap_Q(
                             for (int ig=0; ig<np; ig++)
                             {
 //                              const std::complex<double> local_tmp = lphase * sk[ig] * ylm(lm, ig) * flq[ig];
-                                const std::complex<double> local_tmp = lphase * sk[ig] * ylm(lm, ig) * flq(L,ie,ig) * pow(gk[ig].norm2(),derivative_order);		// Peize Lin add for dpsi 2020.04.23
+                                const std::complex<double> local_tmp = lphase * sk[ig] * ylm(lm, ig) * flq(L,ie,ig) * gpow[ig];		// Peize Lin add for dpsi 2020.04.23
                                 overlap_tmp += conj( local_tmp ) * psi(ib, ig); // psi is bloch orbitals
                             }
                             overlap_Q(ib, this->mu_index[T](I, L, N, m), ie) = overlap_tmp;
@@ -214,7 +216,7 @@ ModuleBase::ComplexArray Numerical_Basis::cal_overlap_Q(
 ModuleBase::ComplexArray Numerical_Basis::cal_overlap_Sq(
     const int &ik,
     const int &np,
-	const int derivative_order) const
+	const double derivative_order) const
 {
     ModuleBase::TITLE("Numerical_Basis","cal_overlap_Sq");
     ModuleBase::timer::tick("Numerical_Basis","cal_overlap_Sq");
@@ -232,6 +234,8 @@ ModuleBase::ComplexArray Numerical_Basis::cal_overlap_Sq(
     for (int ig=0; ig<np; ig++)
         gk[ig] = GlobalC::wf.get_1qvec_cartesian(ik, ig);
 
+    const std::vector<double> gpow = Numerical_Basis::cal_gpow(gk, derivative_order);
+
 	const ModuleBase::realArray flq = this->cal_flq(ik, gk);
 
     const ModuleBase::matrix ylm = Numerical_Basis::cal_ylm(gk);
@@ -281,7 +285,7 @@ ModuleBase::ComplexArray Numerical_Basis::cal_overlap_Sq(
 
                                         std::vector<std::complex<double>> about_ig1(np, std::complex<double>(0.0,0.0));
                                         for (int ig=0; ig<np; ig++)
-                                            about_ig1[ig] = conj( lphase1 * sk1[ig] * ylm(lm1, ig) ) * pow(gk[ig].norm2(),derivative_order);		// Peize Lin add for dpsi 2020.04.23
+                                            about_ig1[ig] = conj( lphase1 * sk1[ig] * ylm(lm1, ig) ) * gpow[ig];		// Peize Lin add for dpsi 2020.04.23
 
                                         for (int m2=0; m2<2*l2+1; m2++) // 2.6
                                         {
@@ -330,18 +334,20 @@ ModuleBase::ComplexArray Numerical_Basis::cal_overlap_Sq(
 // Peize Lin add for dpsi 2020.04.23
 ModuleBase::matrix Numerical_Basis::cal_overlap_V(
 	const ModuleBase::ComplexMatrix *psi,
-	const int derivative_order)
+	const double derivative_order)
 {
 	ModuleBase::matrix overlap_V(GlobalC::kv.nks, GlobalV::NBANDS);
 	for(int ik=0; ik<GlobalC::kv.nks; ++ik)
 	{
+        std::vector<ModuleBase::Vector3<double>> gk(GlobalC::kv.ngk[ik]);
+        for (int ig=0; ig<gk.size(); ig++)
+            gk[ig] = GlobalC::wf.get_1qvec_cartesian(ik, ig);
+
+        const std::vector<double> gpow = Numerical_Basis::cal_gpow(gk, derivative_order);
+
 		for(int ib=0; ib<GlobalV::NBANDS; ++ib)
-		{
 			for(int ig=0; ig<GlobalC::kv.ngk[ik]; ++ig)
-			{
-				overlap_V(ik,ib)+= norm(psi[ik](ib,ig)) * pow(GlobalC::wf.get_1qvec_cartesian(ik,ig).norm2(),derivative_order);
-			}
-		}
+				overlap_V(ik,ib)+= norm(psi[ik](ib,ig)) * gpow[ig];
 	}
 	return overlap_V;
 }
@@ -368,6 +374,25 @@ ModuleBase::matrix Numerical_Basis::cal_ylm(const std::vector<ModuleBase::Vector
     return ylm;
 }
 
+std::vector<double> Numerical_Basis::cal_gpow (const std::vector<ModuleBase::Vector3<double>> &gk, const double derivative_order)
+{
+    constexpr double thr = 1E-12;
+    std::vector<double> gpow(gk.size(), 0.0);
+    for (int ig=0; ig<gpow.size(); ++ig)
+    {
+        if (derivative_order>=0)
+        {
+            gpow[ig] = std::pow(gk[ig].norm2(),derivative_order);
+        }
+        else
+        {
+            if (gk[ig].norm2() >= thr)
+                gpow[ig] = std::pow(gk[ig].norm2(),derivative_order);
+        }
+    }
+    return gpow;
+}
+
 std::vector<ModuleBase::IntArray> Numerical_Basis::init_mu_index(void)
 {
 	GlobalV::ofs_running << " Initialize the mu index" << std::endl;

source/src_io/numerical_basis.h

@@ -34,19 +34,21 @@ class Numerical_Basis
 
 	ModuleBase::ComplexArray cal_overlap_Q(
 		const int &ik, const int &np, const ModuleBase::ComplexMatrix &psi,
-		const int derivative_order ) const;
+		const double derivative_order ) const;
 
 	ModuleBase::ComplexArray cal_overlap_Sq(
 		const int &ik, 
 		const int &np,
-		const int derivative_order ) const;
+		const double derivative_order ) const;
 
-	static ModuleBase::matrix cal_overlap_V(const ModuleBase::ComplexMatrix *psi, const int derivative_order);
+	static ModuleBase::matrix cal_overlap_V(const ModuleBase::ComplexMatrix *psi, const double derivative_order);
 
 	ModuleBase::realArray cal_flq(const int ik, const std::vector<ModuleBase::Vector3<double>> &gk) const;
 
 	static ModuleBase::matrix cal_ylm(const std::vector<ModuleBase::Vector3<double>> &gk);
 
+	static std::vector<double> cal_gpow (const std::vector<ModuleBase::Vector3<double>> &gk, const double derivative_order);
+
 	static void output_info(
 		std::ofstream &ofs,
 		const Bessel_Basis &bessel_basis);

source/src_io/read_txt_input-general.cpp

@@ -0,0 +1,64 @@
+//=======================
+// AUTHOR : Peize Lin
+// DATE :   2021-12-13
+//=======================
+
+#include "read_txt_input_list.h"
+
+#include "src_io/read_txt_tools.h"
+#include <stdexcept>
+
+// for convert
+#include "module_base/global_variable.h"
+#include "src_pw/global.h"
+
+namespace Read_Txt_Input
+{
+	void Input_List::set_items_general()
+	{
+		this->output_labels.push_back("Parameters (1.General)");
+
+		{
+			Input_Item item("nspin");
+			item.default_1(1);
+			item.annotation = "1: single spin; 2: up and down spin; 4: noncollinear spin";
+			item.check_transform = [](Input_Item &self)
+			{
+				if(!Read_Txt_Tools::in_set(self.values[0].geti(), {1,2,4}))
+					throw std::invalid_argument("nspin must be 1,2,4");
+			};
+			item.convert = [](const Input_Item &self)
+			{
+				GlobalV::NSPIN = self.values[0].geti();
+			};
+			this->add_item(item);
+		}
+
+		{
+			Input_Item item("ntype");
+			item.default_1(0);
+			item.annotation = "atom species number";
+			item.check_transform = [](Input_Item &self)
+			{
+				if(self.values[0].geti()<=0)
+					throw std::invalid_argument("ntype must > 0");
+			};
+			item.convert = [](const Input_Item &self)
+			{
+				GlobalC::ucell.ntype = self.values[0].geti();
+			};
+			this->add_item(item);
+		}
+
+		{
+			Input_Item item("symmetry");
+			item.default_1(false);
+			item.annotation = "turn symmetry on or off";
+			item.convert = [](const Input_Item &self)
+			{
+				ModuleSymmetry::Symmetry::symm_flag = self.values[0].getb();
+			};
+			this->add_item(item);
+		}
+	}
+}