Skip to content

Commit b7e91aa

Browse files
YuLiu98mohanchen
YuLiu98
and
mohanchen
authored
Feature: enable uspp in upf100 format (#4151)
* Feature: uspp in upf100 format * Test: add unittests for uspp upf100 * Refactor: update variable initialization --------- Co-authored-by: Mohan Chen <[email protected]>
1 parent fa7d4a4 commit b7e91aa

File tree

12 files changed

+16402
-8388
lines changed

12 files changed

+16402
-8388
lines changed

‎source/module_cell/read_pp_upf100.cpp

+193-101
Original file line numberDiff line numberDiff line change
@@ -4,75 +4,78 @@
44
// Pseudopot_upfpotential Format
55
int Pseudopot_upf::read_pseudo_upf(std::ifstream &ifs)
66
{
7-
std::string dummy;
8-
int ierr = 0;
9-
this->has_so = false;
7+
std::string dummy;
8+
this->has_so = false;
9+
this->q_with_l = false;
1010

11-
//addinfo_loop
12-
ifs.rdstate();
11+
// addinfo_loop
12+
ifs.rdstate();
1313

14-
while (ifs.good())
15-
{
16-
ifs >> dummy;
17-
if(dummy=="<PP_ADDINFO>")
18-
{
19-
ierr = 1;
20-
break;
21-
}
22-
ifs.rdstate();
23-
}
14+
while (ifs.good())
15+
{
16+
ifs >> dummy;
17+
if (dummy == "<PP_ADDINFO>")
18+
{
19+
this->has_so = true;
20+
}
21+
else if (dummy == "<PP_QIJ_WITH_L>")
22+
{
23+
this->q_with_l = true;
24+
}
2425

25-
if (ierr == 1)
26-
{
27-
has_so = true;
28-
}
29-
30-
// Search for Header
31-
// This version doesn't use the new routine SCAN_BEGIN
32-
// because this search must set extra flags for
33-
// compatibility with other pp format reading
26+
if (has_so && q_with_l)
27+
{
28+
break;
29+
}
30+
ifs.rdstate();
31+
}
3432

35-
ierr = 0;
33+
// Search for Header
34+
// This version doesn't use the new routine SCAN_BEGIN
35+
// because this search must set extra flags for
36+
// compatibility with other pp format reading
3637

37-
ifs.clear();
38-
ifs.seekg(0);
39-
ifs.rdstate();
38+
int ierr = 0;
4039

41-
// header_loop:
42-
while (ifs.good())
43-
{
44-
ifs >> dummy;
45-
if(dummy=="<PP_HEADER>")
46-
{
47-
ierr = 1;
48-
//---------------------
49-
// call member function
50-
//---------------------
51-
read_pseudo_header(ifs);
52-
ModuleBase::GlobalFunc::SCAN_END(ifs, "</PP_HEADER>");
53-
break;
54-
}
55-
}
40+
ifs.clear();
41+
ifs.seekg(0);
42+
ifs.rdstate();
5643

57-
if (ierr == 0)
58-
{
59-
// 2: something in pseudopotential file not match.
60-
return 2;
61-
}
44+
// header_loop:
45+
while (ifs.good())
46+
{
47+
ifs >> dummy;
48+
if (dummy == "<PP_HEADER>")
49+
{
50+
ierr = 1;
51+
//---------------------
52+
// call member function
53+
//---------------------
54+
read_pseudo_header(ifs);
55+
ModuleBase::GlobalFunc::SCAN_END(ifs, "</PP_HEADER>");
56+
break;
57+
}
58+
}
6259

63-
// Search for mesh information
64-
if ( ModuleBase::GlobalFunc::SCAN_BEGIN(ifs, "<PP_MESH>") )
65-
{
66-
//---------------------
67-
// call member function
68-
//---------------------
69-
read_pseudo_mesh(ifs);
70-
ModuleBase::GlobalFunc::SCAN_END(ifs, "</PP_MESH>");
71-
}
60+
if (ierr == 0)
61+
{
62+
// 2: something in pseudopotential file not match.
63+
return 2;
64+
}
7265

73-
// If present, search for nlcc
74-
if (this->nlcc)
75-
{
66+
// Search for mesh information
67+
if (ModuleBase::GlobalFunc::SCAN_BEGIN(ifs, "<PP_MESH>"))
68+
{
69+
//---------------------
70+
// call member function
71+
//---------------------
72+
read_pseudo_mesh(ifs);
73+
ModuleBase::GlobalFunc::SCAN_END(ifs, "</PP_MESH>");
74+
}
75+
76+
// If present, search for nlcc
77+
if (this->nlcc)
78+
{
7679
ModuleBase::GlobalFunc::SCAN_BEGIN(ifs, "<PP_NLCC>");
7780
//---------------------
7881
// call member function
@@ -148,29 +151,31 @@ void Pseudopot_upf::read_pseudo_header(std::ifstream &ifs)
148151
if(pp_type=="US")
149152
{
150153
this->tvanp = true;
151-
}
152-
else if(pp_type=="NC")
153-
{
154-
this->tvanp = false;
155-
}
156-
else if(pp_type == "1/r")
157-
{
158-
this->tvanp = false;
159-
this->coulomb_potential = true;
160-
}
161-
else
162-
{
163-
// A bug here!!! can't quit together.
164-
std::cout << " pp_type=" << pp_type << std::endl;
165-
ModuleBase::WARNING_QUIT("Pseudopot_upf::read_pseudo_header","unknown pseudo type");
166-
}
154+
this->coulomb_potential = false;
155+
}
156+
else if (pp_type == "NC")
157+
{
158+
this->tvanp = false;
159+
this->coulomb_potential = false;
160+
}
161+
else if (pp_type == "1/r")
162+
{
163+
this->tvanp = false;
164+
this->coulomb_potential = true;
165+
}
166+
else
167+
{
168+
// A bug here!!! can't quit together.
169+
std::cout << " pp_type=" << pp_type << std::endl;
170+
ModuleBase::WARNING_QUIT("Pseudopot_upf::read_pseudo_header", "unknown pseudo type");
171+
}
167172

168-
// If use nlcc
169-
std::string nlc;
170-
ModuleBase::GlobalFunc::READ_VALUE(ifs, nlc);
173+
// If use nlcc
174+
std::string nlc;
175+
ModuleBase::GlobalFunc::READ_VALUE(ifs, nlc);
171176

172-
if (nlc == "T")
173-
{
177+
if (nlc == "T")
178+
{
174179
this->nlcc = true;
175180
}
176181
else
@@ -282,28 +287,28 @@ void Pseudopot_upf::read_pseudo_local(std::ifstream &ifs)
282287

283288
void Pseudopot_upf::read_pseudo_nl(std::ifstream &ifs)
284289
{
285-
// int nb, mb, n, ir, idum, ldum, lp, i, ikk;
286-
int nb, mb, ir, idum;
287-
288-
if (nbeta == 0)
289-
{
290-
delete[] kbeta;
291-
delete[] lll;
292-
this->kbeta = nullptr;
293-
this->lll = nullptr;
294-
this->beta.create(1, 1);
295-
this->dion.create(1, 1);
296-
kkbeta = 0;
290+
// int nb, mb, n, ir, idum, ldum, lp, i, ikk;
291+
int nb = 0;
292+
int mb = 0;
293+
int ir = 0;
294+
int idum = 0;
295+
int ldum = 0;
296+
297+
if (nbeta == 0)
298+
{
299+
this->nqf = 0;
300+
this->nqlc = 0;
301+
this->kkbeta = 0;
297302
return;
298303
}
299304
else
300305
{
301306
delete[] kbeta;
302307
delete[] lll;
303308
this->kbeta = new int[nbeta];
304-
this->lll = new int[nbeta];
305-
this->beta.create(nbeta , mesh);
306-
this->dion.create(nbeta , nbeta);
309+
this->lll = new int[nbeta];
310+
this->beta.create(nbeta, mesh);
311+
this->dion.create(nbeta, nbeta);
307312
kkbeta = 0;
308313

309314
for (int i = 0; i < nbeta; i++)
@@ -341,11 +346,98 @@ void Pseudopot_upf::read_pseudo_nl(std::ifstream &ifs)
341346
// QIJ
342347
if (tvanp)
343348
{
344-
ModuleBase::WARNING_QUIT("read_pseudo_nl","Ultra Soft Pseudopotential not available yet.");
345-
}
346-
else // not tvanp
347-
{
348-
}
349+
if (!ModuleBase::GlobalFunc::SCAN_BEGIN(ifs, "<PP_QIJ>", false))
350+
{
351+
ModuleBase::GlobalFunc::SCAN_BEGIN(ifs, "<PP_QIJ_WITH_L>", false);
352+
}
353+
// If nqf is not zero, Qij's inside rinner are computed using qfcoef's
354+
ModuleBase::GlobalFunc::READ_VALUE(ifs, this->nqf);
355+
this->nqlc = 2 * this->lmax + 1;
356+
delete[] rinner;
357+
this->rinner = new double[nqlc];
358+
qqq.create(nbeta, nbeta);
359+
if (q_with_l)
360+
{
361+
this->qfuncl.create(2 * lmax + 1, nbeta * (nbeta + 1) / 2, mesh);
362+
}
363+
else
364+
{
365+
this->qfunc.create(nbeta * (nbeta + 1) / 2, mesh);
366+
}
367+
368+
if (nqf <= 0)
369+
{
370+
ModuleBase::GlobalFunc::ZEROS(rinner, nqlc);
371+
this->qfcoef.create(1, 1, 1, 1);
372+
}
373+
else
374+
{
375+
this->qfcoef.create(nbeta, nbeta, nqlc, nqf);
376+
ModuleBase::GlobalFunc::SCAN_BEGIN(ifs, "<PP_RINNER>", false);
377+
for (int i = 0; i < nqlc; i++)
378+
{
379+
ifs >> idum >> rinner[i];
380+
}
381+
ModuleBase::GlobalFunc::SCAN_END(ifs, "</PP_RINNER>");
382+
}
383+
384+
for (int nb = 0; nb < nbeta; nb++)
385+
{
386+
int ln = lll[nb];
387+
for (int mb = nb; mb < nbeta; mb++)
388+
{
389+
int lm = lll[mb];
390+
int nmb = mb * (mb + 1) / 2 + nb;
391+
ifs >> idum >> idum >> ldum; // i j (l(j))
392+
ifs.ignore(75, '\n');
393+
394+
if (ldum != lm)
395+
{
396+
ModuleBase::WARNING_QUIT("Pseudopot_upf::read_pseudo_nl",
397+
"inconsistent angular momentum for Q_ij");
398+
}
399+
400+
ModuleBase::GlobalFunc::READ_VALUE(ifs, this->qqq(nb, mb));
401+
this->qqq(mb, nb) = this->qqq(nb, mb);
402+
403+
if (q_with_l)
404+
{
405+
for (int l = std::abs(ln - lm); l <= ln + lm; l += 2)
406+
{
407+
for (int ir = 0; ir < mesh; ir++)
408+
{
409+
ifs >> qfuncl(l, nmb, ir);
410+
}
411+
}
412+
}
413+
else
414+
{
415+
for (int ir = 0; ir < mesh; ir++)
416+
{
417+
ifs >> qfunc(nmb, ir);
418+
}
419+
}
420+
421+
if (this->nqf > 0)
422+
{
423+
ModuleBase::GlobalFunc::SCAN_BEGIN(ifs, "<PP_QFCOEF>", false);
424+
for (int k = 0; k < nqlc; k++)
425+
{
426+
for (int l = 0; l < nqf; l++)
427+
{
428+
ifs >> qfcoef(nb, mb, k, l);
429+
qfcoef(mb, nb, k, l) = qfcoef(nb, mb, k, l);
430+
}
431+
}
432+
ModuleBase::GlobalFunc::SCAN_END(ifs, "</PP_QFCOEF>");
433+
}
434+
}
435+
}
436+
ModuleBase::GlobalFunc::SCAN_END(ifs, "</PP_QIJ>");
437+
}
438+
else // not tvanp
439+
{
440+
}
349441
}
350442
return;
351443
}

0 commit comments

Comments
 (0)