@@ -204,7 +204,7 @@ void LCAO_Deepks_Interface<TK, TR>::out_deepks_labels(const double& etot,
204204 }
205205
206206 // Bandgap Part
207- if (PARAM.inp .deepks_bandgap )
207+ if (PARAM.inp .deepks_bandgap == 1 )
208208 {
209209 const int nocc = (PARAM.inp .nelec + 1 ) / 2 ;
210210 ModuleBase::matrix o_tot (nks, 1 );
@@ -273,6 +273,268 @@ void LCAO_Deepks_Interface<TK, TR>::out_deepks_labels(const double& etot,
273273 } // end deepks_out_labels == 1
274274 } // end bandgap label
275275
276+ if (PARAM.inp .deepks_bandgap == 2 )
277+ {
278+ const int nocc = (PARAM.inp .nelec + 1 ) / 2 ;
279+ const int range = PARAM.inp .deepks_band_range [1 ] - PARAM.inp .deepks_band_range [0 ];
280+ ModuleBase::matrix o_tot (nks, range);
281+ for (int iks = 0 ; iks < nks; ++iks)
282+ {
283+ // record band gap for each k point (including spin)
284+ for (int ib = 0 ; ib <= range; ++ib)
285+ {
286+ if (ib + PARAM.inp .deepks_band_range [0 ] < -1 )
287+ {
288+ o_tot (iks, ib) = ekb (iks, nocc + ib + PARAM.inp .deepks_band_range [0 ]) - ekb (iks, nocc - 1 );
289+ }
290+ if (ib + PARAM.inp .deepks_band_range [0 ] > -1 )
291+ {
292+ o_tot (iks, ib - 1 ) = ekb (iks, nocc + ib + PARAM.inp .deepks_band_range [0 ]) - ekb (iks, nocc -1 );
293+ }
294+ }
295+ }
296+
297+ const std::string file_otot
298+ = PARAM.globalv .global_out_dir
299+ + (PARAM.inp .deepks_out_labels == 1 ? " deepks_otot.npy" : " deepks_orbital.npy" );
300+ LCAO_deepks_io::save_matrix2npy (file_otot, o_tot, rank); // Unit: Hartree
301+
302+ if (PARAM.inp .deepks_out_labels == 1 ) // don't need these when deepks_out_labels == 2
303+ {
304+ if (PARAM.inp .deepks_scf )
305+ {
306+ std::vector<ModuleBase::matrix> wg_hl_range (range);
307+ std::vector<std::vector<TH>> dm_bandgap_range (range);
308+
309+ // Calculate O_delta
310+ for (int ir = 0 ; ir <= range; ++ir)
311+ {
312+ if (ir + PARAM.inp .deepks_band_range [0 ] < -1 )
313+ {
314+ wg_hl_range[ir].create (nks, PARAM.inp .nbands );
315+ wg_hl_range[ir].zero_out ();
316+ dm_bandgap_range[ir].resize (nks);
317+ for (int iks = 0 ; iks < nks; ++iks)
318+ {
319+ wg_hl_range[ir](iks, nocc - 1 ) = -1.0 ;
320+ wg_hl_range[ir](iks, nocc + ir + PARAM.inp .deepks_band_range [0 ]) = 1.0 ;
321+ }
322+ elecstate::cal_dm (ParaV, wg_hl_range[ir], psi, dm_bandgap_range[ir]);
323+ }
324+ if (ir + PARAM.inp .deepks_band_range [0 ] > -1 )
325+ {
326+ wg_hl_range[ir - 1 ].create (nks, PARAM.inp .nbands );
327+ wg_hl_range[ir - 1 ].zero_out ();
328+ dm_bandgap_range[ir - 1 ].resize (nks);
329+ for (int iks = 0 ; iks < nks; ++iks)
330+ {
331+ wg_hl_range[ir - 1 ](iks, nocc - 1 ) = -1.0 ;
332+ wg_hl_range[ir - 1 ](iks, nocc + ir + PARAM.inp .deepks_band_range [0 ]) = 1.0 ;
333+ }
334+ elecstate::cal_dm (ParaV, wg_hl_range[ir - 1 ], psi, dm_bandgap_range[ir - 1 ]);
335+ }
336+ }
337+
338+ ModuleBase::matrix o_delta (nks, range);
339+ torch::Tensor orbital_precalc_range;
340+
341+ for (int ir = 0 ; ir < range; ++ir)
342+ {
343+ torch::Tensor orbital_precalc_temp;
344+ ModuleBase::matrix o_delta_temp (nks, 1 );
345+ DeePKS_domain::cal_orbital_precalc<TK, TH>(dm_bandgap_range[ir],
346+ lmaxd,
347+ inlmax,
348+ nat,
349+ nks,
350+ inl2l,
351+ kvec_d,
352+ phialpha,
353+ gevdm,
354+ inl_index,
355+ ucell,
356+ orb,
357+ *ParaV,
358+ GridD,
359+ orbital_precalc_temp);
360+ if (ir == 0 )
361+ {
362+ orbital_precalc_range = orbital_precalc_temp;
363+ }
364+ else
365+ {
366+ orbital_precalc_range = torch::cat ({orbital_precalc_range, orbital_precalc_temp}, 0 );
367+ }
368+
369+ DeePKS_domain::cal_o_delta<TK, TH>(dm_bandgap_range[ir], *h_delta, o_delta_temp, *ParaV, nks, nspin);
370+ for (int iks = 0 ; iks < nks; ++iks)
371+ {
372+ o_delta (iks, ir) = o_delta_temp (iks, 0 );
373+ }
374+ }
375+
376+ // save obase and orbital_precalc
377+ const std::string file_orbpre = PARAM.globalv .global_out_dir + " deepks_orbpre.npy" ;
378+ LCAO_deepks_io::save_tensor2npy<double >(file_orbpre, orbital_precalc_range, rank);
379+
380+ const std::string file_obase = PARAM.globalv .global_out_dir + " deepks_obase.npy" ;
381+ LCAO_deepks_io::save_matrix2npy (file_obase, o_tot - o_delta, rank); // Unit: Hartree
382+ } // end deepks_scf == 1
383+ else // deepks_scf == 0
384+ {
385+ const std::string file_obase = PARAM.globalv .global_out_dir + " deepks_obase.npy" ;
386+ LCAO_deepks_io::save_matrix2npy (file_obase, o_tot, rank); // no scf, o_tot=o_base
387+ } // end deepks_scf == 0
388+ } // end deepks_out_labels == 1
389+ } // end bandgap label
390+
391+ if (PARAM.inp .deepks_bandgap == 3 )
392+ {
393+ const int nocc = (PARAM.inp .nelec + 1 ) / 2 ;
394+ ModuleBase::matrix o_tot (nks, 1 );
395+ for (int iks = 0 ; iks < nks; ++iks)
396+ {
397+ // record band gap for each k point (including spin)
398+ o_tot (iks, 0 ) = ekb (iks, nocc + PARAM.inp .deepks_band_range [1 ]) - ekb (iks, nocc + PARAM.inp .deepks_band_range [0 ]);
399+ }
400+
401+ const std::string file_otot
402+ = PARAM.globalv .global_out_dir
403+ + (PARAM.inp .deepks_out_labels == 1 ? " deepks_otot.npy" : " deepks_orbital.npy" );
404+ LCAO_deepks_io::save_matrix2npy (file_otot, o_tot, rank); // Unit: Hartree
405+
406+ if (PARAM.inp .deepks_out_labels == 1 ) // don't need these when deepks_out_labels == 2
407+ {
408+ if (PARAM.inp .deepks_scf )
409+ {
410+ ModuleBase::matrix wg_hl;
411+ std::vector<TH> dm_bandgap;
412+
413+ // Calculate O_delta
414+ wg_hl.create (nks, PARAM.inp .nbands );
415+ dm_bandgap.resize (nks);
416+ wg_hl.zero_out ();
417+ for (int iks = 0 ; iks < nks; ++iks)
418+ {
419+ wg_hl (iks, nocc + PARAM.inp .deepks_band_range [0 ]) = -1.0 ;
420+ wg_hl (iks, nocc + PARAM.inp .deepks_band_range [1 ]) = 1.0 ;
421+ }
422+ elecstate::cal_dm (ParaV, wg_hl, psi, dm_bandgap);
423+
424+ ModuleBase::matrix o_delta (nks, 1 );
425+
426+ // calculate and save orbital_precalc: [nks,NAt,NDscrpt]
427+ torch::Tensor orbital_precalc;
428+ DeePKS_domain::cal_orbital_precalc<TK, TH>(dm_bandgap,
429+ lmaxd,
430+ inlmax,
431+ nat,
432+ nks,
433+ inl2l,
434+ kvec_d,
435+ phialpha,
436+ gevdm,
437+ inl_index,
438+ ucell,
439+ orb,
440+ *ParaV,
441+ GridD,
442+ orbital_precalc);
443+ DeePKS_domain::cal_o_delta<TK, TH>(dm_bandgap, *h_delta, o_delta, *ParaV, nks, nspin);
444+
445+ // save obase and orbital_precalc
446+ const std::string file_orbpre = PARAM.globalv .global_out_dir + " deepks_orbpre.npy" ;
447+ LCAO_deepks_io::save_tensor2npy<double >(file_orbpre, orbital_precalc, rank);
448+
449+ const std::string file_obase = PARAM.globalv .global_out_dir + " deepks_obase.npy" ;
450+ LCAO_deepks_io::save_matrix2npy (file_obase, o_tot - o_delta, rank); // Unit: Hartree
451+ } // end deepks_scf == 1
452+ else // deepks_scf == 0
453+ {
454+ const std::string file_obase = PARAM.globalv .global_out_dir + " deepks_obase.npy" ;
455+ LCAO_deepks_io::save_matrix2npy (file_obase, o_tot, rank); // no scf, o_tot=o_base
456+ } // end deepks_scf == 0
457+ } // end deepks_out_labels == 1
458+ } // end bandgap label
459+
460+ if (PARAM.inp .deepks_bandgap == 4 )
461+ {
462+ int natom_H = 0 ;
463+ for (int it = 0 ; it < ucell.ntype ; it++)
464+ {
465+ if (ucell.atoms [it].label == " H" )
466+ {
467+ natom_H = ucell.atoms [it].na ;
468+ break ;
469+ }
470+ }
471+ const int nocc = (PARAM.inp .nelec - natom_H) / 2 ;
472+ ModuleBase::matrix o_tot (nks, 1 );
473+ for (int iks = 0 ; iks < nks; ++iks)
474+ {
475+ // record band gap for each k point (including spin)
476+ o_tot (iks, 0 ) = ekb (iks, nocc) - ekb (iks, nocc - 1 );
477+ }
478+
479+ const std::string file_otot
480+ = PARAM.globalv .global_out_dir
481+ + (PARAM.inp .deepks_out_labels == 1 ? " deepks_otot.npy" : " deepks_orbital.npy" );
482+ LCAO_deepks_io::save_matrix2npy (file_otot, o_tot, rank); // Unit: Hartree
483+
484+ if (PARAM.inp .deepks_out_labels == 1 ) // don't need these when deepks_out_labels == 2
485+ {
486+ if (PARAM.inp .deepks_scf )
487+ {
488+ ModuleBase::matrix wg_hl;
489+ std::vector<TH> dm_bandgap;
490+
491+ // Calculate O_delta
492+ wg_hl.create (nks, PARAM.inp .nbands );
493+ dm_bandgap.resize (nks);
494+ wg_hl.zero_out ();
495+ for (int iks = 0 ; iks < nks; ++iks)
496+ {
497+ wg_hl (iks, nocc - 1 ) = -1.0 ;
498+ wg_hl (iks, nocc) = 1.0 ;
499+ }
500+ elecstate::cal_dm (ParaV, wg_hl, psi, dm_bandgap);
501+
502+ ModuleBase::matrix o_delta (nks, 1 );
503+
504+ // calculate and save orbital_precalc: [nks,NAt,NDscrpt]
505+ torch::Tensor orbital_precalc;
506+ DeePKS_domain::cal_orbital_precalc<TK, TH>(dm_bandgap,
507+ lmaxd,
508+ inlmax,
509+ nat,
510+ nks,
511+ inl2l,
512+ kvec_d,
513+ phialpha,
514+ gevdm,
515+ inl_index,
516+ ucell,
517+ orb,
518+ *ParaV,
519+ GridD,
520+ orbital_precalc);
521+ DeePKS_domain::cal_o_delta<TK, TH>(dm_bandgap, *h_delta, o_delta, *ParaV, nks, nspin);
522+
523+ // save obase and orbital_precalc
524+ const std::string file_orbpre = PARAM.globalv .global_out_dir + " deepks_orbpre.npy" ;
525+ LCAO_deepks_io::save_tensor2npy<double >(file_orbpre, orbital_precalc, rank);
526+
527+ const std::string file_obase = PARAM.globalv .global_out_dir + " deepks_obase.npy" ;
528+ LCAO_deepks_io::save_matrix2npy (file_obase, o_tot - o_delta, rank); // Unit: Hartree
529+ } // end deepks_scf == 1
530+ else // deepks_scf == 0
531+ {
532+ const std::string file_obase = PARAM.globalv .global_out_dir + " deepks_obase.npy" ;
533+ LCAO_deepks_io::save_matrix2npy (file_obase, o_tot, rank); // no scf, o_tot=o_base
534+ } // end deepks_scf == 0
535+ } // end deepks_out_labels == 1
536+ } // end bandgap label
537+
276538 // not add deepks_out_labels = 2 for HR yet
277539 // H(R) matrix part, for HR, base will not be calculated since they are HContainer objects
278540 if (PARAM.inp .deepks_v_delta < 0 )
0 commit comments