Skip to content

Commit a671dcb

Browse files
Wuming-HUSTWuming-HUST
committed
Feature: add DFT-1/2 and shell DFT-1/2, currently only support PW esolver_ks_pw.
* Added Sep, Sep_Cell, and VSep to organize the self-energy potential of DFT-1/2 * Added a new effective potential pot_sep for calculating the self-energy potential * Added initialization of the self-energy potential in the esolver_ks_pw control flow * Added the keyword SEP_FILES in the STRU file for reading self-energy potential files * Added the dfthalf_type keyword in INPUT to enable DFT-1/2 and shell DFT-1/2
1 parent 9e6bb70 commit a671dcb

29 files changed

+3355
-77
lines changed

‎source/source_cell/CMakeLists.txt‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -26,6 +26,8 @@ add_library(
2626
print_cell.cpp
2727
read_atom_species.cpp
2828
k_vector_utils.cpp
29+
sep.cpp
30+
sep_cell.cpp
2931
)
3032

3133
if(ENABLE_COVERAGE)

‎source/source_cell/sep.cpp‎

Lines changed: 122 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,122 @@
1+
#include "sep.h"
2+
3+
#include "source_base/global_variable.h"
4+
#include "source_base/parallel_common.h"
5+
#include "source_base/tool_title.h"
6+
#include "source_io/output.h"
7+
8+
#include <fstream>
9+
#include <sstream>
10+
#include <string>
11+
12+
SepPot::SepPot()
13+
{
14+
}
15+
16+
SepPot::~SepPot()
17+
{
18+
delete[] r;
19+
r = nullptr;
20+
delete[] rv;
21+
rv = nullptr;
22+
}
23+
24+
int SepPot::read_sep(std::ifstream& ifs)
25+
{
26+
std::string line;
27+
while (std::getline(ifs, line))
28+
{
29+
std::istringstream iss(line);
30+
std::string key;
31+
iss >> key;
32+
33+
if (key == "Sep.Element")
34+
{
35+
iss >> label;
36+
}
37+
else if (key == "Sep.XcType")
38+
{
39+
iss >> xc_type;
40+
}
41+
else if (key == "Sep.Orbital")
42+
{
43+
iss >> orbital;
44+
}
45+
else if (key == "Sep.Points")
46+
{
47+
iss >> mesh;
48+
delete[] r;
49+
r = new double[mesh];
50+
delete[] rv;
51+
rv = new double[mesh];
52+
}
53+
else if (key == "Sep.StripAmount")
54+
{
55+
iss >> strip_elec;
56+
}
57+
else if (key == "<Sep.Potential")
58+
{
59+
double r_val, rv_val;
60+
int idx = 0;
61+
while (std::getline(ifs, line) && line != "Sep.Potential>")
62+
{
63+
std::istringstream data_line(line);
64+
if (data_line >> r_val >> rv_val)
65+
{
66+
r[idx] = r_val;
67+
rv[idx] = rv_val;
68+
idx++;
69+
}
70+
}
71+
break;
72+
}
73+
}
74+
return 0;
75+
}
76+
77+
void SepPot::print_sep_info(std::ofstream& ofs)
78+
{
79+
ofs << "\n sep_vl:";
80+
ofs << "\n sep_info:";
81+
ofs << "\n label " << label;
82+
ofs << "\n xc " << xc_type;
83+
ofs << "\n orbital " << orbital;
84+
ofs << "\n strip electron" << strip_elec;
85+
}
86+
87+
void SepPot::print_sep_vsep(std::ofstream& ofs)
88+
{
89+
ofs << "\n mesh " << mesh;
90+
output::printr1_d(ofs, " r : ", r, mesh);
91+
output::printr1_d(ofs, " vsep : ", rv, mesh);
92+
ofs << "\n -----------------------------";
93+
}
94+
95+
#ifdef __MPI
96+
97+
void SepPot::bcast_sep()
98+
{
99+
ModuleBase::TITLE("SepPot", "bcast_sep");
100+
Parallel_Common::bcast_bool(is_enable);
101+
Parallel_Common::bcast_double(r_in);
102+
Parallel_Common::bcast_double(r_out);
103+
Parallel_Common::bcast_double(r_power);
104+
Parallel_Common::bcast_double(enhence_a);
105+
Parallel_Common::bcast_string(label);
106+
Parallel_Common::bcast_string(xc_type);
107+
Parallel_Common::bcast_string(orbital);
108+
Parallel_Common::bcast_int(strip_elec);
109+
Parallel_Common::bcast_int(mesh);
110+
111+
if (GlobalV::MY_RANK != 0 && mesh > 0)
112+
{
113+
r = new double[mesh];
114+
rv = new double[mesh];
115+
}
116+
117+
Parallel_Common::bcast_double(r, mesh);
118+
Parallel_Common::bcast_double(rv, mesh);
119+
120+
return;
121+
}
122+
#endif // __MPI

