Skip to content

Commit b05a96b

Browse files
dyzhengdyzheng
authored
Merge pull request #1169 from wenfei-li/develop
fix bug in 2nd order charge extrapolation
2 parents 844eeb3 + 1cd7955 commit b05a96b

File tree

8 files changed

+53
-4
lines changed

8 files changed

+53
-4
lines changed

‎source/module_cell/atom_spec.cpp‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -12,6 +12,7 @@ Atom::Atom()
1212
type = 0;
1313
stapos_wf = 0;
1414
tau = new ModuleBase::Vector3<double>[1];
15+
tau_original = new ModuleBase::Vector3<double>[1];
1516
taud = new ModuleBase::Vector3<double>[1];
1617
vel = new ModuleBase::Vector3<double>[1];
1718
mag = new double[1];
@@ -29,6 +30,7 @@ Atom::Atom()
2930
Atom::~Atom()
3031
{
3132
delete[] tau;
33+
delete[] tau_original;
3234
delete[] taud;
3335
delete[] vel;
3436
delete[] mag;
@@ -147,6 +149,7 @@ void Atom::bcast_atom(void)
147149
delete[] angle2;
148150
delete[] m_loc_;
149151
tau = new ModuleBase::Vector3<double>[na];
152+
tau_original = new ModuleBase::Vector3<double>[na];
150153
taud = new ModuleBase::Vector3<double>[na];
151154
vel = new ModuleBase::Vector3<double>[na];
152155
mag = new double[na];

‎source/module_cell/atom_spec.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -31,6 +31,7 @@ class Atom: public Atom_pseudo
3131

3232
std::string label; // atomic symbol
3333
ModuleBase::Vector3<double> *tau;// Cartesian coordinates of each atom in this type.
34+
ModuleBase::Vector3<double> *tau_original;// Cartesian coordinates of each atom in this type, but without periodic adjustment.
3435
ModuleBase::Vector3<double> *taud;// Direct coordinates of each atom in this type.
3536
ModuleBase::Vector3<double> *vel;// velocities of each atom in this type.
3637

‎source/module_cell/read_atoms.cpp‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -575,6 +575,7 @@ bool UnitCell_pseudo::read_atom_positions(std::ifstream &ifpos, std::ofstream &o
575575
delete[] atoms[it].angle2;
576576
delete[] atoms[it].m_loc_;
577577
atoms[it].tau = new ModuleBase::Vector3<double>[na];
578+
atoms[it].tau_original = new ModuleBase::Vector3<double>[na];
578579
atoms[it].taud = new ModuleBase::Vector3<double>[na];
579580
atoms[it].vel = new ModuleBase::Vector3<double>[na];
580581
atoms[it].mbl = new ModuleBase::Vector3<int>[na];
@@ -800,6 +801,7 @@ bool UnitCell_pseudo::read_atom_positions(std::ifstream &ifpos, std::ofstream &o
800801
}
801802

802803
atoms[it].mbl[ia] = mv;
804+
atoms[it].tau_original[ia] = atoms[it].tau[ia];
803805
}//endj
804806
}// end na
805807
}//end for ntype

‎source/module_cell/unitcell.cpp‎

