Skip to content

Commit 2505d4f

Browse files
mohanchenmohanchen
authored
Some small updates for the format (#5972)
* reorganize esolver_ks_lcao.h * delete some useless files in esolver_ks.h * change pelec->charge to this->chr * will delete pelec->charge in near future, we just directly use chr defined in esolver_fp * fix problem in ML_KEDF * update this->chr * fix some format issues * update cal_ux.cpp * update fp * small updates
1 parent 4f72eb8 commit 2505d4f

File tree

7 files changed

+144
-75
lines changed

7 files changed

+144
-75
lines changed

‎source/module_elecstate/cal_ux.cpp‎

Lines changed: 71 additions & 42 deletions
Original file line numberDiff line numberDiff line change
@@ -4,62 +4,91 @@
44
namespace elecstate {
55

66
void cal_ux(UnitCell& ucell) {
7+
78
if (PARAM.inp.nspin != 4)
89
{
910
return;
1011
}
1112

12-
double amag, uxmod;
13+
const double absolute_mag_thr = 1.0e-6;
14+
15+
double amag = 0.0;
16+
double uxmod = 0.0;
17+
1318
int starting_it = 0;
1419
int starting_ia = 0;
15-
bool is_paraller;
20+
bool is_paraller = false;
21+
1622
// do not sign feature in teh general case
1723
ucell.magnet.lsign_ = false;
1824
ModuleBase::GlobalFunc::ZEROS(ucell.magnet.ux_, 3);
1925

20-
for (int it = 0; it < ucell.ntype; it++) {
21-
for (int ia = 0; ia < ucell.atoms[it].na; ia++) {
26+
for (int it = 0; it < ucell.ntype; it++)
27+
{
28+
for (int ia = 0; ia < ucell.atoms[it].na; ia++)
29+
{
30+
// m_loc_: local magnetization vector for each atom
2231
amag = pow(ucell.atoms[it].m_loc_[ia].x, 2)
2332
+ pow(ucell.atoms[it].m_loc_[ia].y, 2)
2433
+ pow(ucell.atoms[it].m_loc_[ia].z, 2);
25-
if (amag > 1e-6) {
26-
ucell.magnet.ux_[0] = ucell.atoms[it].m_loc_[ia].x;
27-
ucell.magnet.ux_[1] = ucell.atoms[it].m_loc_[ia].y;
28-
ucell.magnet.ux_[2] = ucell.atoms[it].m_loc_[ia].z;
29-
starting_it = it;
30-
starting_ia = ia;
31-
ucell.magnet.lsign_ = true;
32-
break;
33-
}
34-
}
35-
if (ucell.magnet.lsign_) {
36-
break;
37-
}
38-
}
34+
35+
// find the first atom (it,ia) whose magnetism is not zero
36+
// compute ux
37+
if (amag > absolute_mag_thr)
38+
{
39+
ucell.magnet.ux_[0] = ucell.atoms[it].m_loc_[ia].x;
40+
ucell.magnet.ux_[1] = ucell.atoms[it].m_loc_[ia].y;
41+
ucell.magnet.ux_[2] = ucell.atoms[it].m_loc_[ia].z;
42+
43+
starting_it = it;
44+
starting_ia = ia;
45+
46+
ucell.magnet.lsign_ = true;
47+
break;
48+
}
49+
}
50+
51+
// if any atom has magnetism, then break the for iteration
52+
if (ucell.magnet.lsign_)
53+
{
54+
break;
55+
}
56+
}
57+
3958
// whether the initial magnetizations is parallel
40-
for (int it = starting_it; it < ucell.ntype; it++) {
41-
for (int ia = 0; ia < ucell.atoms[it].na; ia++) {
42-
if (it > starting_it || ia > starting_ia) {
43-
ucell.magnet.lsign_
44-
= ucell.magnet.lsign_
45-
&& judge_parallel(ucell.magnet.ux_, ucell.atoms[it].m_loc_[ia]);
46-
}
47-
}
48-
}
49-
if (ucell.magnet.lsign_) {
50-
uxmod = pow(ucell.magnet.ux_[0], 2) + pow(ucell.magnet.ux_[1], 2)
51-
+ pow(ucell.magnet.ux_[2], 2);
52-
if (uxmod < 1e-6) {
53-
ModuleBase::WARNING_QUIT("cal_ux", "wrong uxmod");
54-
}
55-
for (int i = 0; i < 3; i++) {
56-
ucell.magnet.ux_[i] *= 1 / sqrt(uxmod);
57-
}
58-
// std::cout<<" Fixed quantization axis for GGA: "
59-
//<<std::setw(10)<<ux[0]<<" "<<std::setw(10)<<ux[1]<<"
60-
//"<<std::setw(10)<<ux[2]<<std::endl;
61-
}
62-
return;
59+
for (int it = starting_it; it < ucell.ntype; it++)
60+
{
61+
for (int ia = 0; ia < ucell.atoms[it].na; ia++)
62+
{
63+
if (it > starting_it || ia > starting_ia)
64+
{
65+
ucell.magnet.lsign_
66+
= ucell.magnet.lsign_
67+
&& judge_parallel(ucell.magnet.ux_, ucell.atoms[it].m_loc_[ia]);
68+
}
69+
}
70+
}
71+
72+
// if all of the atoms have the same parallel magnetism direction,
73+
// then set the direction to a unit vector
74+
if (ucell.magnet.lsign_)
75+
{
76+
uxmod = pow(ucell.magnet.ux_[0], 2)
77+
+ pow(ucell.magnet.ux_[1], 2)
78+
+ pow(ucell.magnet.ux_[2], 2);
79+
80+
if (uxmod < absolute_mag_thr)
81+
{
82+
ModuleBase::WARNING_QUIT("elecstate::cal_ux", "wrong uxmod");
83+
}
84+
85+
// reset the magnetism for each direction
86+
for (int i = 0; i < 3; i++)
87+
{
88+
ucell.magnet.ux_[i] *= 1 / sqrt(uxmod);
89+
}
90+
}
91+
return;
6392
}
6493

6594
bool judge_parallel(double a[3], ModuleBase::Vector3<double> b) {
@@ -71,4 +100,4 @@ bool judge_parallel(double a[3], ModuleBase::Vector3<double> b) {
71100
jp = (fabs(cross) < 1e-6);
72101
return jp;
73102
}
74-
}
103+
}

‎source/module_elecstate/magnetism.cpp‎

Lines changed: 28 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -22,6 +22,8 @@ void Magnetism::compute_magnetization(const double& omega,
2222
const double* const * rho,
2323
double* nelec_spin)
2424
{
25+
assert(nxyz>0);
26+
2527
if (PARAM.inp.nspin==2)
2628
{
2729
this->tot_magnetization = 0.00;
@@ -57,36 +59,53 @@ void Magnetism::compute_magnetization(const double& omega,
5759
// noncolliear :
5860
else if(PARAM.inp.nspin==4)
5961
{
60-
for(int i=0;i<3;i++) {this->tot_magnetization_nc[i] = 0.00;
61-
}
62+
for(int i=0;i<3;i++)
63+
{
64+
this->tot_magnetization_nc[i] = 0.00;
65+
}
66+
6267
this->abs_magnetization = 0.00;
6368
for (int ir=0; ir<nrxx; ir++)
6469
{
6570
double diff = sqrt(pow(rho[1][ir], 2) + pow(rho[2][ir], 2) +pow(rho[3][ir], 2));
6671

67-
for(int i=0;i<3;i++) {this->tot_magnetization_nc[i] += rho[i+1][ir];
68-
}
72+
for(int i=0;i<3;i++)
73+
{
74+
this->tot_magnetization_nc[i] += rho[i+1][ir];
75+
}
6976
this->abs_magnetization += std::abs(diff);
7077
}
7178
#ifdef __MPI
7279
Parallel_Reduce::reduce_pool(this->tot_magnetization_nc, 3);
7380
Parallel_Reduce::reduce_pool(this->abs_magnetization);
7481
#endif
75-
for(int i=0;i<3;i++) {this->tot_magnetization_nc[i] *= omega/ nxyz;
76-
}
82+
for(int i=0;i<3;i++)
83+
{
84+
this->tot_magnetization_nc[i] *= omega/ nxyz;
85+
}
86+
7787
this->abs_magnetization *= omega/ nxyz;
78-
GlobalV::ofs_running<<"total magnetism (Bohr mag/cell)"<<'\t'<<this->tot_magnetization_nc[0]<<'\t'<<this->tot_magnetization_nc[1]<<'\t'<<this->tot_magnetization_nc[2]<<'\n';
88+
89+
GlobalV::ofs_running<<"total magnetism (Bohr mag/cell)"
90+
<<'\t'<<this->tot_magnetization_nc[0]
91+
<<'\t'<<this->tot_magnetization_nc[1]
92+
<<'\t'<<this->tot_magnetization_nc[2]<<'\n';
93+
7994
ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"absolute magnetism (Bohr mag/cell)",this->abs_magnetization);
8095
}
8196

8297
return;
8398
}
8499

100+
85101
bool Magnetism::judge_parallel(double a[3], ModuleBase::Vector3<double> b)
86102
{
87103
bool jp=false;
88-
double cross;
89-
cross = pow((a[1]*b.z-a[2]*b.y),2) + pow((a[2]*b.x-a[0]*b.z),2) + pow((a[0]*b.y-a[1]*b.x),2);
104+
double cross=0.0;
105+
cross = pow((a[1]*b.z-a[2]*b.y),2)
106+
+ pow((a[2]*b.x-a[0]*b.z),2)
107+
+ pow((a[0]*b.y-a[1]*b.x),2);
108+
90109
jp = (fabs(cross)<1e-6);
91110
return jp;
92111
}

‎source/module_elecstate/magnetism.h‎

Lines changed: 4 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -13,7 +13,7 @@ class Magnetism
1313
Magnetism();
1414
~Magnetism();
1515

16-
// notice : becast is done in unitcell
16+
// notice : bcast (MPI operation) is done in unitcell
1717
double *start_magnetization=nullptr;
1818

1919
// tot_magnetization : majority spin - minority spin (nelup - neldw).
@@ -28,10 +28,9 @@ class Magnetism
2828

2929
double* nelec_spin = nullptr);
3030

31-
ModuleBase::Vector3<double> *m_loc_=nullptr;//magnetization for each element along c-axis
32-
double *angle1_=nullptr; //angle between c-axis and real spin std::vector
33-
double *angle2_=nullptr; //angle between a-axis and real spin std::vector projection in ab-plane
34-
//void cal_ux(const int ntype);
31+
ModuleBase::Vector3<double> *m_loc_=nullptr; //magnetization for each element along c-axis
32+
double *angle1_=nullptr; //angle between c-axis and real spin std::vector
33+
double *angle2_=nullptr; //angle between a-axis and real spin std::vector projection in ab-plane
3534
double ux_[3]={0.0};
3635
bool lsign_=false;
3736

‎source/module_esolver/esolver_fp.cpp‎

Lines changed: 18 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -258,6 +258,7 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep)
258258
{
259259
ModuleBase::TITLE("ESolver_FP", "before_scf");
260260

261+
// if the cell has changed
261262
if (ucell.cell_parameter_updated)
262263
{
263264
// only G-vector and K-vector are changed due to the change of lattice
@@ -266,62 +267,71 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep)
266267
this->pw_rho->collect_local_pw();
267268
this->pw_rho->collect_uniqgg();
268269

270+
// if double grid used in USPP, update related quantities in dense grid
269271
if (PARAM.globalv.double_grid)
270272
{
271273
this->pw_rhod->initgrids(ucell.lat0, ucell.latvec, pw_rhod->nx, pw_rhod->ny, pw_rhod->nz);
272274
this->pw_rhod->collect_local_pw();
273275
this->pw_rhod->collect_uniqgg();
274276
}
275277

278+
// reset local pseudopotentials
276279
this->locpp.init_vloc(ucell, this->pw_rhod);
277280
ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "LOCAL POTENTIAL");
278281

