Skip to content

Commit 4c6e7c8

Browse files
kirk0830kirk0830
committed
complete
1 parent 43d1c0e commit 4c6e7c8

File tree

4 files changed

+193
-31
lines changed

4 files changed

+193
-31
lines changed

‎source/module_esolver/esolver_ks_lcao.cpp‎

Lines changed: 7 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1253,9 +1253,14 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(const int istep)
12531253
PARAM.inp.test_grid,
12541254
PARAM.inp.test_atom_input,
12551255
PARAM.inp.search_pbc,
1256-
&GlobalV::ofs_running
1256+
&GlobalV::ofs_running,
1257+
GlobalV::MY_RANK
12571258
);
1258-
1259+
mylcalculator.calculate(PARAM.inp.suffix,
1260+
PARAM.inp.read_file_dir,
1261+
GlobalC::ucell,
1262+
PARAM.inp.out_mat_l[1],
1263+
GlobalV::MY_RANK);
12591264
}
12601265

12611266
// add by jingan in 2018.11.7

‎source/module_io/cal_pLpR.cpp‎

Lines changed: 91 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -13,6 +13,7 @@
1313
#include "module_parameter/parameter.h"
1414
#include "module_io/cal_pLpR.h"
1515
#include "module_base/formatter.h"
16+
#include "module_base/parallel_common.h"
1617
/**
1718
*
1819
* FIXME: the following part will be transfered to TwoCenterIntegrator soon
@@ -96,51 +97,81 @@ ModuleIO::AngularMomentumExpectationCalculator::AngularMomentumExpectationCalcul
9697
const int tgrid,
9798
const int tatom,
9899
const bool searchpbc,
99-
std::ofstream* ptr_log)
100+
std::ofstream* ptr_log,
101+
const int rank)
100102
{
101-
std::vector<std::string> forb(ucell.ntype);
102-
for (int i = 0; i < ucell.ntype; ++i)
103+
104+
// ofs_running
105+
this->ofs_ = ptr_log;
106+
*ofs_ << "\n\n\n\n";
107+
*ofs_ << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>" << std::endl;
108+
*ofs_ << " | |" << std::endl;
109+
*ofs_ << " | Angular momentum expectation value calculation: |" << std::endl;
110+
*ofs_ << " | This is a post-processing step. The expectation value of operator |" << std::endl;
111+
*ofs_ << " | Lx, Ly, Lz (<a|L|b>, in which a and b are ABACUS numerical atomic |" << std::endl;
112+
*ofs_ << " | orbitals) will be calculated. |" << std::endl;
113+
*ofs_ << " | The result will be printed to file with name ${suffix}_Lx/y/z.dat |" << std::endl;
114+
*ofs_ << " | |" << std::endl;
115+
*ofs_ << " <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<" << std::endl;
116+
*ofs_ << "\n\n\n\n";
117+
118+
int ntype_ = ucell.ntype;
119+
#ifdef __MPI
120+
Parallel_Common::bcast_int(ntype_);
121+
#endif
122+
std::vector<std::string> forb(ntype_);
123+
if (rank == 0)
103124
{
104-
forb[i] = orbital_dir + ucell.orbital_fn[i];
125+
for (int i = 0; i < ucell.ntype; ++i)
126+
{
127+
forb[i] = orbital_dir + ucell.orbital_fn[i];
128+
}
105129
}
130+
#ifdef __MPI
131+
Parallel_Common::bcast_string(forb.data(), ntype_);
132+
#endif
133+
106134
this->orb_ = std::unique_ptr<RadialCollection>(new RadialCollection);
107135
this->orb_->build(ucell.ntype, forb.data(), 'o');
108-
136+
109137
ModuleBase::SphericalBesselTransformer sbt(true);
110138
this->orb_->set_transformer(sbt);
111-
139+
112140
const double rcut_max = orb_->rcut_max();
113141
const int ngrid = int(rcut_max / 0.01) + 1;
114142
const double cutoff = 2.0 * rcut_max;
115143
this->orb_->set_uniform_grid(true, ngrid, cutoff, 'i', true);
116-
144+
117145
this->calculator_ = std::unique_ptr<TwoCenterIntegrator>(new TwoCenterIntegrator);
118146
this->calculator_->tabulate(*orb_, *orb_, 'S', ngrid, cutoff);
119-
147+
120148
// Initialize Ylm coefficients
121149
ModuleBase::Ylm::set_coefficients();
122-
123-
// ofs_running
124-
this->ofs_ = ptr_log;
125-
150+
126151
// for neighbor list search
152+
double temp = -1.0;
153+
temp = atom_arrange::set_sr_NL(*ofs_,
154+
PARAM.inp.out_level,
155+
search_radius,
156+
ucell.infoNL.get_rcutmax_Beta(),
157+
PARAM.globalv.gamma_only_local);
158+
temp = std::max(temp, search_radius);
127159
this->neighbor_searcher_ = std::unique_ptr<Grid_Driver>(new Grid_Driver(tdestructor, tgrid));
128-
atom_arrange::search(
129-
searchpbc,
130-
*ofs_,
131-
*neighbor_searcher_,
132-
ucell,
133-
search_radius,
134-
tatom);
160+
atom_arrange::search(searchpbc,
161+
*ofs_,
162+
*neighbor_searcher_,
163+
ucell,
164+
temp,
165+
tatom);
135166
}
136167

137-
void ModuleIO::AngularMomentumExpectationCalculator::calculate(
168+
void ModuleIO::AngularMomentumExpectationCalculator::kernel(
138169
std::ofstream* ofs,
139170
const UnitCell& ucell,
140171
const char dir,
141172
const int precision)
142173
{
143-
if (ofs == nullptr)
174+
if (!ofs->is_open())
144175
{
145176
return;
146177
}
@@ -205,7 +236,8 @@ void ModuleIO::AngularMomentumExpectationCalculator::calculate(
205236
val = cal_LzijR(calculator_,
206237
it, ia, li, iz, mi, jt, ja, lj, jz, mj, dr);
207238
}
208-
*ofs_ << fmt.format(
239+
240+
*ofs << fmt.format(
209241
it, ia, li, iz, mi,
210242
iR.x, iR.y, iR.z,
211243
jt, ja, lj, jz, mj,
@@ -219,4 +251,41 @@ void ModuleIO::AngularMomentumExpectationCalculator::calculate(
219251
}
220252
}
221253
}
254+
}
255+
256+
void ModuleIO::AngularMomentumExpectationCalculator::calculate(
257+
const std::string& prefix,
258+
const std::string& outdir,
259+
const UnitCell& ucell,
260+
const int precision,
261+
const int rank)
262+
{
263+
if (rank != 0)
264+
{
265+
return;
266+
}
267+
std::ofstream ofout;
268+
const std::string dir = "xyz";
269+
const std::string title = "# it ia il iz im iRx iRy iRz jt ja jl jz jm <a|L|b>\n"
270+
"# it: atomtype index of the first atom\n"
271+
"# ia: atomic index of the first atom within the atomtype\n"
272+
"# il: angular momentum index of the first atom\n"
273+
"# iz: zeta function index of the first atom\n"
274+
"# im: magnetic quantum number of the first atom\n"
275+
"# iRx, iRy, iRz: the indices of the supercell\n"
276+
"# jt: atomtype index of the second atom\n"
277+
"# ja: atomic index of the second atom within the atomtype\n"
278+
"# jl: angular momentum index of the second atom\n"
279+
"# jz: zeta function index of the second atom\n"
280+
"# jm: magnetic quantum number of the second atom\n"
281+
"# <a|L|b>: the value of the matrix element\n";
282+
283+
for (char d : dir)
284+
{
285+
std::string fn = outdir + prefix + "_L" + d + ".dat";
286+
ofout.open(fn, std::ios::out);
287+
ofout << title;
288+
this->kernel(&ofout, ucell, d, precision);
289+
ofout.close();
290+
}
222291
}

‎source/module_io/cal_pLpR.h‎

Lines changed: 94 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -71,18 +71,75 @@
7171

7272
namespace ModuleIO
7373
{
74+
/**
75+
* @brief calculate the <phi_i|Lz|phi_j> matrix elements, in which the Lz
76+
* are the angular momentum operators, |phi_i> and |phi_j> are the numerical
77+
* atomic orbitals (NAOs).
78+
*
79+
* @param calculator the std::unique_ptr<TwoCenterIntegrator> instance
80+
* @param it atomtype index of the first atom
81+
* @param ia atomic index of the first atom within the atomtype
82+
* @param il angular momentum index of the first atom
83+
* @param iz zeta function index of the first atom
84+
* @param mi magnetic quantum number of the first atom
85+
* @param jt atomtype index of the second atom
86+
* @param ja atomic index of the second atom within the atomtype
87+
* @param jl angular momentum index of the second atom
88+
* @param jz zeta function index of the second atom
89+
* @param mj magnetic quantum number of the second atom
90+
* @param vR the vector from the first atom to the second atom
91+
* @return std::complex<double>
92+
*/
7493
std::complex<double> cal_LzijR(
7594
const std::unique_ptr<TwoCenterIntegrator>& calculator,
7695
const int it, const int ia, const int il, const int iz, const int mi,
7796
const int jt, const int ja, const int jl, const int jz, const int mj,
7897
const ModuleBase::Vector3<double>& vR);
7998