Lines changed: 41 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -283,14 +283,17 @@ void UnitCell::update_pos_tau(const double* pos)
283283
{
284284
if (atom->mbl[ia].x != 0)
285285
{
286+
atom->tau_original[ia].x += (pos[3 * iat] / this->lat0 - atom->tau[ia].x);
286287
atom->tau[ia].x = pos[3 * iat] / this->lat0;
287288
}
288289
if (atom->mbl[ia].y != 0)
289290
{
291+
atom->tau_original[ia].y += (pos[3 * iat + 1] / this->lat0 - atom->tau[ia].y);
290292
atom->tau[ia].y = pos[3 * iat + 1] / this->lat0;
291293
}
292294
if (atom->mbl[ia].z != 0)
293295
{
296+
atom->tau_original[ia].z += (pos[3 * iat + 2] / this->lat0 - atom->tau[ia].z);
294297
atom->tau[ia].z = pos[3 * iat + 2] / this->lat0;
295298
}
296299

@@ -313,14 +316,17 @@ void UnitCell::update_pos_tau(const ModuleBase::Vector3<double>* posd_in)
313316
{
314317
if (atom->mbl[ia].x != 0)
315318
{
319+
atom->tau_original[ia].x += (posd_in[iat].x / this->lat0 - atom->tau[ia].x);
316320
atom->tau[ia].x = posd_in[iat].x / this->lat0;
317321
}
318322
if (atom->mbl[ia].y != 0)
319323
{
324+
atom->tau_original[ia].y += (posd_in[iat].y / this->lat0 - atom->tau[ia].y);
320325
atom->tau[ia].y = posd_in[iat].y / this->lat0;
321326
}
322327
if (atom->mbl[ia].z != 0)
323328
{
329+
atom->tau_original[ia].z += (posd_in[iat].z / this->lat0 - atom->tau[ia].z);
324330
atom->tau[ia].z = posd_in[iat].z / this->lat0;
325331
}
326332

@@ -444,7 +450,41 @@ void UnitCell::save_cartesian_position(ModuleBase::Vector3<double>* pos) const
444450
Atom* atom = &this->atoms[it];
445451
for (int ia = 0; ia < atom->na; ++ia)
446452
{
447-
pos[iat] = atom->tau[ia] * this->lat0;
453+
pos[iat] = atom->tau_original[ia] * this->lat0;
454+
iat++;
455+
}
456+
}
457+
assert(iat == this->nat);
458+
return;
459+
}
460+
461+
void UnitCell::save_cartesian_position_original(double* pos) const
462+
{
463+
int iat = 0;
464+
for (int it = 0; it < this->ntype; it++)
465+
{
466+
Atom* atom = &this->atoms[it];
467+
for (int ia = 0; ia < atom->na; ia++)
468+
{
469+
pos[3 * iat] = atom->tau_original[ia].x * this->lat0;
470+
pos[3 * iat + 1] = atom->tau_original[ia].y * this->lat0;
471+
pos[3 * iat + 2] = atom->tau_original[ia].z * this->lat0;
472+
iat++;
473+
}
474+
}
475+
assert(iat == this->nat);
476+
return;
477+
}
478+
479+
void UnitCell::save_cartesian_position_original(ModuleBase::Vector3<double>* pos) const
480+
{
481+
int iat = 0;
482+
for (int it = 0; it < this->ntype; ++it)
483+
{
484+
Atom* atom = &this->atoms[it];
485+
for (int ia = 0; ia < atom->na; ++ia)
486+
{
487+
pos[iat] = atom->tau_original[ia] * this->lat0;
448488
iat++;
449489
}
450490
}

‎source/module_cell/unitcell.h‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -89,7 +89,8 @@ class UnitCell
8989
void bcast_atoms_tau();
9090
void save_cartesian_position(double* pos)const;
9191
void save_cartesian_position(ModuleBase::Vector3<double>* pos)const;
92-
92+
void save_cartesian_position_original(double* pos)const;
93+
void save_cartesian_position_original(ModuleBase::Vector3<double>* pos)const;
9394
bool judge_big_cell(void)const;
9495

9596

‎source/module_md/test/setcell.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -56,6 +56,7 @@ class Setcell
5656
delete[] ucell.atoms[0].vel;
5757
delete[] ucell.atoms[0].mbl;
5858
ucell.atoms[0].tau = new ModuleBase::Vector3<double>[4];
59+
ucell.atoms[0].tau_original = new ModuleBase::Vector3<double>[4];
5960
ucell.atoms[0].taud = new ModuleBase::Vector3<double>[4];
6061
ucell.atoms[0].vel = new ModuleBase::Vector3<double>[4];
6162
ucell.atoms[0].mbl = new ModuleBase::Vector3<int>[4];

‎source/src_lcao/test/gamma_rho_mock.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -597,6 +597,7 @@ bool UnitCell_pseudo::read_atom_positions(LCAO_Orbitals &orb,
597597
if (na > 0)
598598
{
599599
atoms[it].tau = new ModuleBase::Vector3<double>[na];
600+
atoms[it].tau_original = new ModuleBase::Vector3<double>[na];
600601
atoms[it].taud = new ModuleBase::Vector3<double>[na];
601602
atoms[it].vel = new ModuleBase::Vector3<double>[na];
602603
for (int ia = 0; ia < na; ia++)

‎source/src_pw/charge_extra.cpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -491,7 +491,7 @@ void Charge_Extra::find_alpha_and_beta(void)
491491

492492
void Charge_Extra::save_pos_next(const UnitCell_pseudo& ucell)
493493
{
494-
ucell.save_cartesian_position(this->pos_next);
494+
ucell.save_cartesian_position_original(this->pos_next);
495495
return;
496496
}
497497

@@ -508,6 +508,6 @@ void Charge_Extra::update_all_pos(const UnitCell_pseudo& ucell)
508508
this->pos_old2[i] = this->pos_old1[i];
509509
this->pos_old1[i] = this->pos_now[i];
510510
}
511-
ucell.save_cartesian_position(this->pos_now);
511+
ucell.save_cartesian_position_original(this->pos_now);
512512
return;
513513
}

0 commit comments

Comments
 (0)