279282
this->pelec->omega = ucell.omega;
280283

284+
// perform symmetry analysis
281285
if (ModuleSymmetry::Symmetry::symm_flag == 1)
282286
{
283287
ucell.symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, GlobalV::ofs_running);
284288
ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SYMMETRY");
285289
}
286290

291+
// reset k-points
287292
kv.set_after_vc(PARAM.inp.nspin, ucell.G, ucell.latvec);
288293
ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT K-POINTS");
289294
}
290295

296+
//----------------------------------------------------------
291297
// charge extrapolation
298+
//----------------------------------------------------------
292299
if (ucell.ionic_position_updated)
293300
{
294301
this->CE.update_all_dis(ucell);
295-
this->CE.extrapolate_charge(&(this->Pgrid),
302+
this->CE.extrapolate_charge(&this->Pgrid,
296303
ucell,
297304
&this->chr,
298-
&(this->sf),
305+
&this->sf,
299306
GlobalV::ofs_running,
300307
GlobalV::ofs_warning);
301308
}
302309

303310
//----------------------------------------------------------
304-
// about vdw, jiyy add vdwd3 and linpz add vdwd2
311+
//! calculate D2 or D3 vdW
305312
//----------------------------------------------------------
306313
auto vdw_solver = vdw::make_vdw(ucell, PARAM.inp, &(GlobalV::ofs_running));
307314
if (vdw_solver != nullptr)
308315
{
309316
this->pelec->f_en.evdw = vdw_solver->get_energy();
310317
}
311318

312-
// calculate ewald energy
319+
//----------------------------------------------------------
320+
//! calculate ewald energy
321+
//----------------------------------------------------------
313322
if (!PARAM.inp.test_skip_ewald)
314323
{
315324
this->pelec->f_en.ewald_energy = H_Ewald_pw::compute_ewald(ucell, this->pw_rhod, this->sf.strucFac);
316325
}
317326

318327
//----------------------------------------------------------
319-
//! cal_ux should be called before init_scf because
320-
//! the direction of ux is used in noncoline_rho
328+
//! set direction of magnetism, used in non-collinear case
321329
//----------------------------------------------------------
322330
elecstate::cal_ux(ucell);
323331

332+
//----------------------------------------------------------
324333
//! output the initial charge density
334+
//----------------------------------------------------------
325335
if (PARAM.inp.out_chg[0] == 2)
326336
{
327337
for (int is = 0; is < PARAM.inp.nspin; is++)
@@ -339,7 +349,9 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep)
339349
}
340350
}
341351

352+
//----------------------------------------------------------
342353
//! output total local potential of the initial charge density
354+
//----------------------------------------------------------
343355
if (PARAM.inp.out_pot == 3)
344356
{
345357
for (int is = 0; is < PARAM.inp.nspin; is++)

‎source/module_esolver/esolver_ks_lcao.cpp‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -160,7 +160,9 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_pa
160160
#ifdef __EXX
161161
// 5) initialize exx
162162
// PLEASE simplify the Exx_Global interface
163-
if (PARAM.inp.calculation == "scf" || PARAM.inp.calculation == "relax" || PARAM.inp.calculation == "cell-relax"
163+
if (PARAM.inp.calculation == "scf"
164+
|| PARAM.inp.calculation == "relax"
165+
|| PARAM.inp.calculation == "cell-relax"
164166
|| PARAM.inp.calculation == "md")
165167
{
166168
if (GlobalC::exx_info.info_global.cal_exx)

0 commit comments

Comments
 (0)