Skip to content

Commit 3c50fbf

Browse files
dyzhengdyzheng
authored
Merge pull request #100 from PeizeLin/develop
1. add functions
2 parents f6ac787 + 7506aba commit 3c50fbf

File tree

14 files changed

+543
-30
lines changed

14 files changed

+543
-30
lines changed

‎source/Makefile.Objects‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -245,6 +245,8 @@ read_rho.o\
245245
read_atoms.o\
246246
read_cell_pseudopots.o\
247247
read_dm.o\
248+
read_txt_tools.o\
249+
read_txt_stru.o\
248250
write_pot.o\
249251
write_rho.o\
250252
write_rho_cube.o\

‎source/src_io/CMakeLists.txt‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -23,6 +23,8 @@ add_library(
2323
print_info.cpp
2424
read_dm.cpp
2525
read_rho.cpp
26+
read_txt_tools.cpp
27+
read_txt_stru.cpp
2628
restart.cpp
2729
rwstream.cpp
2830
to_wannier90.cpp

‎source/src_io/numerical_basis.cpp‎

Lines changed: 38 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -95,18 +95,18 @@ void Numerical_Basis::output_overlap( const ModuleBase::ComplexMatrix *psi)
9595
GlobalV::ofs_running << " --------------------------------------------------------" << std::endl;
9696

9797
// search for all k-points.
98-
overlap_Q[ik] = this->cal_overlap_Q(ik, npw, psi[ik], derivative_order);
98+
overlap_Q[ik] = this->cal_overlap_Q(ik, npw, psi[ik], static_cast<double>(derivative_order));
9999
ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running,"cal_overlap_Q");
100100

101101
// (2) generate Sq matrix if necessary.
102102
if (winput::out_spillage == 2)
103103
{
104-
overlap_Sq[ik] = this->cal_overlap_Sq( ik, npw, derivative_order );
104+
overlap_Sq[ik] = this->cal_overlap_Sq( ik, npw, static_cast<double>(derivative_order) );
105105
ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running,"cal_overlap_Sq");
106106
}
107107
}
108108

109-
const ModuleBase::matrix overlap_V = this->cal_overlap_V(psi, derivative_order); // Peize Lin add 2020.04.23
109+
const ModuleBase::matrix overlap_V = this->cal_overlap_V(psi, static_cast<double>(derivative_order)); // Peize Lin add 2020.04.23
110110

