Skip to content

Commit f2a485a

Browse files
YuLiu98YuLiu98
authored
Fix: fix bug in relax_old (#2150)
1 parent f4d8eda commit f2a485a

File tree

12 files changed

+307
-385
lines changed

12 files changed

+307
-385
lines changed

‎examples/relax/pw_al/log.ref‎

Lines changed: 190 additions & 320 deletions
Large diffs are not rendered by default.

‎source/module_cell/unitcell.cpp‎

Lines changed: 38 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -293,6 +293,34 @@ void UnitCell::set_iat2itia(void)
293293
return;
294294
}
295295

296+
void UnitCell::update_pos_tau(const double* pos)
297+
{
298+
int iat = 0;
299+
for (int it = 0; it < this->ntype; it++)
300+
{
301+
Atom* atom = &this->atoms[it];
302+
for (int ia = 0; ia < atom->na; ia++)
303+
{
304+
for ( int ik = 0; ik < 3; ++ik)
305+
{
306+
if (atom->mbl[ia][ik])
307+
{
308+
atom->dis[ia][ik] = pos[3 * iat + ik] / this->lat0 - atom->tau[ia][ik];
309+
atom->tau[ia][ik] = pos[3 * iat + ik] / this->lat0;
310+
}
311+
}
312+
313+
// the direct coordinates also need to be updated.
314+
atom->dis[ia] = atom->dis[ia] * this->GT;
315+
atom->taud[ia] = atom->tau[ia] * this->GT;
316+
iat++;
317+
}
318+
}
319+
assert(iat == this->nat);
320+
this->periodic_boundary_adjustment();
321+
this->bcast_atoms_tau();
322+
}
323+
296324
void UnitCell::update_pos_taud(double* posd_in)
297325
{
298326
int iat = 0;
@@ -301,9 +329,11 @@ void UnitCell::update_pos_taud(double* posd_in)
301329
Atom* atom = &this->atoms[it];
302330
for (int ia = 0; ia < atom->na; ia++)
303331
{
304-
this->atoms[it].taud[ia].x += posd_in[iat*3];
305-
this->atoms[it].taud[ia].y += posd_in[iat*3 + 1];
306-
this->atoms[it].taud[ia].z += posd_in[iat*3 + 2];
332+
for ( int ik = 0; ik < 3; ++ik)
333+
{
334+
atom->taud[ia][ik] += posd_in[3*iat + ik];
335+
atom->dis[ia][ik] = posd_in[3*iat + ik];
336+
}
307337
iat++;
308338
}
309339
}
@@ -315,18 +345,21 @@ void UnitCell::update_pos_taud(double* posd_in)
315345
// posd_in is atomic displacements here liuyu 2023-03-22
316346
void UnitCell::update_pos_taud(const ModuleBase::Vector3<double>* posd_in)
317347
{
348+
int iat = 0;
318349
for (int it = 0; it < this->ntype; it++)
319350
{
320351
Atom* atom = &this->atoms[it];
321352
for (int ia = 0; ia < atom->na; ia++)
322353
{
323354
for ( int ik = 0; ik < 3; ++ik)
324355
{
325-
atom->taud[ia][ik] += posd_in[ia][ik];
326-
atom->dis[ia][ik] = posd_in[ia][ik];
356+
atom->taud[ia][ik] += posd_in[iat][ik];
357+
atom->dis[ia][ik] = posd_in[iat][ik];
327358
}
359+
iat++;
328360
}
329361
}
362+
assert(iat == this->nat);
330363
this->periodic_boundary_adjustment();
331364
this->bcast_atoms_tau();
332365
}

‎source/module_cell/unitcell.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -180,6 +180,7 @@ class UnitCell
180180
void print_cell_xyz(const std::string &fn)const;
181181
void print_cell_cif(const std::string &fn)const;
182182

183+
void update_pos_tau(const double* pos);
183184
void update_pos_taud(const ModuleBase::Vector3<double>* posd_in);
184185
void update_pos_taud(double* posd_in);
185186
void update_vel(const ModuleBase::Vector3<double>* vel_in);

‎source/module_relax/relax_old/bfgs_basic.cpp‎

