Skip to content

Commit aae122a

Browse files
mohanchenmohanchen
committed
change update_cell.cpp algorithm, when the ion move is larger than the cell length, it is fine to proceed the relaxation calculations
1 parent 8073403 commit aae122a

File tree

5 files changed

+108
-76
lines changed

5 files changed

+108
-76
lines changed

‎source/source_cell/atom_spec.h‎

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -38,8 +38,7 @@ class Atom
3838
std::vector<ModuleBase::Vector3<double>> taud; // Direct coordinates of each atom in this type.
3939
std::vector<ModuleBase::Vector3<double>> vel; // velocities of each atom in this type.
4040
std::vector<ModuleBase::Vector3<double>> force; // force acting on each atom in this type.
41-
std::vector<ModuleBase::Vector3<double>>
42-
lambda; // Lagrange multiplier for each atom in this type. used in deltaspin
41+
std::vector<ModuleBase::Vector3<double>> lambda; // Lagrange multiplier for each atom in this type. used in deltaspin
4342
std::vector<ModuleBase::Vector3<int>> constrain; // constrain for each atom in this type. used in deltaspin
4443
std::string label_orb = "\0"; // atomic Element symbol in the orbital file of lcao
4544

‎source/source_cell/update_cell.cpp‎

Lines changed: 90 additions & 59 deletions
Original file line numberDiff line numberDiff line change
@@ -1,8 +1,10 @@
11
#include "update_cell.h"
22
#include "bcast_cell.h"
33
#include "source_base/global_function.h"
4+
45
namespace unitcell
56
{
7+
68
void remake_cell(Lattice& lat)
79
{
810
ModuleBase::TITLE("Lattice", "rmake_cell");
@@ -13,18 +15,20 @@ void remake_cell(Lattice& lat)
1315
std::string& latName = lat.latName;
1416
ModuleBase::Matrix3& latvec = lat.latvec;
1517

16-
if (latName == "none") {
17-
ModuleBase::WARNING_QUIT("UnitCell",
18-
"to use fixed_ibrav, latname must be provided");
19-
} else if (latName == "sc") // ibrav = 1
18+
if (latName == "none")
2019
{
20+
ModuleBase::WARNING_QUIT("UnitCell", "to use fixed_ibrav, latname must be provided");
21+
}
22+
else if (latName == "sc") // ibrav = 1
23+
{
2124
double celldm = std::sqrt(pow(latvec.e11, 2) + pow(latvec.e12, 2)
2225
+ pow(latvec.e13, 2));
2326

2427
latvec.Zero();
2528
latvec.e11 = latvec.e22 = latvec.e33 = celldm;
26-
} else if (latName == "fcc") // ibrav = 2
27-
{
29+
}
30+
else if (latName == "fcc") // ibrav = 2
31+
{
2832
double celldm = std::sqrt(pow(latvec.e11, 2) + pow(latvec.e12, 2)
2933
+ pow(latvec.e13, 2)) / std::sqrt(2.0);
3034

@@ -37,8 +41,9 @@ void remake_cell(Lattice& lat)
3741
latvec.e31 = -celldm;
3842
latvec.e32 = celldm;
3943
latvec.e33 = 0.0;
40-
} else if (latName == "bcc") // ibrav = 3
41-
{
44+
}
45+
else if (latName == "bcc") // ibrav = 3
46+
{
4247
double celldm = std::sqrt(pow(latvec.e11, 2) + pow(latvec.e12, 2)
4348
+ pow(latvec.e13, 2))
4449
/ std::sqrt(3.0);
@@ -52,9 +57,10 @@ void remake_cell(Lattice& lat)
5257
latvec.e31 = -celldm;
5358
latvec.e32 = -celldm;
5459
latvec.e33 = celldm;
55-
} else if (latName == "hexagonal") // ibrav = 4
56-
{
57-
double celldm1 = std::sqrt(pow(latvec.e11, 2) + pow(latvec.e12, 2)
60+
}
61+
else if (latName == "hexagonal") // ibrav = 4
62+
{
63+
double celldm1 = std::sqrt(pow(latvec.e11, 2) + pow(latvec.e12, 2)
5864
+ pow(latvec.e13, 2));
5965
double celldm3 = std::sqrt(pow(latvec.e31, 2) + pow(latvec.e32, 2)
6066
+ pow(latvec.e33, 2));
@@ -69,19 +75,21 @@ void remake_cell(Lattice& lat)
6975
latvec.e31 = 0.0;
7076
latvec.e32 = 0.0;
7177
latvec.e33 = celldm3;
72-
} else if (latName == "trigonal") // ibrav = 5
73-
{
74-
double celldm1 = std::sqrt(pow(latvec.e11, 2) + pow(latvec.e12, 2)
78+
}
79+
else if (latName == "trigonal") // ibrav = 5
80+
{
81+
double celldm1 = std::sqrt(pow(latvec.e11, 2) + pow(latvec.e12, 2)
7582
+ pow(latvec.e13, 2));
7683
double celldm2 = std::sqrt(pow(latvec.e21, 2) + pow(latvec.e22, 2)
7784
+ pow(latvec.e23, 2));
7885
double celldm12 = (latvec.e11 * latvec.e21 + latvec.e12 * latvec.e22
7986
+ latvec.e13 * latvec.e23);
8087
double cos12 = celldm12 / celldm1 / celldm2;
8188

82-
if (cos12 <= -0.5 || cos12 >= 1.0) {
83-
ModuleBase::WARNING_QUIT("unitcell", "wrong cos12!");
84-
}
89+
if (cos12 <= -0.5 || cos12 >= 1.0)
90+
{
91+
ModuleBase::WARNING_QUIT("unitcell", "wrong cos12!");
92+
}
8593
double t1 = sqrt(1.0 + 2.0 * cos12);
8694
double t2 = sqrt(1.0 - cos12);
8795

@@ -99,8 +107,9 @@ void remake_cell(Lattice& lat)
99107
latvec.e31 = -e11;
100108
latvec.e32 = e12;
101109
latvec.e33 = e13;
102-
} else if (latName == "st") // ibrav = 6
103-
{
110+
}
111+
else if (latName == "st") // ibrav = 6
112+
{
104113
double celldm1 = std::sqrt(pow(latvec.e11, 2) + pow(latvec.e12, 2)
105114
+ pow(latvec.e13, 2));
106115
double celldm3 = std::sqrt(pow(latvec.e31, 2) + pow(latvec.e32, 2)
@@ -114,8 +123,9 @@ void remake_cell(Lattice& lat)
114123
latvec.e31 = 0.0;
115124
latvec.e32 = 0.0;
116125
latvec.e33 = celldm3;
117-
} else if (latName == "bct") // ibrav = 7
118-
{
126+
}
127+
else if (latName == "bct") // ibrav = 7
128+
{
119129
double celldm1 = std::abs(latvec.e11);
120130
double celldm2 = std::abs(latvec.e13);
121131

@@ -128,8 +138,9 @@ void remake_cell(Lattice& lat)
128138
latvec.e31 = -celldm1;
129139
latvec.e32 = -celldm1;
130140
latvec.e33 = celldm2;
131-
} else if (latName == "so") // ibrav = 8
132-
{
141+
}
142+
else if (latName == "so") // ibrav = 8
143+
{
133144
double celldm1 = std::sqrt(pow(latvec.e11, 2) + pow(latvec.e12, 2)
134145
+ pow(latvec.e13, 2));
135146
double celldm2 = std::sqrt(pow(latvec.e21, 2) + pow(latvec.e22, 2)
@@ -146,8 +157,9 @@ void remake_cell(Lattice& lat)
146157
latvec.e31 = 0.0;
147158
latvec.e32 = 0.0;
148159
latvec.e33 = celldm3;
149-
} else if (latName == "baco") // ibrav = 9
150-
{
160+
}
161+
else if (latName == "baco") // ibrav = 9
162+
{
151163
double celldm1 = std::abs(latvec.e11);
152164
double celldm2 = std::abs(latvec.e22);
153165
double celldm3 = std::abs(latvec.e33);
@@ -161,8 +173,9 @@ void remake_cell(Lattice& lat)
161173
latvec.e31 = 0.0;
162174
latvec.e32 = 0.0;
163175
latvec.e33 = celldm3;
164-
} else if (latName == "fco") // ibrav = 10
165-
{
176+
}
177+
else if (latName == "fco") // ibrav = 10
178+
{
166179
double celldm1 = std::abs(latvec.e11);
167180
double celldm2 = std::abs(latvec.e22);
168181
double celldm3 = std::abs(latvec.e33);
@@ -176,8 +189,9 @@ void remake_cell(Lattice& lat)
176189
latvec.e31 = 0.0;
177190
latvec.e32 = celldm2;
178191
latvec.e33 = celldm3;
179-
} else if (latName == "bco") // ibrav = 11
180-
{
192+
}
193+
else if (latName == "bco") // ibrav = 11
194+
{
181195
double celldm1 = std::abs(latvec.e11);
182196
double celldm2 = std::abs(latvec.e12);
183197
double celldm3 = std::abs(latvec.e13);
@@ -191,8 +205,9 @@ void remake_cell(Lattice& lat)
191205
latvec.e31 = -celldm1;
192206
latvec.e32 = -celldm2;
193207
latvec.e33 = celldm3;
194-
} else if (latName == "sm") // ibrav = 12
195-
{
208+
}
209+
else if (latName == "sm") // ibrav = 12
210+
{
196211
double celldm1 = std::sqrt(pow(latvec.e11, 2) + pow(latvec.e12, 2)
197212
+ pow(latvec.e13, 2));
198213
double celldm2 = std::sqrt(pow(latvec.e21, 2) + pow(latvec.e22, 2)
@@ -215,16 +230,18 @@ void remake_cell(Lattice& lat)
215230
latvec.e31 = 0.0;
216231
latvec.e32 = 0.0;
217232
latvec.e33 = celldm3;
218-
} else if (latName == "bacm") // ibrav = 13
219-
{
233+
}
234+
else if (latName == "bacm") // ibrav = 13
235+
{
220236
double celldm1 = std::abs(latvec.e11);
221237
double celldm2 = std::sqrt(pow(latvec.e21, 2) + pow(latvec.e22, 2)
222238
+ pow(latvec.e23, 2));
223239
double celldm3 = std::abs(latvec.e13);
224240

225241
double cos12 = latvec.e21 / celldm2;
226-
if (cos12 >= 1.0) {
227-
ModuleBase::WARNING_QUIT("unitcell", "wrong cos12!");
242+
if (cos12 >= 1.0)
243+
{
244+
ModuleBase::WARNING_QUIT("unitcell", "wrong cos12!");
228245
}
229246

230247
double e21 = celldm2 * cos12;
@@ -239,8 +256,9 @@ void remake_cell(Lattice& lat)
239256
latvec.e31 = celldm1;
240257
latvec.e32 = 0.0;
241258
latvec.e33 = celldm3;
242-
} else if (latName == "triclinic") // ibrav = 14
243-
{
259+
}
260+
else if (latName == "triclinic") // ibrav = 14
261+
{
244262
double celldm1 = std::sqrt(pow(latvec.e11, 2) + pow(latvec.e12, 2)
245263
+ pow(latvec.e13, 2));
246264
double celldm2 = std::sqrt(pow(latvec.e21, 2) + pow(latvec.e22, 2)
@@ -258,8 +276,9 @@ void remake_cell(Lattice& lat)
258276
double cos23 = celldm23 / celldm2 / celldm3;
259277

260278
double sin12 = std::sqrt(1.0 - cos12 * cos12);
261-
if (cos12 >= 1.0) {
262-
ModuleBase::WARNING_QUIT("unitcell", "wrong cos12!");
279+
if (cos12 >= 1.0)
280+
{
281+
ModuleBase::WARNING_QUIT("unitcell", "wrong cos12!");
263282
}
264283

265284
latvec.e11 = celldm1;
@@ -278,15 +297,16 @@ void remake_cell(Lattice& lat)
278297
else
279298
{
280299
std::cout << "latname is : " << latName << std::endl;
281-
ModuleBase::WARNING_QUIT("UnitCell::remake_cell",
282-
"latname not supported!");
300+
ModuleBase::WARNING_QUIT("unitcell::remake_cell",
301+
"latname type not supported!");
283302
}
284303
}
285304

286305
// LiuXh add a new function here,
287306
// 20180515
288-
void setup_cell_after_vc(UnitCell& ucell, std::ofstream& log) {
289-
ModuleBase::TITLE("UnitCell", "setup_cell_after_vc");
307+
void setup_cell_after_vc(UnitCell& ucell, std::ofstream& log)
308+
{
309+
ModuleBase::TITLE("unitcell", "setup_cell_after_vc");
290310
assert(ucell.lat0 > 0.0);
291311
ucell.omega = std::abs(ucell.latvec.Det()) *
292312
pow(ucell.lat0, 3);
@@ -325,9 +345,11 @@ void setup_cell_after_vc(UnitCell& ucell, std::ofstream& log) {
325345
ucell.GGT = ucell.G * ucell.GT;
326346
ucell.invGGT = ucell.GGT.Inverse();
327347

328-
for (int it = 0; it < ucell.ntype; it++) {
329-
Atom* atom = &ucell.atoms[it];
330-
for (int ia = 0; ia < atom->na; ia++) {
348+
for (int it = 0; it < ucell.ntype; it++)
349+
{
350+
Atom* atom = &ucell.atoms[it];
351+
for (int ia = 0; ia < atom->na; ia++)
352+
{
331353
atom->tau[ia] = atom->taud[ia] * ucell.latvec;
332354
}
333355
}
@@ -354,10 +376,13 @@ void update_pos_tau(const Lattice& lat,
354376
Atom* atoms)
355377
{
356378
int iat = 0;
357-
for (int it = 0; it < ntype; it++) {
358-
Atom* atom = &atoms[it];
359-
for (int ia = 0; ia < atom->na; ia++) {
360-
for (int ik = 0; ik < 3; ++ik) {
379+
for (int it = 0; it < ntype; it++)
380+
{
381+
Atom* atom = &atoms[it];
382+
for (int ia = 0; ia < atom->na; ia++)
383+
{
384+
for (int ik = 0; ik < 3; ++ik)
385+
{
361386
if (atom->mbl[ia][ik])
362387
{
363388
atom->dis[ia][ik] = pos[3 * iat + ik] / lat.lat0 - atom->tau[ia][ik];
@@ -454,17 +479,23 @@ void periodic_boundary_adjustment(Atom* atoms,
454479
// first adjust direct coordinates,
455480
// then update them into cartesian coordinates,
456481
//----------------------------------------------
457-
for (int it = 0; it < ntype; it++) {
458-
Atom* atom = &atoms[it];
459-
for (int ia = 0; ia < atom->na; ia++) {
482+
for (int it = 0; it < ntype; it++)
483+
{
484+
Atom* atom = &atoms[it];
485+
for (int ia = 0; ia < atom->na; ia++)
486+
{
460487
// mohan update 2011-03-21
461488
for (int ik = 0; ik < 3; ik++)
462489
{
463-
if (atom->taud[ia][ik] < 0)
490+
// important update: mohan update 2025-07-14, change 'if' to 'while'
491+
// so I guess the following warning will not show up anymore, which
492+
// I am not sure is correct or not. However, I prefer the current
493+
// implementation
494+
while (atom->taud[ia][ik] < 0)
464495
{
465496
atom->taud[ia][ik] += 1.0;
466497
}
467-
if (atom->taud[ia][ik] >= 1.0)
498+
while (atom->taud[ia][ik] >= 1.0)
468499
{
469500
atom->taud[ia][ik] -= 1.0;
470501
}
@@ -476,12 +507,12 @@ void periodic_boundary_adjustment(Atom* atoms,
476507
|| atom->taud[ia].y >= 1.0
477508
|| atom->taud[ia].z >= 1.0)
478509
{
479-
GlobalV::ofs_warning << " it=" << it + 1 << " ia=" << ia + 1 << std::endl;
480-
GlobalV::ofs_warning << "d=" << atom->taud[ia].x << " "
510+
GlobalV::ofs_warning << " atom type=" << it + 1 << " atom index=" << ia + 1 << std::endl;
511+
GlobalV::ofs_warning << " direct coordinate=" << atom->taud[ia].x << " "
481512
<< atom->taud[ia].y << " "
482513
<< atom->taud[ia].z << std::endl;
483-
ModuleBase::WARNING_QUIT("Ions_Move_Basic::move_ions",
484-
"the movement of atom is larger than the length of cell.");
514+
ModuleBase::WARNING_QUIT("unitcell::periodic_boundary_adjustment",
515+
"Movement of atom is larger than the cell length");
485516
}
486517

487518
atom->tau[ia] = atom->taud[ia] * latvec;

‎tests/02_NAO_Gamma/009_NO_GO_CS_CR/INPUT‎

Lines changed: 12 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -1,27 +1,29 @@
11
INPUT_PARAMETERS
22
#Parameters (General)
33
suffix autotest
4-
pseudo_dir ../../PP_ORB
5-
orbital_dir ../../PP_ORB
4+
pseudo_dir ../../PP_ORB
5+
orbital_dir ../../PP_ORB
66

77
nbands 8
88
calculation cell-relax
99
#Parameters (Accuracy)
1010
ecutwfc 20
11-
scf_nmax 20
11+
scf_nmax 20
1212

1313
basis_type lcao
14-
relax_nmax 2
1514

16-
cal_stress 1
17-
stress_thr 1e-6
18-
cal_force 1
19-
force_thr_ev 1.0e-3
15+
cal_stress 1
16+
stress_thr 0.5
17+
cal_force 1
18+
force_thr_ev 0.02
2019

2120
ks_solver scalapack_gvx
2221
mixing_type broyden
2322
mixing_beta 0.7
2423

25-
gamma_only 1
24+
gamma_only 1
2625

27-
relax_new 0
26+
relax_method cg
27+
relax_new 1
28+
relax_nmax 2
29+
relax_scale_force 0.4

‎tests/02_NAO_Gamma/009_NO_GO_CS_CR/STRU‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -13,7 +13,7 @@ LATTICE_VECTORS
1313
0.5 0.5 0.0 #Lattice vector 3
1414

1515
ATOMIC_POSITIONS
16-
Cartesian #Cartesian(Unit is LATTICE_CONSTANT)
16+
Cartesian #Cartesian(Unit is LATTICE_CONSTANT)
1717
Si #Name of element
1818
0.0 #Magnetic for this element.
1919
2 #Number of atoms

‎tests/02_NAO_Gamma/009_NO_GO_CS_CR/result.ref‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,5 @@
1-
etotref -196.6829913753032599
2-
etotperatomref -98.3414956877
1+
etotref -202.4208617343214485
2+
etotperatomref -101.2104308672
33
totalforceref 0.000000
4-
totalstressref 1376.164122
5-
totaltimeref +0.95083
4+
totalstressref 625.031700
5+
totaltimeref 1.81

0 commit comments

Comments
 (0)