99+
/**
100+
* @brief calculate the <phi_i|Ly|phi_j> matrix elements, in which the Lz
101+
* are the angular momentum operators, |phi_i> and |phi_j> are the numerical
102+
* atomic orbitals (NAOs).
103+
*
104+
* @param calculator the std::unique_ptr<TwoCenterIntegrator> instance
105+
* @param it atomtype index of the first atom
106+
* @param ia atomic index of the first atom within the atomtype
107+
* @param il angular momentum index of the first atom
108+
* @param iz zeta function index of the first atom
109+
* @param mi magnetic quantum number of the first atom
110+
* @param jt atomtype index of the second atom
111+
* @param ja atomic index of the second atom within the atomtype
112+
* @param jl angular momentum index of the second atom
113+
* @param jz zeta function index of the second atom
114+
* @param mj magnetic quantum number of the second atom
115+
* @param vR the vector from the first atom to the second atom
116+
* @return std::complex<double>
117+
*/
80118
std::complex<double> cal_LyijR(
81119
const std::unique_ptr<TwoCenterIntegrator>& calculator,
82120
const int it, const int ia, const int il, const int iz, const int im,
83121
const int jt, const int ja, const int jl, const int jz, const int jm,
84122
const ModuleBase::Vector3<double>& vR);
85123