111111
#ifdef __MPI
112112
for (int ik=0; ik<GlobalC::kv.nks; ik++)
@@ -139,7 +139,7 @@ ModuleBase::ComplexArray Numerical_Basis::cal_overlap_Q(
139139
const int &ik,
140140
const int &np,
141141
const ModuleBase::ComplexMatrix &psi,
142-
const int derivative_order) const
142+
const double derivative_order) const
143143
{
144144
ModuleBase::TITLE("Numerical_Basis","cal_overlap_Q");
145145
ModuleBase::timer::tick("Numerical_Basis","cal_overlap_Q");
@@ -156,6 +156,8 @@ ModuleBase::ComplexArray Numerical_Basis::cal_overlap_Q(
156156
for (int ig=0; ig<np; ig++)
157157
gk[ig] = GlobalC::wf.get_1qvec_cartesian(ik, ig);
158158

159+
const std::vector<double> gpow = Numerical_Basis::cal_gpow(gk, derivative_order);
160+
159161
const ModuleBase::realArray flq = this->cal_flq(ik, gk);
160162

161163
const ModuleBase::matrix ylm = Numerical_Basis::cal_ylm(gk);
@@ -195,7 +197,7 @@ ModuleBase::ComplexArray Numerical_Basis::cal_overlap_Q(
195197
for (int ig=0; ig<np; ig++)
196198
{
197199
// const std::complex<double> local_tmp = lphase * sk[ig] * ylm(lm, ig) * flq[ig];
198-
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
200+
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
199201
overlap_tmp += conj( local_tmp ) * psi(ib, ig); // psi is bloch orbitals
200202
}
201203
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(
214216
ModuleBase::ComplexArray Numerical_Basis::cal_overlap_Sq(
215217
const int &ik,
216218
const int &np,
217-
const int derivative_order) const
219+
const double derivative_order) const
218220
{
219221
ModuleBase::TITLE("Numerical_Basis","cal_overlap_Sq");
220222
ModuleBase::timer::tick("Numerical_Basis","cal_overlap_Sq");
@@ -232,6 +234,8 @@ ModuleBase::ComplexArray Numerical_Basis::cal_overlap_Sq(
232234
for (int ig=0; ig<np; ig++)
233235
gk[ig] = GlobalC::wf.get_1qvec_cartesian(ik, ig);
234236

237+
const std::vector<double> gpow = Numerical_Basis::cal_gpow(gk, derivative_order);
238+
235239
const ModuleBase::realArray flq = this->cal_flq(ik, gk);
236240

237241
const ModuleBase::matrix ylm = Numerical_Basis::cal_ylm(gk);
@@ -281,7 +285,7 @@ ModuleBase::ComplexArray Numerical_Basis::cal_overlap_Sq(
281285

282286
std::vector<std::complex<double>> about_ig1(np, std::complex<double>(0.0,0.0));
283287
for (int ig=0; ig<np; ig++)
284-
about_ig1[ig] = conj( lphase1 * sk1[ig] * ylm(lm1, ig) ) * pow(gk[ig].norm2(),derivative_order); // Peize Lin add for dpsi 2020.04.23
288+
about_ig1[ig] = conj( lphase1 * sk1[ig] * ylm(lm1, ig) ) * gpow[ig]; // Peize Lin add for dpsi 2020.04.23
285289

286290
for (int m2=0; m2<2*l2+1; m2++) // 2.6
287291
{
@@ -330,18 +334,20 @@ ModuleBase::ComplexArray Numerical_Basis::cal_overlap_Sq(
330334
// Peize Lin add for dpsi 2020.04.23
331335
ModuleBase::matrix Numerical_Basis::cal_overlap_V(
332336
const ModuleBase::ComplexMatrix *psi,
333-
const int derivative_order)
337+
const double derivative_order)
334338
{
335339
ModuleBase::matrix overlap_V(GlobalC::kv.nks, GlobalV::NBANDS);
336340
for(int ik=0; ik<GlobalC::kv.nks; ++ik)
337341
{
342+
std::vector<ModuleBase::Vector3<double>> gk(GlobalC::kv.ngk[ik]);
343+
for (int ig=0; ig<gk.size(); ig++)
344+
gk[ig] = GlobalC::wf.get_1qvec_cartesian(ik, ig);
345+
346+
const std::vector<double> gpow = Numerical_Basis::cal_gpow(gk, derivative_order);
347+
338348
for(int ib=0; ib<GlobalV::NBANDS; ++ib)
339-
{
340349
for(int ig=0; ig<GlobalC::kv.ngk[ik]; ++ig)
341-
{
342-
overlap_V(ik,ib)+= norm(psi[ik](ib,ig)) * pow(GlobalC::wf.get_1qvec_cartesian(ik,ig).norm2(),derivative_order);
343-
}
344-
}
350+
overlap_V(ik,ib)+= norm(psi[ik](ib,ig)) * gpow[ig];
345351
}
346352
return overlap_V;
347353
}
@@ -368,6 +374,25 @@ ModuleBase::matrix Numerical_Basis::cal_ylm(const std::vector<ModuleBase::Vector
368374
return ylm;
369375
}
370376

377+
std::vector<double> Numerical_Basis::cal_gpow (const std::vector<ModuleBase::Vector3<double>> &gk, const double derivative_order)
378+
{
379+
constexpr double thr = 1E-12;
380+
std::vector<double> gpow(gk.size(), 0.0);
381+
for (int ig=0; ig<gpow.size(); ++ig)
382+
{
383+
if (derivative_order>=0)
384+
{
385+
gpow[ig] = std::pow(gk[ig].norm2(),derivative_order);
386+
}
387+
else
388+
{
389+
if (gk[ig].norm2() >= thr)
390+
gpow[ig] = std::pow(gk[ig].norm2(),derivative_order);
391+
}
392+
}
393+
return gpow;
394+
}
395+
371396
std::vector<ModuleBase::IntArray> Numerical_Basis::init_mu_index(void)
372397
{
373398
GlobalV::ofs_running << " Initialize the mu index" << std::endl;

‎source/src_io/numerical_basis.h‎

Lines changed: 5 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -34,19 +34,21 @@ class Numerical_Basis
3434

3535
ModuleBase::ComplexArray cal_overlap_Q(
3636
const int &ik, const int &np, const ModuleBase::ComplexMatrix &psi,
37-
const int derivative_order ) const;
37+
const double derivative_order ) const;
3838

3939
ModuleBase::ComplexArray cal_overlap_Sq(
4040
const int &ik,
4141
const int &np,
42-
const int derivative_order ) const;
42+
const double derivative_order ) const;
4343

44-
static ModuleBase::matrix cal_overlap_V(const ModuleBase::ComplexMatrix *psi, const int derivative_order);
44+
static ModuleBase::matrix cal_overlap_V(const ModuleBase::ComplexMatrix *psi, const double derivative_order);
4545

4646
ModuleBase::realArray cal_flq(const int ik, const std::vector<ModuleBase::Vector3<double>> &gk) const;
4747

4848
static ModuleBase::matrix cal_ylm(const std::vector<ModuleBase::Vector3<double>> &gk);
4949

50+
static std::vector<double> cal_gpow (const std::vector<ModuleBase::Vector3<double>> &gk, const double derivative_order);
51+
5052
static void output_info(
5153
std::ofstream &ofs,
5254
const Bessel_Basis &bessel_basis);

‎source/src_io/read_txt_stru.cpp‎

Lines changed: 127 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,127 @@
1+
//=======================
2+
// AUTHOR : Peize Lin
3+
// DATE : 2021-12-08
4+
//=======================
5+
6+
#include "read_txt_stru.h"
7+
#include "read_txt_tools.h"
8+
9+
#include <algorithm>
10+
11+
namespace Read_Txt_Stru
12+
{
13+
//lat_info:
14+
// {
15+
// Lattice_Constant: [30.0, Bohr],
16+
// Lattice_Vectors:
17+
// [1.0, 0.0, 0.0,
18+
// 0.0, 1.0, 0.0,
19+
// 0.0, 0.0, 1.0],
20+
// Atomic_Positions: [Cartesian, Angstrom]
21+
// }
22+
//elements_info:
23+
// {
24+
// O: {
25+
// Pseudo_Potantial: [O_ONCV_PBE-1.0.upf],
26+
// Numerical_Orbital: [orb_O.dat]
27+
// },
28+
// C: {
29+
// Pseudo_Potantial: [C_ONCV_PBE-1.0.upf]
30+
// }
31+
// }
32+
//atoms:
33+
// [
34+
// [O, 0, 0, 0],
35+
// [H, 0, 0, 3, Force, 1, 1, 1]
36+
// ]
37+
std::tuple<
38+
std::map<std::string, std::vector<std::string>>,
39+
std::map<std::string, std::map<std::string, std::vector<std::string>>>,
40+
std::vector<std::vector<std::string>>>
41+
read_stru(const std::string &file_name)
42+
{
43+
const std::set<std::string> labels_lat = {"Lattice_Constant", "Lattice_Vectors", "Atomic_Positions"};
44+
std::set<std::string> labels = {"Element", "Atoms"};
45+
labels.insert(labels_lat.begin(), labels_lat.end());
46+
47+
std::vector<std::vector<std::string>> stru = Read_Txt_Tools::read_file_to_vector(file_name, {"#","\\"});
48+
49+
std::map<std::string, std::vector<std::string>> lat_info;
50+
std::map<std::string, std::map<std::string, std::vector<std::string>>> elements_info;
51+
std::vector<std::vector<std::string>> atoms_info;
52+
53+
while(!stru.empty())
54+
{
55+
std::vector<std::vector<std::string>> stru_one = Read_Txt_Tools::cut_paragraph(stru, labels);
56+
if(stru_one[0][0]=="Element")
57+
{
58+
std::map<std::string, std::vector<std::string>> &element_info = elements_info[stru_one[0][1]];
59+
for(auto ptr=stru_one.begin()+1; ptr<stru_one.end(); ++ptr)
60+
{
61+
element_info[(*ptr)[0]] = std::vector<std::string>(ptr->begin()+1, ptr->end());
62+
}
63+
}
64+
else if(stru_one[0][0]=="Atoms")
65+
{
66+
// [
67+
// [O, 0, 0, 0],
68+
// [H, 0, 0, 3],
69+
// [Force, 1, 1, 1]
70+
// ]
71+
if(stru_one[0].size()==1)
72+
{
73+
atoms_info.insert(atoms_info.end(), stru_one.begin()+1, stru_one.end());
74+
}
75+
else
76+
{
77+
std::vector<std::string> atom0 = std::vector<std::string>(stru_one[0].begin()+1, stru_one[0].end());
78+
atoms_info.push_back(std::move(atom0));
79+
atoms_info.insert(atoms_info.end(), stru_one.begin()+1, stru_one.end());
80+
}
81+
82+
}
83+
else
84+
{
85+
const auto ptr = std::find(labels_lat.begin(), labels_lat.end(), stru_one[0][0]);
86+
if(ptr!=labels_lat.end())
87+
{
88+
const std::vector<std::string> data = Read_Txt_Tools::chain(stru_one);
89+
lat_info[data[0]] = std::vector<std::string>(data.begin()+1, data.end());
90+
}
91+
else
92+
throw std::invalid_argument(stru_one[0][0]);
93+
}
94+
}
95+
96+
std::set<std::string> elements_label;
97+
for(const auto & i : elements_info)
98+
elements_label.insert(i.first);
99+
std::vector<std::vector<std::string>> atoms_info_new = Read_Txt_Stru::organize_atom_info(std::move(atoms_info), elements_label);
100+
101+
return make_tuple(std::move(lat_info), std::move(elements_info), std::move(atoms_info_new));
102+
}
103+
104+
// atoms_info_old:
105+
// [
106+
// [H, 0, 0, 0],
107+
// [O, 0, 0, 3],
108+
// [Force, 1, 1, 1]
109+
// ]
110+
// atoms_info_new:
111+
// [
112+
// [H, 0, 0, 0],
113+
// [O, 0, 0, 3, Force, 1, 1, 1]
114+
// ]
115+
std::vector<std::vector<std::string>> organize_atom_info(
116+
std::vector<std::vector<std::string>> &&atoms_info_old,
117+
const std::set<std::string> &elements_label)
118+
{
119+
std::vector<std::vector<std::string>> atoms_info_new;
120+
while(!atoms_info_old.empty())
121+
{
122+
const std::vector<std::vector<std::string>> atom_one = Read_Txt_Tools::cut_paragraph(atoms_info_old, elements_label);
123+
atoms_info_new.push_back(Read_Txt_Tools::chain(atom_one));
124+
}
125+
return atoms_info_new;
126+
}
127+
}

‎source/src_io/read_txt_stru.h‎

Lines changed: 65 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,65 @@
1+
//=======================
2+
// AUTHOR : Peize Lin
3+
// DATE : 2021-12-08
4+
//=======================
5+
6+
#ifndef READ_TXT_STRU_H
7+
#define READ_TXT_STRU_H
8+
9+
#include <tuple>
10+
#include <vector>
11+
#include <map>
12+
#include <set>
13+
#include <string>
14+
15+
namespace Read_Txt_Stru
16+
{
17+
//lat_info:
18+
// {
19+
// Lattice_Constant: [30.0, Bohr],
20+
// Lattice_Vectors:
21+
// [1.0, 0.0, 0.0,
22+
// 0.0, 1.0, 0.0,
23+
// 0.0, 0.0, 1.0],
24+
// Atomic_Positions: [Cartesian, Angstrom]
25+
// }
26+
//elements_info:
27+
// {
28+
// O: {
29+
// Pseudo_Potantial: [O_ONCV_PBE-1.0.upf],
30+
// Numerical_Orbital: [orb_O.dat]
31+
// },
32+
// C: {
33+
// Pseudo_Potantial: [C_ONCV_PBE-1.0.upf]
34+
// }
35+
// }
36+
//atoms:
37+
// [
38+
// [O, 0, 0, 0],
39+
// [H, 0, 0, 3, Force, 1, 1, 1]
40+
// ]
41+
std::tuple<
42+
std::map<std::string, std::vector<std::string>>,
43+
std::map<std::string, std::map<std::string, std::vector<std::string>>>,
44+
std::vector<std::vector<std::string>>>
45+
read_stru(const std::string &file_name);
46+
47+
48+
49+
// atoms_info_old:
50+
// [
51+
// [H, 0, 0, 0],
52+
// [O, 0, 0, 3],
53+
// [Force, 1, 1, 1]
54+
// ]
55+
// atoms_info_new:
56+
// [
57+
// [H, 0, 0, 0],
58+
// [O, 0, 0, 3, Force, 1, 1, 1]
59+
// ]
60+
std::vector<std::vector<std::string>> organize_atom_info(
61+
std::vector<std::vector<std::string>> &&atoms_info_old,
62+
const std::set<std::string> &elements_label);
63+
}
64+
65+
#endif

0 commit comments

Comments
 (0)