Lines changed: 15 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -8,6 +8,8 @@ double BFGS_Basic::relax_bfgs_w2 = -1.0; // defalut is 0.05
88

99
BFGS_Basic::BFGS_Basic()
1010
{
11+
pos = nullptr;
12+
pos_p = nullptr;
1113
grad = nullptr;
1214
grad_p = nullptr;
1315
move = nullptr;
@@ -18,6 +20,8 @@ BFGS_Basic::BFGS_Basic()
1820

1921
BFGS_Basic::~BFGS_Basic()
2022
{
23+
delete[] pos;
24+
delete[] pos_p;
2125
delete[] grad;
2226
delete[] grad_p;
2327
delete[] move;
@@ -28,17 +32,23 @@ void BFGS_Basic::allocate_basic(void)
2832
{
2933
assert(dim>0);
3034

35+
delete[] pos;
36+
delete[] pos_p;
3137
delete[] grad;
3238
delete[] grad_p;
3339
delete[] move;
3440
delete[] move_p;
3541

42+
pos = new double[dim];
43+
pos_p = new double [dim];
3644
grad = new double[dim];
3745
grad_p = new double [dim];
3846
move = new double [dim];
3947
move_p = new double [dim];
4048

49+
ModuleBase::GlobalFunc::ZEROS(pos, dim);
4150
ModuleBase::GlobalFunc::ZEROS(grad, dim);
51+
ModuleBase::GlobalFunc::ZEROS(pos_p, dim);
4252
ModuleBase::GlobalFunc::ZEROS(grad_p, dim);
4353
ModuleBase::GlobalFunc::ZEROS(move, dim);
4454
ModuleBase::GlobalFunc::ZEROS(move_p, dim);
@@ -61,8 +71,9 @@ void BFGS_Basic::update_inverse_hessian(void)
6171

6272
for(int i=0;i<dim;i++)
6373
{
74+
// s[i] = this->pos[i] - this->pos_p[i];
6475
// mohan update 2010-07-27
65-
s[i] = this->check_move( move_p[i] );
76+
s[i] = this->check_move( pos[i], pos_p[i] );
6677
s[i] *= GlobalC::ucell.lat0;
6778

6879
y[i] = this->grad[i] - this->grad_p[i];
@@ -190,6 +201,7 @@ void BFGS_Basic::save_bfgs(void)
190201
this->save_flag = true;
191202
for(int i=0;i<dim;i++)
192203
{
204+
this->pos_p[i] = this->pos[i];
193205
this->grad_p[i] = this->grad[i];
194206
this->move_p[i] = this->move[i];
195207
}
@@ -378,12 +390,12 @@ void BFGS_Basic::compute_trust_radius(void)
378390
return;
379391
}
380392

381-
double BFGS_Basic::check_move(const double &move_p)
393+
double BFGS_Basic::check_move(const double &pos, const double &pos_p)
382394
{
383395
// this must be careful.
384396
// unit is GlobalC::ucell.lat0.
385397
assert(GlobalC::ucell.lat0>0.0);
386-
const double direct_move = move_p / GlobalC::ucell.lat0;
398+
const double direct_move = (pos - pos_p)/GlobalC::ucell.lat0;
387399
double shortest_move = direct_move;
388400
for(int cell=-1; cell<=1; ++cell)
389401
{

‎source/module_relax/relax_old/bfgs_basic.h‎

Lines changed: 5 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -29,9 +29,11 @@ class BFGS_Basic
2929
void reset_hessian(void);
3030
void save_bfgs(void);
3131

32-
double* grad; // std::vector containing 3N components of ( grad( V(x) ) )
33-
double* move; // movements
32+
double* pos; // std::vector containing 3N coordinates of the system ( x )
33+
double* grad; //std::vector containing 3N components of ( grad( V(x) ) )
34+
double* move; // pos = pos_p + move.
3435

36+
double* pos_p; // p: previous
3537
double* grad_p; // p: previous
3638
double* move_p;
3739

@@ -47,7 +49,7 @@ class BFGS_Basic
4749
// to the minimum value at the previous step
4850

4951
// mohan add 2010-07-27
50-
double check_move(const double &move_p);
52+
double check_move(const double &pos, const double &pos_p);
5153
private:
5254

5355
bool wolfe_flag;

‎source/module_relax/relax_old/ions_move_basic.cpp‎

Lines changed: 25 additions & 43 deletions
Original file line numberDiff line numberDiff line change
@@ -20,51 +20,49 @@ double Ions_Move_Basic::best_xxx = 1.0;
2020

2121
int Ions_Move_Basic::out_stru=0;
2222

23-
void Ions_Move_Basic::setup_gradient(double *grad, const ModuleBase::matrix &force)
23+
void Ions_Move_Basic::setup_gradient(double* pos, double *grad, const ModuleBase::matrix &force)
2424
{
2525
ModuleBase::TITLE("Ions_Move_Basic","setup_gradient");
2626

2727
assert(GlobalC::ucell.ntype>0);
28+
assert(pos!=NULL);
2829
assert(grad!=NULL);
2930
assert(dim == 3*GlobalC::ucell.nat);
3031

32+
ModuleBase::GlobalFunc::ZEROS(pos, dim);
3133
ModuleBase::GlobalFunc::ZEROS(grad, dim);
3234

3335
// (1) init gradient
36+
// the unit of pos: Bohr.
3437
// the unit of force: Ry/Bohr.
3538
// the unit of gradient:
3639
int iat=0;
3740
for(int it = 0;it < GlobalC::ucell.ntype;it++)
3841
{
3942
Atom* atom = &GlobalC::ucell.atoms[it];
4043
for(int ia =0;ia< GlobalC::ucell.atoms[it].na;ia++)
41-
{
42-
if(atom->mbl[ia].x == 1)
43-
{
44-
grad[3*iat ] = -force(iat, 0)*GlobalC::ucell.lat0;
45-
//this->grad[3*iat ] = -force(iat, 0);
46-
}
47-
if(atom->mbl[ia].y == 1)
48-
{
49-
grad[3*iat+1] = -force(iat, 1)*GlobalC::ucell.lat0;
50-
}
51-
if(atom->mbl[ia].z == 1)
52-
{
53-
grad[3*iat+2] = -force(iat, 2)*GlobalC::ucell.lat0;
54-
//std::cout << " grad=" << grad[3*iat+2] << std::endl;
55-
}
44+
{
45+
for ( int ik = 0; ik < 3; ++ik)
46+
{
47+
pos[3*iat + ik] = atom->tau[ia][ik] * GlobalC::ucell.lat0;
48+
if (atom->mbl[ia][ik])
49+
{
50+
grad[3*iat + ik] = - force(iat, ik) * GlobalC::ucell.lat0;
51+
}
52+
}
5653
++iat;
5754
}
5855
}
5956

6057
return;
6158
}
6259

63-
void Ions_Move_Basic::move_atoms(double *move)
60+
void Ions_Move_Basic::move_atoms(double *move, double *pos)
6461
{
6562
ModuleBase::TITLE("Ions_Move_Basic","move_atoms");
6663

6764
assert(move!=NULL);
65+
assert(pos!=NULL);
6866

6967
//------------------------
7068
// for test only
@@ -91,31 +89,16 @@ void Ions_Move_Basic::move_atoms(double *move)
9189
assert( iat == GlobalC::ucell.nat );
9290
}
9391

94-
int iat = 0;
95-
const double move_threshold = 1.0e-10;
96-
ModuleBase::Vector3<double> *move_ion = new ModuleBase::Vector3<double> [GlobalC::ucell.nat];
97-
for(int it = 0; it < GlobalC::ucell.ntype; it++)
98-
{
99-
Atom* atom = &GlobalC::ucell.atoms[it];
100-
for(int ia = 0; ia < atom->na; ia++)
101-
{
102-
for (int k = 0; k < 3; ++k)
103-
{
104-
if( abs(move[3*iat + k]) > move_threshold && atom->mbl[ia][k])
105-
{
106-
move_ion[iat][k] = move[3*iat + k] / GlobalC::ucell.lat0;
107-
}
108-
else
109-
{
110-
move_ion[iat][k] = 0;
111-
}
112-
}
113-
move_ion[iat] = move_ion[iat] * GlobalC::ucell.GT;
114-
iat++;
115-
}
116-
}
117-
assert( iat == GlobalC::ucell.nat );
118-
GlobalC::ucell.update_pos_taud(move_ion);
92+
const double move_threshold = 1.0e-10;
93+
const int total_freedom = GlobalC::ucell.nat * 3;
94+
for(int i =0;i<total_freedom;i++)
95+
{
96+
if( abs(move[i]) > move_threshold )
97+
{
98+
pos[i] += move[i];
99+
}
100+
}
101+
GlobalC::ucell.update_pos_tau(pos);
119102

120103
//--------------------------------------------
121104
// Print out the structure file.
@@ -128,7 +111,6 @@ void Ions_Move_Basic::move_atoms(double *move)
128111
{
129112
GlobalC::ucell.print_cell_cif("STRU_NOW.cif");
130113
}
131-
delete[] move_ion;
132114
return;
133115
}
134116

‎source/module_relax/relax_old/ions_move_basic.h‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -26,12 +26,12 @@ namespace Ions_Move_Basic
2626
//----------------------------------------------------------------------------
2727
// setup the gradient, all the same for any geometry optimization methods.
2828
//----------------------------------------------------------------------------
29-
void setup_gradient(double* grad, const ModuleBase::matrix &force);
29+
void setup_gradient(double *pos, double* grad, const ModuleBase::matrix &force);
3030

3131
//----------------------------------------------------------------------------
3232
// move the atom positions, considering the periodic boundary condition.
3333
//----------------------------------------------------------------------------
34-
void move_atoms(double* move);
34+
void move_atoms(double* move, double *pos);
3535

3636
//----------------------------------------------------------------------------
3737
// check the converged conditions ( if largest gradient is smaller than

‎source/module_relax/relax_old/ions_move_bfgs.cpp‎

Lines changed: 5 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -40,7 +40,7 @@ void Ions_Move_BFGS::start(const ModuleBase::matrix& force, const double &energy
4040
// istep must be set eariler.
4141

4242
// use force to setup gradient.
43-
Ions_Move_Basic::setup_gradient(this->grad, force);
43+
Ions_Move_Basic::setup_gradient(this->pos, this->grad, force);
4444
// use energy_in and istep to setup etot and etot_old.
4545
Ions_Move_Basic::setup_etot(energy_in, 0);
4646
// use gradient and etot and etot_old to check
@@ -72,7 +72,7 @@ void Ions_Move_BFGS::start(const ModuleBase::matrix& force, const double &energy
7272
// even if the energy is higher, we save the information.
7373
this->save_bfgs();
7474

75-
Ions_Move_Basic::move_atoms(move);
75+
Ions_Move_Basic::move_atoms(move, pos);
7676
}
7777
return;
7878
}
@@ -109,13 +109,14 @@ void Ions_Move_BFGS::restart_bfgs(void)
109109
{
110110
// mohan add 2010-07-26.
111111
// there must be one of the two has the correct sign and value.
112-
this->move_p[i] = this->check_move(move_p[i])/trust_radius_old;
112+
this->move_p[i] = this->check_move(pos[i], pos_p[i])/trust_radius_old;
113113
//std::cout << " " << std::setw(20) << move_p[i] << std::setw(20) << dpmin << std::endl;
114114
}
115115
}
116116
else
117117
{
118118
// bfgs initialization
119+
ModuleBase::GlobalFunc::ZEROS(pos_p, dim);
119120
ModuleBase::GlobalFunc::ZEROS(grad_p, dim);
120121
ModuleBase::GlobalFunc::ZEROS(move_p, dim);
121122

@@ -233,6 +234,7 @@ void Ions_Move_BFGS::bfgs_routine(void)
233234
etot = etot_p;
234235
for(int i=0;i<dim;i++)
235236
{
237+
this->pos[i] = pos_p[i];
236238
this->grad[i] = grad_p[i];
237239
}
238240

0 commit comments

Comments
 (0)