1
1
#include " update_cell.h"
2
2
#include " bcast_cell.h"
3
3
#include " source_base/global_function.h"
4
+
4
5
namespace unitcell
5
6
{
7
+
6
8
void remake_cell (Lattice& lat)
7
9
{
8
10
ModuleBase::TITLE (" Lattice" , " rmake_cell" );
@@ -13,18 +15,20 @@ void remake_cell(Lattice& lat)
13
15
std::string& latName = lat.latName ;
14
16
ModuleBase::Matrix3& latvec = lat.latvec ;
15
17
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" )
20
19
{
20
+ ModuleBase::WARNING_QUIT (" UnitCell" , " to use fixed_ibrav, latname must be provided" );
21
+ }
22
+ else if (latName == " sc" ) // ibrav = 1
23
+ {
21
24
double celldm = std::sqrt (pow (latvec.e11 , 2 ) + pow (latvec.e12 , 2 )
22
25
+ pow (latvec.e13 , 2 ));
23
26
24
27
latvec.Zero ();
25
28
latvec.e11 = latvec.e22 = latvec.e33 = celldm;
26
- } else if (latName == " fcc" ) // ibrav = 2
27
- {
29
+ }
30
+ else if (latName == " fcc" ) // ibrav = 2
31
+ {
28
32
double celldm = std::sqrt (pow (latvec.e11 , 2 ) + pow (latvec.e12 , 2 )
29
33
+ pow (latvec.e13 , 2 )) / std::sqrt (2.0 );
30
34
@@ -37,8 +41,9 @@ void remake_cell(Lattice& lat)
37
41
latvec.e31 = -celldm;
38
42
latvec.e32 = celldm;
39
43
latvec.e33 = 0.0 ;
40
- } else if (latName == " bcc" ) // ibrav = 3
41
- {
44
+ }
45
+ else if (latName == " bcc" ) // ibrav = 3
46
+ {
42
47
double celldm = std::sqrt (pow (latvec.e11 , 2 ) + pow (latvec.e12 , 2 )
43
48
+ pow (latvec.e13 , 2 ))
44
49
/ std::sqrt (3.0 );
@@ -52,9 +57,10 @@ void remake_cell(Lattice& lat)
52
57
latvec.e31 = -celldm;
53
58
latvec.e32 = -celldm;
54
59
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 )
58
64
+ pow (latvec.e13 , 2 ));
59
65
double celldm3 = std::sqrt (pow (latvec.e31 , 2 ) + pow (latvec.e32 , 2 )
60
66
+ pow (latvec.e33 , 2 ));
@@ -69,19 +75,21 @@ void remake_cell(Lattice& lat)
69
75
latvec.e31 = 0.0 ;
70
76
latvec.e32 = 0.0 ;
71
77
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 )
75
82
+ pow (latvec.e13 , 2 ));
76
83
double celldm2 = std::sqrt (pow (latvec.e21 , 2 ) + pow (latvec.e22 , 2 )
77
84
+ pow (latvec.e23 , 2 ));
78
85
double celldm12 = (latvec.e11 * latvec.e21 + latvec.e12 * latvec.e22
79
86
+ latvec.e13 * latvec.e23 );
80
87
double cos12 = celldm12 / celldm1 / celldm2;
81
88
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
+ }
85
93
double t1 = sqrt (1.0 + 2.0 * cos12);
86
94
double t2 = sqrt (1.0 - cos12);
87
95
@@ -99,8 +107,9 @@ void remake_cell(Lattice& lat)
99
107
latvec.e31 = -e11 ;
100
108
latvec.e32 = e12 ;
101
109
latvec.e33 = e13 ;
102
- } else if (latName == " st" ) // ibrav = 6
103
- {
110
+ }
111
+ else if (latName == " st" ) // ibrav = 6
112
+ {
104
113
double celldm1 = std::sqrt (pow (latvec.e11 , 2 ) + pow (latvec.e12 , 2 )
105
114
+ pow (latvec.e13 , 2 ));
106
115
double celldm3 = std::sqrt (pow (latvec.e31 , 2 ) + pow (latvec.e32 , 2 )
@@ -114,8 +123,9 @@ void remake_cell(Lattice& lat)
114
123
latvec.e31 = 0.0 ;
115
124
latvec.e32 = 0.0 ;
116
125
latvec.e33 = celldm3;
117
- } else if (latName == " bct" ) // ibrav = 7
118
- {
126
+ }
127
+ else if (latName == " bct" ) // ibrav = 7
128
+ {
119
129
double celldm1 = std::abs (latvec.e11 );
120
130
double celldm2 = std::abs (latvec.e13 );
121
131
@@ -128,8 +138,9 @@ void remake_cell(Lattice& lat)
128
138
latvec.e31 = -celldm1;
129
139
latvec.e32 = -celldm1;
130
140
latvec.e33 = celldm2;
131
- } else if (latName == " so" ) // ibrav = 8
132
- {
141
+ }
142
+ else if (latName == " so" ) // ibrav = 8
143
+ {
133
144
double celldm1 = std::sqrt (pow (latvec.e11 , 2 ) + pow (latvec.e12 , 2 )
134
145
+ pow (latvec.e13 , 2 ));
135
146
double celldm2 = std::sqrt (pow (latvec.e21 , 2 ) + pow (latvec.e22 , 2 )
@@ -146,8 +157,9 @@ void remake_cell(Lattice& lat)
146
157
latvec.e31 = 0.0 ;
147
158
latvec.e32 = 0.0 ;
148
159
latvec.e33 = celldm3;
149
- } else if (latName == " baco" ) // ibrav = 9
150
- {
160
+ }
161
+ else if (latName == " baco" ) // ibrav = 9
162
+ {
151
163
double celldm1 = std::abs (latvec.e11 );
152
164
double celldm2 = std::abs (latvec.e22 );
153
165
double celldm3 = std::abs (latvec.e33 );
@@ -161,8 +173,9 @@ void remake_cell(Lattice& lat)
161
173
latvec.e31 = 0.0 ;
162
174
latvec.e32 = 0.0 ;
163
175
latvec.e33 = celldm3;
164
- } else if (latName == " fco" ) // ibrav = 10
165
- {
176
+ }
177
+ else if (latName == " fco" ) // ibrav = 10
178
+ {
166
179
double celldm1 = std::abs (latvec.e11 );
167
180
double celldm2 = std::abs (latvec.e22 );
168
181
double celldm3 = std::abs (latvec.e33 );
@@ -176,8 +189,9 @@ void remake_cell(Lattice& lat)
176
189
latvec.e31 = 0.0 ;
177
190
latvec.e32 = celldm2;
178
191
latvec.e33 = celldm3;
179
- } else if (latName == " bco" ) // ibrav = 11
180
- {
192
+ }
193
+ else if (latName == " bco" ) // ibrav = 11
194
+ {
181
195
double celldm1 = std::abs (latvec.e11 );
182
196
double celldm2 = std::abs (latvec.e12 );
183
197
double celldm3 = std::abs (latvec.e13 );
@@ -191,8 +205,9 @@ void remake_cell(Lattice& lat)
191
205
latvec.e31 = -celldm1;
192
206
latvec.e32 = -celldm2;
193
207
latvec.e33 = celldm3;
194
- } else if (latName == " sm" ) // ibrav = 12
195
- {
208
+ }
209
+ else if (latName == " sm" ) // ibrav = 12
210
+ {
196
211
double celldm1 = std::sqrt (pow (latvec.e11 , 2 ) + pow (latvec.e12 , 2 )
197
212
+ pow (latvec.e13 , 2 ));
198
213
double celldm2 = std::sqrt (pow (latvec.e21 , 2 ) + pow (latvec.e22 , 2 )
@@ -215,16 +230,18 @@ void remake_cell(Lattice& lat)
215
230
latvec.e31 = 0.0 ;
216
231
latvec.e32 = 0.0 ;
217
232
latvec.e33 = celldm3;
218
- } else if (latName == " bacm" ) // ibrav = 13
219
- {
233
+ }
234
+ else if (latName == " bacm" ) // ibrav = 13
235
+ {
220
236
double celldm1 = std::abs (latvec.e11 );
221
237
double celldm2 = std::sqrt (pow (latvec.e21 , 2 ) + pow (latvec.e22 , 2 )
222
238
+ pow (latvec.e23 , 2 ));
223
239
double celldm3 = std::abs (latvec.e13 );
224
240
225
241
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!" );
228
245
}
229
246
230
247
double e21 = celldm2 * cos12;
@@ -239,8 +256,9 @@ void remake_cell(Lattice& lat)
239
256
latvec.e31 = celldm1;
240
257
latvec.e32 = 0.0 ;
241
258
latvec.e33 = celldm3;
242
- } else if (latName == " triclinic" ) // ibrav = 14
243
- {
259
+ }
260
+ else if (latName == " triclinic" ) // ibrav = 14
261
+ {
244
262
double celldm1 = std::sqrt (pow (latvec.e11 , 2 ) + pow (latvec.e12 , 2 )
245
263
+ pow (latvec.e13 , 2 ));
246
264
double celldm2 = std::sqrt (pow (latvec.e21 , 2 ) + pow (latvec.e22 , 2 )
@@ -258,8 +276,9 @@ void remake_cell(Lattice& lat)
258
276
double cos23 = celldm23 / celldm2 / celldm3;
259
277
260
278
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!" );
263
282
}
264
283
265
284
latvec.e11 = celldm1;
@@ -278,15 +297,16 @@ void remake_cell(Lattice& lat)
278
297
else
279
298
{
280
299
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!" );
283
302
}
284
303
}
285
304
286
305
// LiuXh add a new function here,
287
306
// 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" );
290
310
assert (ucell.lat0 > 0.0 );
291
311
ucell.omega = std::abs (ucell.latvec .Det ()) *
292
312
pow (ucell.lat0 , 3 );
@@ -325,9 +345,11 @@ void setup_cell_after_vc(UnitCell& ucell, std::ofstream& log) {
325
345
ucell.GGT = ucell.G * ucell.GT ;
326
346
ucell.invGGT = ucell.GGT .Inverse ();
327
347
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
+ {
331
353
atom->tau [ia] = atom->taud [ia] * ucell.latvec ;
332
354
}
333
355
}
@@ -354,10 +376,13 @@ void update_pos_tau(const Lattice& lat,
354
376
Atom* atoms)
355
377
{
356
378
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
+ {
361
386
if (atom->mbl [ia][ik])
362
387
{
363
388
atom->dis [ia][ik] = pos[3 * iat + ik] / lat.lat0 - atom->tau [ia][ik];
@@ -454,9 +479,11 @@ void periodic_boundary_adjustment(Atom* atoms,
454
479
// first adjust direct coordinates,
455
480
// then update them into cartesian coordinates,
456
481
// ----------------------------------------------
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
+ {
460
487
// mohan update 2011-03-21
461
488
for (int ik = 0 ; ik < 3 ; ik++)
462
489
{
@@ -476,12 +503,12 @@ void periodic_boundary_adjustment(Atom* atoms,
476
503
|| atom->taud [ia].y >= 1.0
477
504
|| atom->taud [ia].z >= 1.0 )
478
505
{
479
- GlobalV::ofs_warning << " it =" << it + 1 << " ia =" << ia + 1 << std::endl;
480
- GlobalV::ofs_warning << " d =" << atom->taud [ia].x << " "
506
+ GlobalV::ofs_warning << " atom type =" << it + 1 << " atom index =" << ia + 1 << std::endl;
507
+ GlobalV::ofs_warning << " direct coordinate =" << atom->taud [ia].x << " "
481
508
<< atom->taud [ia].y << " "
482
509
<< 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. " );
510
+ ModuleBase::WARNING_QUIT (" unitcell::periodic_boundary_adjustment " ,
511
+ " Movement of atom is larger than the cell length " );
485
512
}
486
513
487
514
atom->tau [ia] = atom->taud [ia] * latvec;
0 commit comments