124+
/**
125+
* @brief calculate the <phi_i|Lx|phi_j> matrix elements, in which the Lz
126+
* are the angular momentum operators, |phi_i> and |phi_j> are the numerical
127+
* atomic orbitals (NAOs).
128+
*
129+
* @param calculator the std::unique_ptr<TwoCenterIntegrator> instance
130+
* @param it atomtype index of the first atom
131+
* @param ia atomic index of the first atom within the atomtype
132+
* @param il angular momentum index of the first atom
133+
* @param iz zeta function index of the first atom
134+
* @param mi magnetic quantum number of the first atom
135+
* @param jt atomtype index of the second atom
136+
* @param ja atomic index of the second atom within the atomtype
137+
* @param jl angular momentum index of the second atom
138+
* @param jz zeta function index of the second atom
139+
* @param mj magnetic quantum number of the second atom
140+
* @param vR the vector from the first atom to the second atom
141+
* @return std::complex<double>
142+
*/
86143
std::complex<double> cal_LxijR(
87144
const std::unique_ptr<TwoCenterIntegrator>& calculator,
88145
const int it, const int ia, const int il, const int iz, const int im,
@@ -103,7 +160,20 @@ namespace ModuleIO
103160
class AngularMomentumExpectationCalculator
104161
{
105162
public:
106-
AngularMomentumExpectationCalculator() = delete; // the default constructor is meaningless
163+
// the default constructor is meaningless
164+
AngularMomentumExpectationCalculator() = delete;
165+
/**
166+
* @brief Construct a new Angular Momentum Expectation Calculator object
167+
*
168+
* @param orbital_dir the directory of the orbital file
169+
* @param ucell the unit cell object
170+
* @param search_radius the search radius for the neighboring atoms
171+
* @param tdestructor test flag, for destructor
172+
* @param tgrid test flag, for grid
173+
* @param tatom test flag, for atom input
174+
* @param searchpbc
175+
* @param ptr_log pointer to the ofstream object for logging
176+
*/
107177
AngularMomentumExpectationCalculator(
108178
const std::string& orbital_dir,
109179
const UnitCell& ucell,
@@ -112,13 +182,15 @@ namespace ModuleIO
112182
const int tgrid,
113183
const int tatom,
114184
const bool searchpbc,
115-
std::ofstream* ptr_log = nullptr);
185+
std::ofstream* ptr_log = nullptr,
186+
const int rank = 0);
116187
~AngularMomentumExpectationCalculator() = default;
117188

118-
void calculate(std::ofstream* ofs,
119-
const UnitCell& ucell,
120-
const char dir = 'x',
121-
const int precision = 10);
189+
void calculate(const std::string& prefix,
190+
const std::string& outdir,
191+
const UnitCell& ucell,
192+
const int precision = 10,
193+
const int rank = 0);
122194

123195
private:
124196
// ofsrunning
@@ -132,5 +204,21 @@ namespace ModuleIO
132204

133205
// neighboring searcher
134206
std::unique_ptr<Grid_Driver> neighbor_searcher_;
207+
208+
/**
209+
* @brief calculate the <phi_i|Lx/Ly/Lz|phi_j> matrix elements. Due to
210+
* the large size of the matrix, the result will be printed to file
211+
* directly.
212+
*
213+
* @param ofs pointer to the ofstream object for printing, if nullptr,
214+
* the result will not be printed
215+
* @param ucell the unit cell object
216+
* @param dir the direction of the angular momentum operator, 'x', 'y' or 'z'
217+
* @param precision the precision of the output, default is 10
218+
*/
219+
void kernel(std::ofstream* ofs,
220+
const UnitCell& ucell,
221+
const char dir = 'x',
222+
const int precision = 10);
135223
};
136224
} // namespace ModuleIO

‎source/module_io/test/CMakeLists.txt‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -242,7 +242,7 @@ add_test(NAME orb_io_test_parallel
242242
if(ENABLE_LCAO)
243243
AddTest(
244244
TARGET cal_pLpR_test
245-
LIBS parameter base ${math_libs} device
245+
LIBS parameter base ${math_libs} device neighbor
246246
SOURCES
247247
cal_pLpR_test.cpp
248248
../cal_pLpR.cpp

0 commit comments

Comments
 (0)