‎source/source_cell/sep.h‎

Lines changed: 39 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,39 @@
1+
#ifndef SEP_H
2+
#define SEP_H
3+
4+
#include <fstream>
5+
#include <string>
6+
7+
/**
8+
* Sep Potential for DFT-1/2 etc.
9+
*
10+
* Sep Potential
11+
*/
12+
class SepPot
13+
{
14+
public:
15+
SepPot();
16+
~SepPot();
17+
18+
bool is_enable = false;
19+
double r_in = 0.0; /**< cut-off radius inner */
20+
double r_out = 0.0; /**< cut-off radius outter */
21+
double r_power = 20.0; /**< shell function exp factor */
22+
double enhence_a = 1.0; /**< scale sep potential */
23+
std::string label; /**< element nameof sep */
24+
std::string xc_type; /**< Exch-Corr type */
25+
std::string orbital; /** atomic angular moment s,p,d,f */
26+
int mesh = 0; /**< number of points in radial mesh */
27+
int strip_elec = 0; /**< strip electron amount 1->0.01 50->0.5 */
28+
double* r = nullptr; /**< ridial mesh */
29+
double* rv = nullptr; /**< sep potential, but rV, unit: Ry */
30+
31+
int read_sep(std::ifstream& is);
32+
void print_sep_info(std::ofstream& ofs);
33+
void print_sep_vsep(std::ofstream& ofs);
34+
#ifdef __MPI
35+
void bcast_sep();
36+
#endif /* ifdef __MPI */
37+
};
38+
39+
#endif /* ifndef SEP_H */

‎source/source_cell/sep_cell.cpp‎

Lines changed: 127 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,127 @@
1+
#include "sep_cell.h"
2+
3+
#include "source_base/global_function.h"
4+
#include "source_base/global_variable.h"
5+
#include "source_base/parallel_common.h"
6+
#include "source_base/tool_title.h"
7+
#include "source_cell/unitcell.h"
8+
9+
#include <algorithm>
10+
#include <string>
11+
12+
namespace GlobalC
13+
{
14+
Sep_Cell sep_cell;
15+
}
16+
17+
Sep_Cell::Sep_Cell() noexcept : ntype(0), omega(0.0), tpiba2(0.0)
18+
{
19+
}
20+
21+
Sep_Cell::~Sep_Cell() noexcept = default;
22+
23+
void Sep_Cell::init(const int ntype_in)
24+
{
25+
this->ntype = ntype_in;
26+
this->seps.resize(ntype);
27+
this->sep_enable.resize(ntype);
28+
std::fill(this->sep_enable.begin(), this->sep_enable.end(), false);
29+
}
30+
31+
void Sep_Cell::set_omega(const double omega_in, const double tpiba2_in)
32+
{
33+
this->omega = omega_in;
34+
this->tpiba2 = tpiba2_in;
35+
}
36+
37+
/**
38+
* read sep potential files
39+
*
40+
* need to add following lines in STRU file, and order of elements must match ATOMIC_SPECIES.
41+
* SEP_FILES
42+
* symbol is_enable r_in r_out r_power enhence_a
43+
*
44+
* example
45+
* Li 0
46+
* F 1 F_pbe_50.sep 0.0 2.0 20.0 1.0
47+
*/
48+
int Sep_Cell::read_sep_potentials(std::ifstream& ifpos,
49+
const std::string& pp_dir,
50+
std::ofstream& ofs_running,
51+
UnitCell& ucell)
52+
{
53+
ModuleBase::TITLE("Sep_Cell", "read_sep_potentials");
54+
55+
if (!ModuleBase::GlobalFunc::SCAN_BEGIN(ifpos, "SEP_FILES"))
56+
{
57+
GlobalV::ofs_running << "Cannot find SEP_FILES section in STRU" << std::endl;
58+
return false;
59+
}
60+
61+
ifpos.ignore(300, '\n');
62+
63+
for (int i = 0; i < this->ntype; ++i)
64+
{
65+
std::string one_line, atom_label;
66+
std::getline(ifpos, one_line);
67+
std::stringstream ss(one_line);
68+
69+
// read the label of the atom
70+
bool enable_tmp;
71+
ss >> atom_label >> enable_tmp;
72+
73+
// Validate atom label
74+
if (atom_label != ucell.atom_label[i])
75+
{
76+
GlobalV::ofs_running << "Sep potential and atom order do not match. "
77+
<< "Expected: " << ucell.atoms[i].label << ", Got: " << atom_label << std::endl;
78+
return false;
79+
}
80+
this->sep_enable[i] = enable_tmp;
81+
if (this->sep_enable[i])
82+
{
83+
this->seps[i].is_enable = this->sep_enable[i];
84+
std::string sep_filename;
85+
ss >> sep_filename;
86+
ss >> this->seps[i].r_in >> this->seps[i].r_out >> this->seps[i].r_power >> this->seps[i].enhence_a;
87+
std::string sep_addr = pp_dir + sep_filename;
88+
std::ifstream sep_ifs(sep_addr.c_str(), std::ios::in);
89+
if (!sep_ifs)
90+
{
91+
GlobalV::ofs_running << "Cannot find sep potential file: " << sep_addr << std::endl;
92+
return false;
93+
}
94+
this->seps[i].read_sep(sep_ifs);
95+
}
96+
}
97+
98+
return true;
99+
}
100+
101+
#ifdef __MPI
102+
void Sep_Cell::bcast_sep_cell()
103+
{
104+
ModuleBase::TITLE("Sep_Cell", "bcast_sep_cell");
105+
Parallel_Common::bcast_int(this->ntype);
106+
107+
if (GlobalV::MY_RANK != 0)
108+
{
109+
this->seps.resize(this->ntype);
110+
this->sep_enable.resize(this->ntype);
111+
}
112+
for (int i = 0; i < this->ntype; ++i)
113+
{
114+
bool tmp = false;
115+
if (GlobalV::MY_RANK == 0)
116+
{
117+
tmp = this->sep_enable[i];
118+
}
119+
Parallel_Common::bcast_bool(tmp);
120+
if (GlobalV::MY_RANK != 0)
121+
{
122+
this->sep_enable[i] = tmp;
123+
}
124+
this->seps[i].bcast_sep();
125+
}
126+
}
127+
#endif // __MPI

