@@ -278,4 +278,69 @@ void remake_cell(Lattice& lat)
278278 ModuleBase::WARNING_QUIT (" UnitCell::read_atom_species" ,
279279 " latname not supported!" );
280280 }
281- }
281+ }
282+
283+ // LiuXh add a new function here,
284+ // 20180515
285+ void setup_cell_after_vc (UnitCell& ucell, std::ofstream& log) {
286+ ModuleBase::TITLE (" UnitCell" , " setup_cell_after_vc" );
287+ assert (ucell.lat0 > 0.0 );
288+ ucell.omega = std::abs (ucell.latvec .Det ()) *
289+ ucell.lat0 * ucell.lat0 * ucell.lat0 ;
290+ if (ucell.omega <= 0 ) {
291+ ModuleBase::WARNING_QUIT (" setup_cell_after_vc" , " omega <= 0 ." );
292+ } else {
293+ log << std::endl;
294+ ModuleBase::GlobalFunc::OUT (log, " Volume (Bohr^3)" , ucell.omega );
295+ ModuleBase::GlobalFunc::OUT (log, " Volume (A^3)" ,
296+ ucell.omega * pow (ModuleBase::BOHR_TO_A, 3 ));
297+ }
298+
299+ ucell.lat0_angstrom = ucell.lat0 * 0.529177 ;
300+ ucell.tpiba = ModuleBase::TWO_PI / ucell.lat0 ;
301+ ucell.tpiba2 = ucell.tpiba * ucell.tpiba ;
302+
303+ // lattice vectors in another form.
304+ ucell.a1 .x = ucell.latvec .e11 ;
305+ ucell.a1 .y = ucell.latvec .e12 ;
306+ ucell.a1 .z = ucell.latvec .e13 ;
307+
308+ ucell.a2 .x = ucell.latvec .e21 ;
309+ ucell.a2 .y = ucell.latvec .e22 ;
310+ ucell.a2 .z = ucell.latvec .e23 ;
311+
312+ ucell.a3 .x = ucell.latvec .e31 ;
313+ ucell.a3 .y = ucell.latvec .e32 ;
314+ ucell.a3 .z = ucell.latvec .e33 ;
315+
316+ // ==========================================================
317+ // Calculate recip. lattice vectors and dot products
318+ // latvec has the unit of lat0, but G has the unit 2Pi/lat0
319+ // ==========================================================
320+ ucell.GT = ucell.latvec .Inverse ();
321+ ucell.G = ucell.GT .Transpose ();
322+ ucell.GGT = ucell.G * ucell.GT ;
323+ ucell.invGGT = ucell.GGT .Inverse ();
324+
325+ for (int it = 0 ; it < ucell.ntype ; it++) {
326+ Atom* atom = &ucell.atoms [it];
327+ for (int ia = 0 ; ia < atom->na ; ia++) {
328+ atom->tau [ia] = atom->taud [ia] * ucell.latvec ;
329+ }
330+ }
331+
332+ #ifdef __MPI
333+ ucell.bcast_unitcell ();
334+ #endif
335+
336+ log << std::endl;
337+ output::printM3 (log,
338+ " Lattice vectors: (Cartesian coordinate: in unit of a_0)" ,
339+ ucell.latvec );
340+ output::printM3 (
341+ log,
342+ " Reciprocal vectors: (Cartesian coordinate: in unit of 2 pi/a_0)" ,
343+ ucell.G );
344+
345+ return ;
346+ }
0 commit comments