1+ #include < cmath>
2+ #include < vector>
3+ #include < map>
4+ #include < tuple>
5+ #include < complex>
6+ #include < fstream>
7+ #include < memory>
8+ #include " module_cell/unitcell.h"
9+ #include " module_base/spherical_bessel_transformer.h"
10+ #include " module_basis/module_nao/two_center_integrator.h"
11+ #include " module_cell/module_neighbor/sltk_grid_driver.h"
12+ #include " module_cell/module_neighbor/sltk_atom_arrange.h"
13+ #include " module_parameter/parameter.h"
14+ #include " module_io/cal_pLpR.h"
15+ #include " module_base/formatter.h"
16+ #include " module_base/parallel_common.h"
17+ /* *
18+ *
19+ * FIXME: the following part will be transfered to TwoCenterIntegrator soon
20+ *
21+ */
22+
23+ // L+|l, m> = sqrt((l-m)(l+m+1))|l, m+1>, return the sqrt((l-m)(l+m+1))
24+ double _lplus_on_ylm (const int l, const int m)
25+ {
26+ return std::sqrt ((l - m) * (l + m + 1 ));
27+ }
28+
29+ // L-|l, m> = sqrt((l+m)(l-m+1))|l, m-1>, return the sqrt((l+m)(l-m+1))
30+ double _lminus_on_ylm (const int l, const int m)
31+ {
32+ return std::sqrt ((l + m) * (l - m + 1 ));
33+ }
34+
35+ std::complex <double > ModuleIO::cal_LzijR (
36+ const std::unique_ptr<TwoCenterIntegrator>& calculator,
37+ const int it, const int ia, const int il, const int iz, const int mi,
38+ const int jt, const int ja, const int jl, const int jz, const int mj,
39+ const ModuleBase::Vector3<double >& vR)
40+ {
41+ double val_ = 0 ;
42+ calculator->calculate (it, il, iz, mi, jt, jl, jz, mj, vR, &val_);
43+ return std::complex <double >(mi) * val_;
44+ }
45+
46+ std::complex <double > ModuleIO::cal_LyijR (
47+ const std::unique_ptr<TwoCenterIntegrator>& calculator,
48+ const int it, const int ia, const int il, const int iz, const int im,
49+ const int jt, const int ja, const int jl, const int jz, const int jm,
50+ const ModuleBase::Vector3<double >& vR)
51+ {
52+ // Ly = -i/2 * (L+ - L-)
53+ const double plus_ = _lplus_on_ylm (jl, jm);
54+ const double minus_ = _lminus_on_ylm (jl, jm);
55+ double val_plus = 0 , val_minus = 0 ;
56+ if (plus_ != 0 )
57+ {
58+ calculator->calculate (it, il, iz, im, jt, jl, jz, jm + 1 , vR, &val_plus);
59+ val_plus *= plus_;
60+ }
61+ if (minus_ != 0 )
62+ {
63+ calculator->calculate (it, il, iz, im, jt, jl, jz, jm - 1 , vR, &val_minus);
64+ val_minus *= minus_;
65+ }
66+ return std::complex <double >(0 , -0.5 ) * (val_plus - val_minus);
67+ }
68+
69+ std::complex <double > ModuleIO::cal_LxijR (
70+ const std::unique_ptr<TwoCenterIntegrator>& calculator,
71+ const int it, const int ia, const int il, const int iz, const int im,
72+ const int jt, const int ja, const int jl, const int jz, const int jm,
73+ const ModuleBase::Vector3<double >& vR)
74+ {
75+ // Lx = 1/2 * (L+ + L-)
76+ const double plus_ = _lplus_on_ylm (jl, jm);
77+ const double minus_ = _lminus_on_ylm (jl, jm);
78+ double val_plus = 0 , val_minus = 0 ;
79+ if (plus_ != 0 )
80+ {
81+ calculator->calculate (it, il, iz, im, jt, jl, jz, jm + 1 , vR, &val_plus);
82+ val_plus *= plus_;
83+ }
84+ if (minus_ != 0 )
85+ {
86+ calculator->calculate (it, il, iz, im, jt, jl, jz, jm - 1 , vR, &val_minus);
87+ val_minus *= minus_;
88+ }
89+ return std::complex <double >(0.5 ) * (val_plus + val_minus);
90+ }
91+
92+ ModuleIO::AngularMomentumCalculator::AngularMomentumCalculator (
93+ const std::string& orbital_dir,
94+ const UnitCell& ucell,
95+ const double & search_radius,
96+ const int tdestructor,
97+ const int tgrid,
98+ const int tatom,
99+ const bool searchpbc,
100+ std::ofstream* ptr_log,
101+ const int rank)
102+ {
103+
104+ // ofs_running
105+ this ->ofs_ = ptr_log;
106+ *ofs_ << " \n\n\n\n "
107+ *ofs_ << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>"
108+ *ofs_ << " | |"
109+ *ofs_ << " | Angular momentum expectation value calculation: |"
110+ *ofs_ << " | This is a post-processing step. The expectation value of operator |"
111+ *ofs_ << " | Lx, Ly, Lz (<a|L|b>, in which a and b are ABACUS numerical atomic |"
112+ *ofs_ << " | orbitals) will be calculated. |"
113+ *ofs_ << " | The result will be printed to file with name ${suffix}_Lx/y/z.dat |"
114+ *ofs_ << " | |"
115+ *ofs_ << " <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<"
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 )
124+ {
125+ for (int i = 0 ; i < ucell.ntype ; ++i)
126+ {
127+ forb[i] = orbital_dir + ucell.orbital_fn [i];
128+ }
129+ }
130+ #ifdef __MPI
131+ Parallel_Common::bcast_string (forb.data (), ntype_);
132+ #endif
133+
134+ this ->orb_ = std::unique_ptr<RadialCollection>(new RadialCollection);
135+ this ->orb_ ->build (ucell.ntype , forb.data (), ' o'
136+
137+ ModuleBase::SphericalBesselTransformer sbt (true );
138+ this ->orb_ ->set_transformer (sbt);
139+
140+ const double rcut_max = orb_->rcut_max ();
141+ const int ngrid = int (rcut_max / 0.01 ) + 1 ;
142+ const double cutoff = 2.0 * rcut_max;
143+ this ->orb_ ->set_uniform_grid (true , ngrid, cutoff, ' i' true );
144+
145+ this ->calculator_ = std::unique_ptr<TwoCenterIntegrator>(new TwoCenterIntegrator);
146+ this ->calculator_ ->tabulate (*orb_, *orb_, ' S'
147+
148+ // Initialize Ylm coefficients
149+ ModuleBase::Ylm::set_coefficients ();
150+
151+ // 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);
159+ this ->neighbor_searcher_ = std::unique_ptr<Grid_Driver>(new Grid_Driver (tdestructor, tgrid));
160+ atom_arrange::search (searchpbc,
161+ *ofs_,
162+ *neighbor_searcher_,
163+ ucell,
164+ temp,
165+ tatom);
166+ }
167+
168+ void ModuleIO::AngularMomentumCalculator::kernel (
169+ std::ofstream* ofs,
170+ const UnitCell& ucell,
171+ const char dir,
172+ const int precision)
173+ {
174+ if (!ofs->is_open ())
175+ {
176+ return ;
177+ }
178+ // an easy sanity check
179+ assert (dir == ' x' ' y' ' z'
180+
181+ // it, ia, il, iz, im, iRx, iRy, iRz, jt, ja, jl, jz, jm
182+ // the iRx, iRy, iRz are the indices of the supercell in which the two-center-integral
183+ // it and jt are indexes of atomtypes,
184+ // ia and ja are indexes of atoms within the atomtypes,
185+ // il and jl are indexes of the angular momentum,
186+ // iz and jz are indexes of the zeta functions
187+ // im and jm are indexes of the magnetic quantum numbers.
188+ std::string fmtstr = " %4d%4d%4d%4d%4d%4d%4d%4d%4d%4d%4d%4d%4d"
189+ fmtstr += " %" std::to_string (precision*2 ) + " ." std::to_string (precision) + " e"
190+ fmtstr += " %" std::to_string (precision*2 ) + " ." std::to_string (precision) + " e\n "
191+ FmtCore fmt (fmtstr);
192+
193+ ModuleBase::Vector3<double > ri, rj, dr;
194+ for (int it = 0 ; it < ucell.ntype ; it++)
195+ {
196+ const Atom& atyp_i = ucell.atoms [it];
197+ for (int ia = 0 ; ia < atyp_i.na ; ia++)
198+ {
199+ ri = atyp_i.tau [ia];
200+ neighbor_searcher_->Find_atom (ucell, ri, it, ia);
201+ for (int ia_adj = 0 ; ia_adj < neighbor_searcher_->getAdjacentNum (); ia_adj++)
202+ {
203+ rj = neighbor_searcher_->getAdjacentTau (ia_adj);
204+ int jt = neighbor_searcher_->getType (ia_adj);
205+ const Atom& atyp_j = ucell.atoms [jt];
206+ int ja = neighbor_searcher_->getNatom (ia_adj);
207+ dr = (ri - rj) * ucell.lat0 ;
208+ const ModuleBase::Vector3<int > iR = neighbor_searcher_->getBox (ia_adj);
209+ // the two-center-integral
210+
211+ for (int li = 0 ; li < atyp_i.nwl + 1 ; li++)
212+ {
213+ for (int iz = 0 ; iz < atyp_i.l_nchi [li]; iz++)
214+ {
215+ for (int mi = -li; mi <= li; mi++)
216+ {
217+ for (int lj = 0 ; lj < atyp_j.nwl + 1 ; lj++)
218+ {
219+ for (int jz = 0 ; jz < atyp_j.l_nchi [lj]; jz++)
220+ {
221+ for (int mj = -lj; mj <= lj; mj++)
222+ {
223+ std::complex <double > val = 0 ;
224+ if (dir == ' x'
225+ {
226+ val = cal_LxijR (calculator_,
227+ it, ia, li, iz, mi, jt, ja, lj, jz, mj, dr);
228+ }
229+ else if (dir == ' y'
230+ {
231+ val = cal_LyijR (calculator_,
232+ it, ia, li, iz, mi, jt, ja, lj, jz, mj, dr);
233+ }
234+ else if (dir == ' z'
235+ {
236+ val = cal_LzijR (calculator_,
237+ it, ia, li, iz, mi, jt, ja, lj, jz, mj, dr);
238+ }
239+
240+ *ofs << fmt.format (
241+ it, ia, li, iz, mi,
242+ iR.x , iR.y , iR.z ,
243+ jt, ja, lj, jz, mj,
244+ val.real (), val.imag ());
245+ }
246+ }
247+ }
248+ }
249+ }
250+ }
251+ }
252+ }
253+ }
254+ }
255+
256+ void ModuleIO::AngularMomentumCalculator::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" " .dat"
286+ ofout.open (fn, std::ios::out);
287+ ofout << title;
288+ this ->kernel (&ofout, ucell, d, precision);
289+ ofout.close ();
290+ }
291+ }
0 commit comments