‎source/source_cell/sep_cell.h‎

Lines changed: 74 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,74 @@
1+
// The Sep_Cell class is container for Sep potential.
2+
3+
#ifndef SEP_CELL
4+
#define SEP_CELL
5+
6+
#include "source_cell/sep.h"
7+
#include "source_cell/unitcell.h"
8+
9+
#include <fstream>
10+
#include <vector>
11+
12+
class Sep_Cell
13+
{
14+
public:
15+
Sep_Cell() noexcept;
16+
~Sep_Cell() noexcept;
17+
18+
// Sets the number of atom types and initializes internal vectors
19+
void init(const int ntype_in);
20+
21+
void set_omega(const double omega_in, const double tpiba2_in);
22+
23+
// Reads self potentials from STRU file and xx.sep files
24+
// Returns true if successful, false otherwise
25+
int read_sep_potentials(std::ifstream& ifpos,
26+
const std::string& pp_dir,
27+
std::ofstream& ofs_running,
28+
UnitCell& ucell);
29+
30+
#ifdef __MPI
31+
// Broadcasts the Sep_Cell object to all processes
32+
void bcast_sep_cell();
33+
#endif // __MPI
34+
35+
// Getter methods
36+
const std::vector<SepPot>& get_seps() const
37+
{
38+
return seps;
39+
}
40+
int get_ntype() const
41+
{
42+
return ntype;
43+
}
44+
const std::vector<bool>& get_sep_enable() const
45+
{
46+
return sep_enable;
47+
}
48+
49+
double get_omega() const
50+
{
51+
return omega;
52+
}
53+
54+
double get_tpiba2() const
55+
{
56+
return tpiba2;
57+
}
58+
59+
private:
60+
std::vector<SepPot> seps; // Self potentials for each atom type
61+
int ntype; // number of atom types
62+
std::vector<bool> sep_enable; // Whether self potential is enabled for each atom type
63+
64+
// unit cell data for VSep
65+
double omega; // unit cell Volume
66+
double tpiba2; // tpiba ^ 2
67+
};
68+
69+
namespace GlobalC
70+
{
71+
extern Sep_Cell sep_cell;
72+
}
73+
74+
#endif // SEP_CEll

0 commit comments

Comments
 (0)