diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index 4b9fadc8e3..f392740f83 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -2027,10 +2027,13 @@ Warning: this function is not robust enough for the current version. Please try ### deepks_out_labels -- **Type**: Boolean +- **Type**: Integer - **Availability**: numerical atomic orbital basis -- **Description**: Print labels and descriptors for DeePKS training in OUT.${suffix}. The names of these files start with "deepks". -- **Note**: In `LCAO` calculation, the path of a numerical descriptor (an `orb` file) is needed to be specified under the `NUMERICAL_DESCRIPTOR` tag in the `STRU` file. For example: +- **Description**: Print labels and descriptors for DeePKS in OUT.${suffix}. The names of these files start with "deepks". + - 0 : No output. + - 1 : Output intermediate files needed during DeePKS training. + - 2 : Output target labels for label preperation. The label files are named as `deepks_.npy`, where the units and formats are the same as label files `.npy` required for training, except that the first dimension (`nframes`) is excluded. System structrue files are also given in `deepks_atom.npy` and `deepks_box.npy` in the unit of *Bohr*, which means `lattice_constant` should be set to 1 when training. +- **Note**: When `deepks_out_labels` equals **1**, the path of a numerical descriptor (an `orb` file) is needed to be specified under the `NUMERICAL_DESCRIPTOR` tag in the `STRU` file. For example: ```text NUMERICAL_ORBITAL @@ -2040,8 +2043,8 @@ Warning: this function is not robust enough for the current version. Please try NUMERICAL_DESCRIPTOR jle.orb ``` - -- **Default**: False + This is not needed when `deepks_out_labels` equals 2. +- **Default**: 0 ### deepks_scf diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp index 683420fff3..75e0feea1e 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp @@ -504,17 +504,23 @@ void Force_Stress_LCAO::getForceStress(UnitCell& ucell, // DeePKS force if (PARAM.inp.deepks_out_labels) // not parallelized yet { - const std::string file_ftot = PARAM.globalv.global_out_dir + "deepks_ftot.npy"; - LCAO_deepks_io::save_matrix2npy(file_ftot, fcs, GlobalV::MY_RANK); // Ry/Bohr, F_tot + const std::string file_ftot = PARAM.globalv.global_out_dir + + (PARAM.inp.deepks_out_labels == 1 ? "deepks_ftot.npy" : "deepks_force.npy"); + LCAO_deepks_io::save_matrix2npy(file_ftot, fcs, GlobalV::MY_RANK); // Hartree/Bohr, F_tot - const std::string file_fbase = PARAM.globalv.global_out_dir + "deepks_fbase.npy"; - if (PARAM.inp.deepks_scf) + if (PARAM.inp.deepks_out_labels == 1) { - LCAO_deepks_io::save_matrix2npy(file_fbase, fcs - fvnl_dalpha, GlobalV::MY_RANK); // Ry/Bohr, F_base - } - else - { - LCAO_deepks_io::save_matrix2npy(file_fbase, fcs, GlobalV::MY_RANK); // no scf, F_base=F_tot + const std::string file_fbase = PARAM.globalv.global_out_dir + "deepks_fbase.npy"; + if (PARAM.inp.deepks_scf) + { + LCAO_deepks_io::save_matrix2npy(file_fbase, + fcs - fvnl_dalpha, + GlobalV::MY_RANK); // Hartree/Bohr, F_base + } + else + { + LCAO_deepks_io::save_matrix2npy(file_fbase, fcs, GlobalV::MY_RANK); // no scf, F_base=F_tot + } } } #endif @@ -686,25 +692,38 @@ void Force_Stress_LCAO::getForceStress(UnitCell& ucell, #ifdef __DEEPKS if (PARAM.inp.deepks_out_labels) // not parallelized yet { - const std::string file_stot = PARAM.globalv.global_out_dir + "deepks_stot.npy"; - LCAO_deepks_io::save_matrix2npy(file_stot, - scs, - GlobalV::MY_RANK, - ucell.omega, - 'U'); // change to energy unit Ry when printing, S_tot; - - const std::string file_sbase = PARAM.globalv.global_out_dir + "deepks_sbase.npy"; - if (PARAM.inp.deepks_scf) + if (PARAM.inp.deepks_out_labels == 1) { - LCAO_deepks_io::save_matrix2npy(file_sbase, - scs - svnl_dalpha, + const std::string file_stot = PARAM.globalv.global_out_dir + "deepks_stot.npy"; + LCAO_deepks_io::save_matrix2npy(file_stot, + scs, GlobalV::MY_RANK, ucell.omega, - 'U'); // change to energy unit Ry when printing, S_base; + 'U'); // change to energy unit Ry when printing, S_tot; + + const std::string file_sbase = PARAM.globalv.global_out_dir + "deepks_sbase.npy"; + if (PARAM.inp.deepks_scf) + { + LCAO_deepks_io::save_matrix2npy(file_sbase, + scs - svnl_dalpha, + GlobalV::MY_RANK, + ucell.omega, + 'U'); // change to energy unit Ry when printing, S_base; + } + else + { + LCAO_deepks_io::save_matrix2npy(file_sbase, + scs, + GlobalV::MY_RANK, + ucell.omega, + 'U'); // sbase = stot + } } - else + else if (PARAM.inp.deepks_out_labels == 2) { - LCAO_deepks_io::save_matrix2npy(file_sbase, scs, GlobalV::MY_RANK, ucell.omega, 'U'); // sbase = stot + const std::string file_stot = PARAM.globalv.global_out_dir + "deepks_stress.npy"; + LCAO_deepks_io::save_matrix2npy(file_stot, scs, GlobalV::MY_RANK, ucell.omega, + 'F'); // flat mode } } #endif diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_allocate.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_allocate.cpp index ed25d09bcf..81c2b9a77b 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_allocate.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_allocate.cpp @@ -18,7 +18,7 @@ void DeePKS_init(const UnitCell& ucell, { ModuleBase::TITLE("LCAO_domain", "DeePKS_init"); // preparation for DeePKS - if (PARAM.inp.deepks_out_labels || PARAM.inp.deepks_scf) + if (PARAM.inp.deepks_out_labels == 1 || PARAM.inp.deepks_scf) { // allocate relevant data structures for calculating descriptors std::vector na; diff --git a/source/module_hamilt_lcao/module_deepks/LCAO_deepks_interface.cpp b/source/module_hamilt_lcao/module_deepks/LCAO_deepks_interface.cpp index 8cc1d03418..c103cf6ecc 100644 --- a/source/module_hamilt_lcao/module_deepks/LCAO_deepks_interface.cpp +++ b/source/module_hamilt_lcao/module_deepks/LCAO_deepks_interface.cpp @@ -60,7 +60,7 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, // Note : update PDM and all other quantities with the current dm // DeePKS PDM and descriptor - if (PARAM.inp.deepks_out_labels || PARAM.inp.deepks_scf) + if (PARAM.inp.deepks_out_labels == 1 || PARAM.inp.deepks_scf) { // this part is for integrated test of deepks // so it is printed no matter even if deepks_out_labels is not used @@ -80,17 +80,14 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, descriptor, rank); - if (PARAM.inp.deepks_out_labels) - { - LCAO_deepks_io::save_npy_d(nat, - des_per_atom, - inlmax, - inl2l, - PARAM.inp.deepks_equiv, - descriptor, - PARAM.globalv.global_out_dir, - rank); // libnpy needed - } + LCAO_deepks_io::save_npy_d(nat, + des_per_atom, + inlmax, + inl2l, + PARAM.inp.deepks_equiv, + descriptor, + PARAM.globalv.global_out_dir, + rank); // libnpy needed if (PARAM.inp.deepks_scf) { @@ -138,26 +135,29 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, } // Energy Part - const std::string file_etot = PARAM.globalv.global_out_dir + "deepks_etot.npy"; - const std::string file_ebase = PARAM.globalv.global_out_dir + "deepks_ebase.npy"; - + const std::string file_etot = PARAM.globalv.global_out_dir + + (PARAM.inp.deepks_out_labels == 1 ? "deepks_etot.npy" : "deepks_energy.npy"); LCAO_deepks_io::save_npy_e(etot, file_etot, rank); - if (PARAM.inp.deepks_scf) - { - /// ebase :no deepks E_delta including - LCAO_deepks_io::save_npy_e(etot - E_delta, file_ebase, rank); - } - else // deepks_scf = 0; base calculation + if (PARAM.inp.deepks_out_labels == 1) { - /// no scf, e_tot=e_base - LCAO_deepks_io::save_npy_e(etot, file_ebase, rank); + const std::string file_ebase = PARAM.globalv.global_out_dir + "deepks_ebase.npy"; + if (PARAM.inp.deepks_scf) + { + /// ebase :no deepks E_delta including + LCAO_deepks_io::save_npy_e(etot - E_delta, file_ebase, rank); + } + else // deepks_scf = 0; base calculation + { + /// no scf, e_tot=e_base + LCAO_deepks_io::save_npy_e(etot, file_ebase, rank); + } } // Force Part if (PARAM.inp.cal_force) { - if (PARAM.inp.deepks_scf + if (PARAM.inp.deepks_scf && PARAM.inp.deepks_out_labels == 1 // don't need these when deepks_out_labels == 2 && !PARAM.inp.deepks_equiv) // training with force label not supported by equivariant version now { torch::Tensor gdmx; @@ -180,7 +180,7 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, // Stress Part if (PARAM.inp.cal_stress) { - if (PARAM.inp.deepks_scf + if (PARAM.inp.deepks_scf && PARAM.inp.deepks_out_labels == 1 // don't need these when deepks_out_labels == 2 && !PARAM.inp.deepks_equiv) // training with stress label not supported by equivariant version now { torch::Tensor gdmepsl; @@ -211,60 +211,66 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, o_tot(iks, 0) = ekb(iks, nocc) - ekb(iks, nocc - 1); } - const std::string file_otot = PARAM.globalv.global_out_dir + "deepks_otot.npy"; - LCAO_deepks_io::save_matrix2npy(file_otot, o_tot, rank); // Unit: Ry + const std::string file_otot + = PARAM.globalv.global_out_dir + + (PARAM.inp.deepks_out_labels == 1 ? "deepks_otot.npy" : "deepks_orbital.npy"); + LCAO_deepks_io::save_matrix2npy(file_otot, o_tot, rank); // Unit: Hartree - if (PARAM.inp.deepks_scf) + if (PARAM.inp.deepks_out_labels == 1) // don't need these when deepks_out_labels == 2 { - ModuleBase::matrix wg_hl; - std::vector dm_bandgap; - - // Calculate O_delta - wg_hl.create(nks, PARAM.inp.nbands); - dm_bandgap.resize(nks); - wg_hl.zero_out(); - for (int iks = 0; iks < nks; ++iks) + if (PARAM.inp.deepks_scf) { - wg_hl(iks, nocc - 1) = -1.0; - wg_hl(iks, nocc) = 1.0; - } - elecstate::cal_dm(ParaV, wg_hl, psi, dm_bandgap); - - ModuleBase::matrix o_delta(nks, 1); - - // calculate and save orbital_precalc: [nks,NAt,NDscrpt] - torch::Tensor orbital_precalc; - DeePKS_domain::cal_orbital_precalc(dm_bandgap, - lmaxd, - inlmax, - nat, - nks, - inl2l, - kvec_d, - phialpha, - gevdm, - inl_index, - ucell, - orb, - *ParaV, - GridD, - orbital_precalc); - DeePKS_domain::cal_o_delta(dm_bandgap, *h_delta, o_delta, *ParaV, nks, nspin); - - // save obase and orbital_precalc - const std::string file_orbpre = PARAM.globalv.global_out_dir + "deepks_orbpre.npy"; - LCAO_deepks_io::save_tensor2npy(file_orbpre, orbital_precalc, rank); - - const std::string file_obase = PARAM.globalv.global_out_dir + "deepks_obase.npy"; - LCAO_deepks_io::save_matrix2npy(file_obase, o_tot - o_delta, rank); // Unit: Ry - } // end deepks_scf == 1 - else // deepks_scf == 0 - { - const std::string file_obase = PARAM.globalv.global_out_dir + "deepks_obase.npy"; - LCAO_deepks_io::save_matrix2npy(file_obase, o_tot, rank); // no scf, o_tot=o_base - } // end deepks_scf == 0 - } // end bandgap label + ModuleBase::matrix wg_hl; + std::vector dm_bandgap; + + // Calculate O_delta + wg_hl.create(nks, PARAM.inp.nbands); + dm_bandgap.resize(nks); + wg_hl.zero_out(); + for (int iks = 0; iks < nks; ++iks) + { + wg_hl(iks, nocc - 1) = -1.0; + wg_hl(iks, nocc) = 1.0; + } + elecstate::cal_dm(ParaV, wg_hl, psi, dm_bandgap); + + ModuleBase::matrix o_delta(nks, 1); + + // calculate and save orbital_precalc: [nks,NAt,NDscrpt] + torch::Tensor orbital_precalc; + DeePKS_domain::cal_orbital_precalc(dm_bandgap, + lmaxd, + inlmax, + nat, + nks, + inl2l, + kvec_d, + phialpha, + gevdm, + inl_index, + ucell, + orb, + *ParaV, + GridD, + orbital_precalc); + DeePKS_domain::cal_o_delta(dm_bandgap, *h_delta, o_delta, *ParaV, nks, nspin); + + // save obase and orbital_precalc + const std::string file_orbpre = PARAM.globalv.global_out_dir + "deepks_orbpre.npy"; + LCAO_deepks_io::save_tensor2npy(file_orbpre, orbital_precalc, rank); + + const std::string file_obase = PARAM.globalv.global_out_dir + "deepks_obase.npy"; + LCAO_deepks_io::save_matrix2npy(file_obase, o_tot - o_delta, rank); // Unit: Hartree + } // end deepks_scf == 1 + else // deepks_scf == 0 + { + const std::string file_obase = PARAM.globalv.global_out_dir + "deepks_obase.npy"; + LCAO_deepks_io::save_matrix2npy(file_obase, o_tot, rank); // no scf, o_tot=o_base + } // end deepks_scf == 0 + } // end deepks_out_labels == 1 + } // end bandgap label + // not add deepks_out_labels = 2 for HR yet // H(R) matrix part, for HR, base will not be calculated since they are HContainer objects if (PARAM.inp.deepks_v_delta < 0) { @@ -328,92 +334,120 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, if (PARAM.inp.deepks_v_delta > 0) { std::vector h_tot(nks); - std::vector> h_mat(nks, std::vector(ParaV->nloc)); - for (int ik = 0; ik < nks; ik++) - { - h_tot[ik].create(nlocal, nlocal); - - p_ham->updateHk(ik); - - const TK* hk_ptr = p_ham->getHk(); - - for (int i = 0; i < ParaV->nloc; i++) - { - h_mat[ik][i] = hk_ptr[i]; - } - } - - DeePKS_domain::collect_h_mat(*ParaV, h_mat, h_tot, nlocal, nks); + DeePKS_domain::get_h_tot(*ParaV, p_ham, h_tot, nlocal, nks, 'H'); - const std::string file_htot = PARAM.globalv.global_out_dir + "deepks_htot.npy"; + const std::string file_htot + = PARAM.globalv.global_out_dir + + (PARAM.inp.deepks_out_labels == 1 ? "deepks_htot.npy" : "deepks_hamiltonian.npy"); LCAO_deepks_io::save_npy_h(h_tot, file_htot, nlocal, nks, rank); - if (PARAM.inp.deepks_scf) + if (PARAM.inp.deepks_out_labels == 1) // don't need these when deepks_out_labels == 2 { - std::vector v_delta(nks); - std::vector h_base(nks); - for (int ik = 0; ik < nks; ik++) + if (PARAM.inp.deepks_scf) { - v_delta[ik].create(nlocal, nlocal); - h_base[ik].create(nlocal, nlocal); + std::vector v_delta(nks); + std::vector h_base(nks); + for (int ik = 0; ik < nks; ik++) + { + v_delta[ik].create(nlocal, nlocal); + h_base[ik].create(nlocal, nlocal); + } + DeePKS_domain::collect_h_mat(*ParaV, *h_delta, v_delta, nlocal, nks); + + // save v_delta and h_base + const std::string file_hbase = PARAM.globalv.global_out_dir + "deepks_hbase.npy"; + for (int ik = 0; ik < nks; ik++) + { + h_base[ik] = h_tot[ik] - v_delta[ik]; + } + LCAO_deepks_io::save_npy_h(h_base, file_hbase, nlocal, nks, rank); + + const std::string file_vdelta = PARAM.globalv.global_out_dir + "deepks_vdelta.npy"; + LCAO_deepks_io::save_npy_h(v_delta, file_vdelta, nlocal, nks, rank); + + if (PARAM.inp.deepks_v_delta == 1) // v_delta_precalc storage method 1 + { + torch::Tensor v_delta_precalc; + DeePKS_domain::cal_v_delta_precalc(nlocal, + lmaxd, + inlmax, + nat, + nks, + inl2l, + kvec_d, + phialpha, + gevdm, + inl_index, + ucell, + orb, + *ParaV, + GridD, + v_delta_precalc); + + const std::string file_vdpre = PARAM.globalv.global_out_dir + "deepks_vdpre.npy"; + LCAO_deepks_io::save_tensor2npy(file_vdpre, v_delta_precalc, rank); + } + else if (PARAM.inp.deepks_v_delta == 2) // v_delta_precalc storage method 2 + { + torch::Tensor phialpha_out; + DeePKS_domain::prepare_phialpha(nlocal, + lmaxd, + inlmax, + nat, + nks, + kvec_d, + phialpha, + ucell, + orb, + *ParaV, + GridD, + phialpha_out); + const std::string file_phialpha = PARAM.globalv.global_out_dir + "deepks_phialpha.npy"; + LCAO_deepks_io::save_tensor2npy(file_phialpha, phialpha_out, rank); + + torch::Tensor gevdm_out; + DeePKS_domain::prepare_gevdm(nat, lmaxd, inlmax, orb, gevdm, gevdm_out); + const std::string file_gevdm = PARAM.globalv.global_out_dir + "deepks_gevdm.npy"; + LCAO_deepks_io::save_tensor2npy(file_gevdm, gevdm_out, rank); + } } - DeePKS_domain::collect_h_mat(*ParaV, *h_delta, v_delta, nlocal, nks); - - // save v_delta and h_base - const std::string file_hbase = PARAM.globalv.global_out_dir + "deepks_hbase.npy"; - for (int ik = 0; ik < nks; ik++) + else // deepks_scf == 0 { - h_base[ik] = h_tot[ik] - v_delta[ik]; + const std::string file_hbase = PARAM.globalv.global_out_dir + "deepks_hbase.npy"; + LCAO_deepks_io::save_npy_h(h_tot, file_hbase, nlocal, nks, rank); } - LCAO_deepks_io::save_npy_h(h_base, file_hbase, nlocal, nks, rank); + } // end deepks_out_labels == 1 + } // end v_delta label - const std::string file_vdelta = PARAM.globalv.global_out_dir + "deepks_vdelta.npy"; - LCAO_deepks_io::save_npy_h(v_delta, file_vdelta, nlocal, nks, rank); + } // end deepks_out_labels - if (PARAM.inp.deepks_v_delta == 1) // v_delta_precalc storage method 1 - { - torch::Tensor v_delta_precalc; - DeePKS_domain::cal_v_delta_precalc(nlocal, - lmaxd, - inlmax, - nat, - nks, - inl2l, - kvec_d, - phialpha, - gevdm, - inl_index, - ucell, - orb, - *ParaV, - GridD, - v_delta_precalc); - - const std::string file_vdpre = PARAM.globalv.global_out_dir + "deepks_vdpre.npy"; - LCAO_deepks_io::save_tensor2npy(file_vdpre, v_delta_precalc, rank); - } - else if (PARAM.inp.deepks_v_delta == 2) // v_delta_precalc storage method 2 - { - torch::Tensor phialpha_out; - DeePKS_domain::prepare_phialpha< - TK>(nlocal, lmaxd, inlmax, nat, nks, kvec_d, phialpha, ucell, orb, *ParaV, GridD, phialpha_out); - const std::string file_phialpha = PARAM.globalv.global_out_dir + "deepks_phialpha.npy"; - LCAO_deepks_io::save_tensor2npy(file_phialpha, phialpha_out, rank); - - torch::Tensor gevdm_out; - DeePKS_domain::prepare_gevdm(nat, lmaxd, inlmax, orb, gevdm, gevdm_out); - const std::string file_gevdm = PARAM.globalv.global_out_dir + "deepks_gevdm.npy"; - LCAO_deepks_io::save_tensor2npy(file_gevdm, gevdm_out, rank); - } - } - else // deepks_scf == 0 - { - const std::string file_hbase = PARAM.globalv.global_out_dir + "deepks_hbase.npy"; - LCAO_deepks_io::save_npy_h(h_tot, file_hbase, nlocal, nks, rank); - } - } // end v_delta label + if (PARAM.inp.deepks_out_labels == 2) + { + // output atom.npy and box.npy + torch::Tensor atom_out; + DeePKS_domain::prepare_atom(ucell, atom_out); + const std::string file_atom = PARAM.globalv.global_out_dir + "deepks_atom.npy"; + LCAO_deepks_io::save_tensor2npy(file_atom, atom_out, rank); - } // end deepks_out_labels + torch::Tensor box_out; + DeePKS_domain::prepare_box(ucell, box_out); + const std::string file_box = PARAM.globalv.global_out_dir + "deepks_box.npy"; + LCAO_deepks_io::save_tensor2npy(file_box, box_out, rank); + + if (PARAM.inp.deepks_v_delta > 0) + { + // prepare for overlap.npy, very much like h_tot except for p_ham->getSk() + std::vector s_tot(nks); + DeePKS_domain::get_h_tot(*ParaV, p_ham, s_tot, nlocal, nks, 'S'); + const std::string file_stot = PARAM.globalv.global_out_dir + "deepks_overlap.npy"; + LCAO_deepks_io::save_npy_h(s_tot, + file_stot, + nlocal, + nks, + rank, + 1.0); // don't need unit_scale for overlap + } + } /// print out deepks information to the screen if (PARAM.inp.deepks_scf) diff --git a/source/module_hamilt_lcao/module_deepks/LCAO_deepks_io.cpp b/source/module_hamilt_lcao/module_deepks/LCAO_deepks_io.cpp index 2620b9a0d6..07a116b9d5 100644 --- a/source/module_hamilt_lcao/module_deepks/LCAO_deepks_io.cpp +++ b/source/module_hamilt_lcao/module_deepks/LCAO_deepks_io.cpp @@ -136,7 +136,7 @@ void LCAO_deepks_io::save_npy_d(const int nat, } // saves energy in numpy format -void LCAO_deepks_io::save_npy_e(const double& e, const std::string& e_file, const int rank) +void LCAO_deepks_io::save_npy_e(const double& e, const std::string& e_file, const int rank, const double unit_scale) { ModuleBase::TITLE("LCAO_deepks_io", "save_npy_e"); if (rank != 0) @@ -146,7 +146,8 @@ void LCAO_deepks_io::save_npy_e(const double& e, const std::string& e_file, cons // save energy in .npy format const long unsigned eshape[] = {1}; - npy::SaveArrayAsNumpy(e_file, false, 1, eshape, &e); + double e_hartree = e * unit_scale; + npy::SaveArrayAsNumpy(e_file, false, 1, eshape, &e_hartree); return; } @@ -155,7 +156,8 @@ void LCAO_deepks_io::save_npy_h(const std::vector& hamilt, const std::string& h_file, const int nlocal, const int nks, - const int rank) + const int rank, + const double unit_scale) { ModuleBase::TITLE("LCAO_deepks_io", "save_npy_h"); if (rank != 0) @@ -173,7 +175,7 @@ void LCAO_deepks_io::save_npy_h(const std::vector& hamilt, { for (int j = 0; j < nlocal; j++) { - npy_h.push_back(hamilt[k](i, j)); + npy_h.push_back(hamilt[k](i, j) * unit_scale); } } } @@ -186,7 +188,8 @@ void LCAO_deepks_io::save_matrix2npy(const std::string& file_name, const ModuleBase::matrix& matrix, const int rank, const double& scale, - const char mode) + const char mode, + const double unit_scale) { ModuleBase::TITLE("LCAO_deepks_io", "save_matrix2npy"); @@ -213,6 +216,12 @@ void LCAO_deepks_io::save_matrix2npy(const std::string& file_name, shape[0] = nr; shape[1] = nc; } + else if (mode == 'F') // flat + { + size = nr * nc; + shape.resize(1); + shape[0] = size; + } else { ModuleBase::WARNING_QUIT("save_matrix2npy", "Invalid mode! Support only 'U', 'L', 'N'."); @@ -226,7 +235,7 @@ void LCAO_deepks_io::save_matrix2npy(const std::string& file_name, { for (int j = i; j < nc; ++j) { - scaled_data[index] = matrix(i, j) * scale; + scaled_data[index] = (matrix(i, j) * scale) * unit_scale; index++; } } @@ -238,18 +247,18 @@ void LCAO_deepks_io::save_matrix2npy(const std::string& file_name, { for (int j = 0; j <= i; ++j) { - scaled_data[index] = matrix(i, j) * scale; + scaled_data[index] = (matrix(i, j) * scale) * unit_scale; index++; } } } - else // normal + else // normal or flat { for (int i = 0; i < nr; ++i) { for (int j = 0; j < nc; ++j) { - scaled_data[i * nc + j] = matrix(i, j) * scale; + scaled_data[i * nc + j] = (matrix(i, j) * scale) * unit_scale; } } } @@ -266,7 +275,8 @@ void LCAO_deepks_io::save_tensor2npy(const std::string& file_name, const torch:: { return; } - using T_tensor = typename std::conditional>::value, c10::complex, T>::type; + using T_tensor = + typename std::conditional>::value, c10::complex, T>::type; const int dim = tensor.dim(); std::vector shape(dim); for (int i = 0; i < dim; i++) @@ -300,13 +310,15 @@ template void LCAO_deepks_io::save_npy_h(const std::vector>(const std::vector& hamilt, const std::string& h_file, const int nlocal, const int nks, - const int rank); + const int rank, + const double unit_scale); template void LCAO_deepks_io::save_tensor2npy(const std::string& file_name, const torch::Tensor& tensor, diff --git a/source/module_hamilt_lcao/module_deepks/LCAO_deepks_io.h b/source/module_hamilt_lcao/module_deepks/LCAO_deepks_io.h index 4c536a4999..5c33816961 100644 --- a/source/module_hamilt_lcao/module_deepks/LCAO_deepks_io.h +++ b/source/module_hamilt_lcao/module_deepks/LCAO_deepks_io.h @@ -32,6 +32,10 @@ namespace LCAO_deepks_io /// 4. save_matrix2npy : ModuleBase::matrix -> .npy, for force, stress and orbital /// 5. save_tensor2npy : torch::Tensor -> .npy, for precalculation variables +/// Ry2Hartree : convert Ry to Hartree, for energy, force, stress, orbital and Hamiltonian, which is consistent with +/// deepks-kit. Used in save_npy_e, save_npy_h and save_matrix2npy by default +constexpr double Ry2Hartree = 0.5; + /// print density matrices template void print_dm(const int nks, const int nlocal, const int nrow, const std::vector>& dm); @@ -51,7 +55,8 @@ void save_npy_d(const int nat, // save energy void save_npy_e(const double& e, /**<[in] \f$E_{base}\f$ or \f$E_{tot}\f$, in Ry*/ const std::string& e_file, - const int rank); + const int rank, + const double unit_scale = Ry2Hartree); // save Hamiltonian template @@ -59,13 +64,15 @@ void save_npy_h(const std::vector& hamilt, const std::string& h_file, const int nlocal, const int nks, - const int rank); + const int rank, + const double unit_scale = Ry2Hartree); void save_matrix2npy(const std::string& file_name, const ModuleBase::matrix& matrix, const int rank, const double& scale = 1.0, - const char mode = 'N'); + const char mode = 'N', + const double unit_scale = Ry2Hartree); template void save_tensor2npy(const std::string& file_name, const torch::Tensor& tensor, const int rank); diff --git a/source/module_hamilt_lcao/module_deepks/deepks_basic.cpp b/source/module_hamilt_lcao/module_deepks/deepks_basic.cpp index 7be7f58fc2..e5a4547f35 100644 --- a/source/module_hamilt_lcao/module_deepks/deepks_basic.cpp +++ b/source/module_hamilt_lcao/module_deepks/deepks_basic.cpp @@ -5,6 +5,7 @@ #ifdef __DEEPKS #include "deepks_basic.h" +#include "module_base/atom_in.h" #include "module_base/timer.h" #include "module_parameter/parameter.h" @@ -236,4 +237,42 @@ void DeePKS_domain::check_gedm(const int inlmax, const std::vector& inl2l, } } +void DeePKS_domain::prepare_atom(const UnitCell& ucell, torch::Tensor& atom_out) +{ + int nat = ucell.nat; + atom_out = torch::zeros({nat, 4}, torch::TensorOptions().dtype(torch::kFloat64)); + + // get atom information + atom_in AtomInfo; + int index = 0; + for (int it = 0; it < ucell.ntype; ++it) + { + for (int ia = 0; ia < ucell.atoms[it].na; ++ia) + { + atom_out[index][0] = AtomInfo.atom_Z[ucell.atom_label[it]]; + + // use bohr as unit + atom_out[index][1] = ucell.atoms[it].tau[ia].x * ucell.lat0; + atom_out[index][2] = ucell.atoms[it].tau[ia].y * ucell.lat0; + atom_out[index][3] = ucell.atoms[it].tau[ia].z * ucell.lat0; + index++; + } + } +} +void DeePKS_domain::prepare_box(const UnitCell& ucell, torch::Tensor& box_out) +{ + box_out = torch::zeros({9}, torch::TensorOptions().dtype(torch::kFloat64)); + + // use bohr as unit + box_out[0] = ucell.latvec.e11 * ucell.lat0; + box_out[1] = ucell.latvec.e12 * ucell.lat0; + box_out[2] = ucell.latvec.e13 * ucell.lat0; + box_out[3] = ucell.latvec.e21 * ucell.lat0; + box_out[4] = ucell.latvec.e22 * ucell.lat0; + box_out[5] = ucell.latvec.e23 * ucell.lat0; + box_out[6] = ucell.latvec.e31 * ucell.lat0; + box_out[7] = ucell.latvec.e32 * ucell.lat0; + box_out[8] = ucell.latvec.e33 * ucell.lat0; +} + #endif diff --git a/source/module_hamilt_lcao/module_deepks/deepks_basic.h b/source/module_hamilt_lcao/module_deepks/deepks_basic.h index bd8ff9ff08..db3f88ef4b 100644 --- a/source/module_hamilt_lcao/module_deepks/deepks_basic.h +++ b/source/module_hamilt_lcao/module_deepks/deepks_basic.h @@ -5,6 +5,7 @@ #include "LCAO_deepks_io.h" #include "module_base/parallel_reduce.h" #include "module_base/tool_title.h" +#include "module_cell/unitcell.h" #include #include @@ -23,6 +24,8 @@ namespace DeePKS_domain // caculated using torch::autograd::grad // 4. check_gedm : prints gedm for checking // 5. cal_edelta_gedm_equiv : calculates E_delta and d(E_delta)/d(pdm) for equivariant version +// 6. prepare_atom : prepares atom tensor for output as deepks_out_labels = 2 +// 7. prepare_box : prepares box tensor for output as deepks_out_labels = 2 // load the trained neural network models void load_model(const std::string& model_file, torch::jit::script::Module& model); @@ -56,6 +59,8 @@ void cal_edelta_gedm_equiv(const int nat, double& E_delta, const int rank); +void prepare_atom(const UnitCell& ucell, torch::Tensor& atom_out); +void prepare_box(const UnitCell& ucell, torch::Tensor& box_out); } // namespace DeePKS_domain #endif #endif \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_deepks/deepks_vdelta.cpp b/source/module_hamilt_lcao/module_deepks/deepks_vdelta.cpp index d01a755d6b..fee7d684b0 100644 --- a/source/module_hamilt_lcao/module_deepks/deepks_vdelta.cpp +++ b/source/module_hamilt_lcao/module_deepks/deepks_vdelta.cpp @@ -15,6 +15,7 @@ #include "module_base/parallel_reduce.h" #include "module_base/timer.h" #include "module_base/tool_title.h" +#include "module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h" // calculating sum of correction band energies template @@ -128,6 +129,42 @@ void DeePKS_domain::collect_h_mat(const Parallel_Orbitals& pv, } } +template +void DeePKS_domain::get_h_tot(const Parallel_Orbitals& pv, + hamilt::HamiltLCAO* p_ham, + std::vector& h_tot, + const int nlocal, + const int nks, + const char matrix_type) // 'H' for H(k), 'S' for S(k) +{ + ModuleBase::TITLE("DeePKS_domain", "get_h_tot"); + TK* (hamilt::HamiltLCAO::*getMatrixFunc)() const = nullptr; + if (matrix_type == 'H') + { + getMatrixFunc = &hamilt::HamiltLCAO::getHk; + } + else if (matrix_type == 'S') + { + getMatrixFunc = &hamilt::HamiltLCAO::getSk; + } + else + { + throw std::invalid_argument("Invalid matrix_type. Use 'H' for H(k) or 'S' for S(k)."); + } + std::vector> h_mat(nks, std::vector(pv.nloc)); + for (int ik = 0; ik < nks; ik++) + { + h_tot[ik].create(nlocal, nlocal); + p_ham->updateHk(ik); + const TK* hk_ptr = (p_ham->*getMatrixFunc)(); + for (int i = 0; i < pv.nloc; i++) + { + h_mat[ik][i] = hk_ptr[i]; + } + } + DeePKS_domain::collect_h_mat(pv, h_mat, h_tot, nlocal, nks); +} + template void DeePKS_domain::cal_e_delta_band(const std::vector>& dm, const std::vector>& V_delta, const int nks, @@ -155,4 +192,26 @@ template void DeePKS_domain::collect_h_mat, ModuleBase::Com const int nlocal, const int nks); +template void DeePKS_domain::get_h_tot(const Parallel_Orbitals& pv, + hamilt::HamiltLCAO* p_ham, + std::vector& h_tot, + const int nlocal, + const int nks, + const char matrix_type); + +template void DeePKS_domain::get_h_tot, ModuleBase::ComplexMatrix, double>( + const Parallel_Orbitals& pv, + hamilt::HamiltLCAO, double>* p_ham, + std::vector& h_tot, + const int nlocal, + const int nks, + const char matrix_type); + +template void DeePKS_domain::get_h_tot, ModuleBase::ComplexMatrix, std::complex>( + const Parallel_Orbitals& pv, + hamilt::HamiltLCAO, std::complex>* p_ham, + std::vector& h_tot, + const int nlocal, + const int nks, + const char matrix_type); #endif diff --git a/source/module_hamilt_lcao/module_deepks/deepks_vdelta.h b/source/module_hamilt_lcao/module_deepks/deepks_vdelta.h index 650119d0c3..bef38af7a7 100644 --- a/source/module_hamilt_lcao/module_deepks/deepks_vdelta.h +++ b/source/module_hamilt_lcao/module_deepks/deepks_vdelta.h @@ -6,6 +6,12 @@ #include "module_base/matrix.h" #include "module_basis/module_ao/parallel_orbitals.h" +// break the circular dependency of HamiltLCAO +namespace hamilt +{ +template +class HamiltLCAO; +} namespace DeePKS_domain { //------------------------ @@ -28,10 +34,19 @@ void cal_e_delta_band(const std::vector>& dm, // Collect data in h_in to matrix h_out. Note that left lower trianger in h_out is filled template void collect_h_mat(const Parallel_Orbitals& pv, - const std::vector>& h_in, - std::vector& h_out, - const int nlocal, - const int nks); + const std::vector>& h_in, + std::vector& h_out, + const int nlocal, + const int nks); + +// Get H(k) or S(k) matrix from p_hamilt and store it in h_tot +template +void get_h_tot(const Parallel_Orbitals& pv, + hamilt::HamiltLCAO* p_ham, + std::vector& h_tot, + const int nlocal, + const int nks, + const char matrix_type); // 'H' for H(k), 'S' for S(k) } // namespace DeePKS_domain #endif #endif \ No newline at end of file diff --git a/source/module_io/read_input_item_deepks.cpp b/source/module_io/read_input_item_deepks.cpp index 6d4268a502..82652ca1c5 100644 --- a/source/module_io/read_input_item_deepks.cpp +++ b/source/module_io/read_input_item_deepks.cpp @@ -1,6 +1,6 @@ #include "module_base/constants.h" -#include "module_parameter/parameter.h" #include "module_base/tool_quit.h" +#include "module_parameter/parameter.h" #include "read_input.h" #include "read_input_tool.h" namespace ModuleIO @@ -9,8 +9,8 @@ void ReadInput::item_deepks() { { Input_Item item("deepks_out_labels"); - item.annotation = ">0 compute descriptor for deepks"; - read_sync_bool(input.deepks_out_labels); + item.annotation = ">0 compute descriptor for deepks. 1 used during training, 2 used for label production"; + read_sync_int(input.deepks_out_labels); this->add_item(item); } { @@ -18,13 +18,13 @@ void ReadInput::item_deepks() item.annotation = ">0 add V_delta to Hamiltonian"; read_sync_bool(input.deepks_scf); item.check_value = [](const Input_Item& item, const Parameter& para) { - #ifndef __DEEPKS - if (PARAM.inp.deepks_scf || PARAM.inp.deepks_out_labels || - PARAM.inp.deepks_bandgap || PARAM.inp.deepks_v_delta) +#ifndef __DEEPKS + if (PARAM.inp.deepks_scf || PARAM.inp.deepks_out_labels || PARAM.inp.deepks_bandgap + || PARAM.inp.deepks_v_delta) { ModuleBase::WARNING_QUIT("Input_conv", "please compile with DeePKS"); } - #endif +#endif }; this->add_item(item); } @@ -48,7 +48,8 @@ void ReadInput::item_deepks() } { Input_Item item("deepks_v_delta"); - item.annotation = ">0 for v_delta label. when output, 1 for v_delta_precalc, 2 for phialpha and grad_evdm ( can save memory )"; + item.annotation = ">0 for v_delta label. when output, 1 for v_delta_precalc, 2 for phialpha and grad_evdm ( " + "can save memory )"; read_sync_int(input.deepks_v_delta); this->add_item(item); } @@ -60,16 +61,17 @@ void ReadInput::item_deepks() item.reset_value = [](const Input_Item& item, Parameter& para) { if (para.input.deepks_out_unittest) { - para.input.deepks_out_labels = true; + para.input.deepks_out_labels = 1; para.input.deepks_scf = true; - } }; item.check_value = [](const Input_Item& item, const Parameter& para) { - if (para.input.deepks_out_unittest){ + if (para.input.deepks_out_unittest) + { if (para.input.cal_force != 1 || para.input.cal_stress != 1) { - ModuleBase::WARNING_QUIT("ReadInput", "force and stress are required in generating deepks unittest"); + ModuleBase::WARNING_QUIT("ReadInput", + "force and stress are required in generating deepks unittest"); } } }; diff --git a/source/module_io/read_set_globalv.cpp b/source/module_io/read_set_globalv.cpp index ca840f65e5..83bde5c62d 100644 --- a/source/module_io/read_set_globalv.cpp +++ b/source/module_io/read_set_globalv.cpp @@ -23,7 +23,7 @@ void ReadInput::set_globalv(const Input_para& inp, System_para& sys) } } /// set deepks_setorb - if (inp.deepks_scf || inp.deepks_out_labels) + if (inp.deepks_scf || inp.deepks_out_labels == 1) { sys.deepks_setorb = true; } @@ -75,8 +75,8 @@ void ReadInput::set_globalv(const Input_para& inp, System_para& sys) sys.has_float_data = (inp.precision == "float") || (inp.precision == "mixing") || float_cond; } -/// @note Here para.inp has not been synchronized of all ranks. -/// Only para.inp in rank 0 is right. +/// @note Here para.inp has not been synchronized of all ranks. +/// Only para.inp in rank 0 is right. /// So we need to broadcast the results to all ranks. void ReadInput::set_global_dir(const Input_para& inp, System_para& sys) { diff --git a/source/module_parameter/input_parameter.h b/source/module_parameter/input_parameter.h index b3b08150c7..1a82b79d4d 100644 --- a/source/module_parameter/input_parameter.h +++ b/source/module_parameter/input_parameter.h @@ -33,7 +33,7 @@ struct Input_para int kpar = 1; ///< ecch pool is for one k point int bndpar = 1; ///< parallel for stochastic/deterministic bands std::string latname = "none"; ///< lattice name - double ecutwfc = 0; ///< energy cutoff for wavefunctions + double ecutwfc = 0; ///< energy cutoff for wavefunctions double ecutrho = 0; ///< energy cutoff for charge/potential int nx = 0, ny = 0, nz = 0; ///< three dimension of FFT wavefunc @@ -50,7 +50,7 @@ struct Input_para bool dm_to_rho = false; ///< read density matrix from npz format and calculate charge density std::string chg_extrap = "default"; ///< xiaohui modify 2015-02-01 bool init_vel = false; ///< read velocity from STRU or not liuyu 2021-07-14 - + std::string input_file = "INPUT"; ///< input file name std::string stru_file = "STRU"; ///< file contains atomic positions -- ///< xiaohui modify 2015-02-01 @@ -84,11 +84,11 @@ struct Input_para ///< use mesh, which is used in QE. int nspin = 1; ///< LDA ; LSDA ; non-linear spin int pw_diag_nmax = 50; - double pw_diag_thr = 0.01; ///< used in cg method + double pw_diag_thr = 0.01; ///< used in cg method bool diago_smooth_ethr = false; ///< smooth ethr for iter methods - int pw_diag_ndim = 4; ///< dimension of workspace for Davidson diagonalization - int diago_cg_prec = 1; ///< mohan add 2012-03-31 - int diag_subspace = 0; // 0: Lapack, 1: elpa, 2: scalapack + int pw_diag_ndim = 4; ///< dimension of workspace for Davidson diagonalization + int diago_cg_prec = 1; ///< mohan add 2012-03-31 + int diag_subspace = 0; // 0: Lapack, 1: elpa, 2: scalapack std::string smearing_method = "gauss"; ///< "gauss", ///< "mp","methfessel-paxton" @@ -117,9 +117,9 @@ struct Input_para int scf_thr_type = -1; ///< type of the criterion of scf_thr, 1: reci drho, 2: real drho bool final_scf = false; ///< whether to do final scf bool scf_os_stop = false; ///< whether to stop scf when oscillation is detected - double scf_os_thr = -0.01; ///< drho threshold for oscillation + double scf_os_thr = -0.01; ///< drho threshold for oscillation int scf_os_ndim = 0; ///< number of old iterations used for oscillation detection - int sc_os_ndim = 5; ///< number of old iterations used for oscillation detection in Spin-Constrained DFT + int sc_os_ndim = 5; ///< number of old iterations used for oscillation detection in Spin-Constrained DFT bool lspinorb = false; ///< consider the spin-orbit interaction bool noncolin = false; ///< using non-collinear-spin @@ -209,40 +209,41 @@ struct Input_para std::string of_kernel_file = "WTkernel.txt"; ///< The name of WT kernel file. // ML KEDF, sunliang added on 2022-11-07 - bool of_ml_gene_data = false; ///< Generate training data or not + bool of_ml_gene_data = false; ///< Generate training data or not // device - std::string of_ml_device = "cpu"; ///< Run NN on GPU or CPU - int of_ml_feg = 0; ///< The Free Electron Gas limit: 0: no, 3: yes + std::string of_ml_device = "cpu"; ///< Run NN on GPU or CPU + int of_ml_feg = 0; ///< The Free Electron Gas limit: 0: no, 3: yes // kernel - int of_ml_nkernel = 1; ///< Number of kernels - std::vector of_ml_kernel = {1}; ///< Type of kernel, 1 for wt, 2 for yukawa, and 3 for TKK - std::vector of_ml_kernel_scaling = {1.0}; ///< Scaling parameter of kernel, w(r-r') = lambda^3 * w(lambda (r-r')), lambda = 1/scaling + int of_ml_nkernel = 1; ///< Number of kernels + std::vector of_ml_kernel = {1}; ///< Type of kernel, 1 for wt, 2 for yukawa, and 3 for TKK + std::vector of_ml_kernel_scaling + = {1.0}; ///< Scaling parameter of kernel, w(r-r') = lambda^3 * w(lambda (r-r')), lambda = 1/scaling std::vector of_ml_yukawa_alpha = {1.0}; ///< Parameter alpha of yukawa kernel std::vector of_ml_kernel_file = {"none"}; ///< The file of TKK // semi-local descriptors - bool of_ml_gamma = false; ///< Descriptor: gamma = (rho / rho0)^(1/3) - bool of_ml_p = false; ///< Descriptor: p = |nabla rho|^2 / [2 (3 pi^2)^(1/3) rho^(4/3)]^2 - bool of_ml_q = false; ///< Descriptor: q = nabla^2 rho / [4 (3 pi^2)^(2/3) rho^(5/3)] - bool of_ml_tanhp = false; ///< Descriptor: tanhp = tanh(chi_p * p) - bool of_ml_tanhq = false; ///< Descriptor: tanhq = tanh(chi_q * q) - double of_ml_chi_p = 1.0; ///< Hyperparameter: tanhp = tanh(chi_p * p) - double of_ml_chi_q = 1.0; ///< Hyperparameter: tanhq = tanh(chi_q * q) + bool of_ml_gamma = false; ///< Descriptor: gamma = (rho / rho0)^(1/3) + bool of_ml_p = false; ///< Descriptor: p = |nabla rho|^2 / [2 (3 pi^2)^(1/3) rho^(4/3)]^2 + bool of_ml_q = false; ///< Descriptor: q = nabla^2 rho / [4 (3 pi^2)^(2/3) rho^(5/3)] + bool of_ml_tanhp = false; ///< Descriptor: tanhp = tanh(chi_p * p) + bool of_ml_tanhq = false; ///< Descriptor: tanhq = tanh(chi_q * q) + double of_ml_chi_p = 1.0; ///< Hyperparameter: tanhp = tanh(chi_p * p) + double of_ml_chi_q = 1.0; ///< Hyperparameter: tanhq = tanh(chi_q * q) // non-local descriptors // of_ml_gammanl should be a vector of bool, but here we use a vector of int for convinience - std::vector of_ml_gammanl = {0}; ///< Descriptor: gammanl = int{gamma(r') * w(r-r') dr'} - std::vector of_ml_pnl = {0}; ///< Descriptor: pnl = int{p(r') * w(r-r') dr'} - std::vector of_ml_qnl = {0}; ///< Descriptor: qnl = int{q(r') * w(r-r') dr'} - std::vector of_ml_xi = {0}; ///< Descriptor: xi = int{rho(r')^(1/3) * w(r-r') dr'} / rho^(1/3) - std::vector of_ml_tanhxi = {0}; ///< Descriptor: tanhxi = tanh(chi_xi * xi) - std::vector of_ml_tanhxi_nl = {0}; ///< Descriptor: tanhxi_nl = int{tanhxi(r') * w(r-r') dr'} - std::vector of_ml_tanh_pnl = {0}; ///< Descriptor: tanh_pnl = tanh(chi_pnl * pnl) - std::vector of_ml_tanh_qnl = {0}; ///< Descriptor: tanh_qnl = tanh(chi_qnl * qnl) - std::vector of_ml_tanhp_nl = {0}; ///< Descriptor: tanhp_nl = int{tanhp(r') * w(r-r') dr'} - std::vector of_ml_tanhq_nl = {0}; ///< Descriptor: tanhq_nl = int{tanhq(r') * w(r-r') dr'} - std::vector of_ml_chi_xi = {1.0}; ///< Hyperparameter: tanhpxi = tanh(chi_xi * xi) - std::vector of_ml_chi_pnl = {1.0}; ///< Hyperparameter: tanh_pnl = tanh(chi_pnl * pnl) - std::vector of_ml_chi_qnl = {1.0}; ///< Hyperparameter: tanh_qnl = tanh(chi_qnl * qnl) - bool of_ml_local_test = false; ///< Test: read in the density, and output the F and Pauli potential + std::vector of_ml_gammanl = {0}; ///< Descriptor: gammanl = int{gamma(r') * w(r-r') dr'} + std::vector of_ml_pnl = {0}; ///< Descriptor: pnl = int{p(r') * w(r-r') dr'} + std::vector of_ml_qnl = {0}; ///< Descriptor: qnl = int{q(r') * w(r-r') dr'} + std::vector of_ml_xi = {0}; ///< Descriptor: xi = int{rho(r')^(1/3) * w(r-r') dr'} / rho^(1/3) + std::vector of_ml_tanhxi = {0}; ///< Descriptor: tanhxi = tanh(chi_xi * xi) + std::vector of_ml_tanhxi_nl = {0}; ///< Descriptor: tanhxi_nl = int{tanhxi(r') * w(r-r') dr'} + std::vector of_ml_tanh_pnl = {0}; ///< Descriptor: tanh_pnl = tanh(chi_pnl * pnl) + std::vector of_ml_tanh_qnl = {0}; ///< Descriptor: tanh_qnl = tanh(chi_qnl * qnl) + std::vector of_ml_tanhp_nl = {0}; ///< Descriptor: tanhp_nl = int{tanhp(r') * w(r-r') dr'} + std::vector of_ml_tanhq_nl = {0}; ///< Descriptor: tanhq_nl = int{tanhq(r') * w(r-r') dr'} + std::vector of_ml_chi_xi = {1.0}; ///< Hyperparameter: tanhpxi = tanh(chi_xi * xi) + std::vector of_ml_chi_pnl = {1.0}; ///< Hyperparameter: tanh_pnl = tanh(chi_pnl * pnl) + std::vector of_ml_chi_qnl = {1.0}; ///< Hyperparameter: tanh_qnl = tanh(chi_qnl * qnl) + bool of_ml_local_test = false; ///< Test: read in the density, and output the F and Pauli potential // ============== #Parameters (7.stochastic DFT) =========================== int method_sto = 2; ///< different methods for sdft, 1: slow, less memory 2: @@ -261,17 +262,17 @@ struct Input_para //========================================================== // DeepKS -- added by caoyu and mohan //========================================================== - bool deepks_out_labels = false; ///< (need libnpy) prints energy and force labels and - ///< descriptors for training, wenfei 2022-1-12 - bool deepks_scf = false; ///< (need libnpy and libtorch) if set to true, a trained model - ///< would be needed to calculate V_delta and F_delta - bool deepks_bandgap = false; ///< for bandgap label. QO added 2021-12-15 - int deepks_v_delta = 0; ///< for v_delta label. xuan added - bool deepks_equiv = false; ///< whether to use equivariant version of DeePKS - bool deepks_out_unittest = false; ///< if set to true, prints intermediate quantities that shall - ///< be used for making unit test - std::string deepks_model = "None"; ///< needed when deepks_scf=1 - + int deepks_out_labels = 0; ///< (need libnpy) prints energy and force labels and + ///< descriptors for training, wenfei 2022-1-12 + bool deepks_scf = false; ///< (need libnpy and libtorch) if set to true, a trained model + ///< would be needed to calculate V_delta and F_delta + bool deepks_bandgap = false; ///< for bandgap label. QO added 2021-12-15 + int deepks_v_delta = 0; ///< for v_delta label. xuan added + bool deepks_equiv = false; ///< whether to use equivariant version of DeePKS + bool deepks_out_unittest = false; ///< if set to true, prints intermediate quantities that shall + ///< be used for making unit test + std::string deepks_model = "None"; ///< needed when deepks_scf=1 + int bessel_descriptor_lmax = 2; ///< lmax used in descriptor std::string bessel_descriptor_ecut = "default"; ///< energy cutoff for spherical bessel functions(Ry) double bessel_descriptor_tolerence = 1e-12; ///< tolerance for spherical bessel root @@ -280,8 +281,8 @@ struct Input_para double bessel_descriptor_sigma = 0.1; ///< spherical bessel smearing_sigma // ============== #Parameters (9.rt-tddft) =========================== - double td_force_dt = 0.02; ///<"fs" - bool td_vext = false; ///< add extern potential or not + double td_force_dt = 0.02; ///<"fs" + bool td_vext = false; ///< add extern potential or not // std::string td_vext_dire = "1"; ///< vext direction std::vector td_vext_dire = {1}; ///< vector of vext direction @@ -337,20 +338,23 @@ struct Input_para std::vector ocp_kb = {}; ///< OCP kb values // ============== #Parameters (10.lr-tddft) =========================== - int lr_nstates = 1; ///< the number of 2-particle states to be solved - std::vector lr_init_xc_kernel = {}; ///< The method to initalize the xc kernel - int nocc = -1; ///< the number of occupied orbitals to form the 2-particle basis - int nvirt = 1; ///< the number of virtual orbitals to form the 2-particle basis (nocc + nvirt <= nbands) + int lr_nstates = 1; ///< the number of 2-particle states to be solved + std::vector lr_init_xc_kernel = {}; ///< The method to initalize the xc kernel + int nocc = -1; ///< the number of occupied orbitals to form the 2-particle basis + int nvirt = 1; ///< the number of virtual orbitals to form the 2-particle basis (nocc + nvirt <= nbands) std::string xc_kernel = "LDA"; ///< exchange correlation (XC) kernel for LR-TDDFT std::string lr_solver = "dav"; ///< the eigensolver for LR-TDDFT double lr_thr = 1e-2; ///< convergence threshold of the LR-TDDFT eigensolver bool out_wfc_lr = false; ///< whether to output the eigenvectors (excitation amplitudes) in the particle-hole basis - bool lr_unrestricted = false; ///< whether to use the unrestricted construction for LR-TDDFT + bool lr_unrestricted = false; ///< whether to use the unrestricted construction for LR-TDDFT std::vector abs_wavelen_range = {}; ///< the range of wavelength(nm) to output the absorption spectrum - double abs_broadening = 0.01; ///< the broadening (eta) for LR-TDDFT absorption spectrum - std::string abs_gauge = "length"; ///< whether to use length or velocity gauge to calculate the absorption spectrum in LR-TDDFT - std::string ri_hartree_benchmark = "none"; ///< whether to use the RI approximation for the Hartree potential in LR-TDDFT for benchmark (with FHI-aims/ABACUS read-in style) - std::vector aims_nbasis = {}; ///< the number of basis functions for each atom type used in FHI-aims (for benchmark) + double abs_broadening = 0.01; ///< the broadening (eta) for LR-TDDFT absorption spectrum + std::string abs_gauge + = "length"; ///< whether to use length or velocity gauge to calculate the absorption spectrum in LR-TDDFT + std::string ri_hartree_benchmark = "none"; ///< whether to use the RI approximation for the Hartree potential in + ///< LR-TDDFT for benchmark (with FHI-aims/ABACUS read-in style) + std::vector aims_nbasis + = {}; ///< the number of basis functions for each atom type used in FHI-aims (for benchmark) // ============== #Parameters (11.Output) =========================== bool out_stru = false; ///< outut stru file each ion step int out_freq_elec = 0; ///< the frequency of electronic iter to output charge and wavefunction @@ -372,13 +376,13 @@ struct Input_para bool out_bandgap = false; ///< QO added for bandgap printing std::vector out_mat_hs = {0, 8}; ///< output H matrix and S matrix in local basis. std::vector out_mat_tk = {0, 8}; ///< output T(k) matrix in local basis. - std::vector out_mat_l = {0, 8}; ///< output L matrix in local basis. + std::vector out_mat_l = {0, 8}; ///< output L matrix in local basis. bool out_mat_hs2 = false; ///< LiuXh add 2019-07-16, output H(R) matrix and ///< S(R) matrix in local basis. bool out_mat_dh = false; bool out_mat_xc = false; ///< output exchange-correlation matrix in ///< KS-orbital representation. - bool out_mat_xc2 = false; ///< output exchange-correlation matrix Vxc(R) in NAO representation. + bool out_mat_xc2 = false; ///< output exchange-correlation matrix Vxc(R) in NAO representation. bool out_eband_terms = false; ///< output the band energy terms separately bool out_hr_npz = false; ///< output exchange-correlation matrix in ///< KS-orbital representation. @@ -405,7 +409,7 @@ struct Input_para std::vector out_wfc_norm = {}; ///< specify the bands to be calculated for norm of wfc std::vector out_wfc_re_im = {}; ///< specify the bands to be calculated for real and imaginary parts of wfc bool if_separate_k = false; ///< whether to write partial charge for all k-points to individual files or merge them - std::vector out_elf = {0, 3}; ///< output the electron localization function (ELF). 0: no; 1: yes + std::vector out_elf = {0, 3}; ///< output the electron localization function (ELF). 0: no; 1: yes // ============== #Parameters (12.Postprocess) =========================== double dos_emin_ev = -15.0; @@ -508,42 +512,43 @@ struct Input_para // exx // Peize Lin add 2018-06-20 // ========================================================== - std::string exx_hybrid_alpha = "default"; ///< fraction of Fock exchange in hybrid functionals - double exx_hse_omega = 0.11; ///< range-separation parameter in HSE functional - bool exx_separate_loop = true; ///< if 1, a two-step method is employed, else it will start - ///< with a GGA-Loop, and then Hybrid-Loop - int exx_hybrid_step = 100; ///< the maximal electronic iteration number in - ///< the evaluation of Fock exchange - double exx_mixing_beta = 1.0; ///< mixing_beta for outer-loop when exx_separate_loop=1 - double exx_lambda = 0.3; ///< used to compensate for divergence points at G=0 in the - ///< evaluation of Fock exchange using lcao_in_pw method - std::string exx_real_number = "default"; ///< exx calculated in real or complex - double exx_pca_threshold = 0.0001; ///< threshold to screen on-site ABFs in exx - double exx_c_threshold = 0.0001; ///< threshold to screen C matrix in exx - double exx_v_threshold = 0.1; ///< threshold to screen C matrix in exx - double exx_dm_threshold = 0.0001; ///< threshold to screen density matrix in exx - double exx_schwarz_threshold = 0; ///< threshold to screen exx using Cauchy-Schwartz inequality - double exx_cauchy_threshold = 1e-07; ///< threshold to screen exx using Cauchy-Schwartz inequality - double exx_c_grad_threshold = 0.0001; ///< threshold to screen nabla C matrix in exx - double exx_v_grad_threshold = 0.1; ///< threshold to screen nabla V matrix in exx - double exx_c_grad_r_threshold = 0.0001; ///< threshold to screen nabla C matrix in exx - double exx_v_grad_r_threshold = 0.1; ///< threshold to screen nabla V matrix in exx - double exx_cauchy_force_threshold = 1e-07; ///< threshold to screen exx force using Cauchy-Schwartz - ///< inequality - double exx_cauchy_stress_threshold = 1e-07; ///< threshold to screen exx stress using Cauchy-Schwartz - ///< inequality - std::string exx_ccp_rmesh_times = "default"; ///< how many times larger the radial mesh required for - ///< calculating Columb potential is to that of atomic orbitals - std::string exx_distribute_type = "htime"; ///< distribute type (assuming default as no specific value - ///< provided) - int exx_opt_orb_lmax = 0; ///< the maximum l of the spherical Bessel functions for opt ABFs - double exx_opt_orb_ecut = 0.0; ///< the cut-off of plane wave expansion for opt ABFs - double exx_opt_orb_tolerence = 0.0; ///< the threshold when solving for the zeros of spherical Bessel - ///< functions for opt ABFs - bool exx_symmetry_realspace = true; ///< whether to reduce the real-space sector in when using symmetry=1 in EXX calculation - double rpa_ccp_rmesh_times = 10.0; ///< how many times larger the radial mesh required for - ///< calculating Columb potential is to that of atomic orbitals - bool out_ri_cv = false; /// des_tmp.txt + total_des=`sum_file des_tmp.txt 5` + rm des_tmp.txt + echo "totaldes $total_des" >>$1 + if ! test -z "$deepks_scf" && [ $deepks_scf == 1 ]; then + deepks_e_dm=`python3 get_dm_eig.py` + echo "deepks_e_dm $deepks_e_dm" >>$1 + fi + if ! test -z "$has_force" && [ $has_force == 1 ]; then + deepks_f_label=`python3 get_grad_vx.py` + echo "deepks_f_label $deepks_f_label" >>$1 + fi + if ! test -z "$has_stress" && [ $has_stress == 1 ]; then + deepks_s_label=`python3 get_grad_vepsl.py` + echo "deepks_s_label $deepks_s_label" >>$1 + fi +fi +if ! test -z "$deepks_out_labels" && [ $deepks_out_labels == 2 ]; then + deepks_atom=`python3 get_sum_numpy.py OUT.autotest/deepks_atom.npy ` + echo "deepks_atom $deepks_atom" >> $1 + deepks_box=`python3 get_sum_numpy.py OUT.autotest/deepks_box.npy ` + echo "deepks_box $deepks_box" >> $1 + deepks_energy=`python3 get_sum_numpy.py OUT.autotest/deepks_energy.npy ` + echo "deepks_energy $deepks_energy" >> $1 + if ! test -z "$has_force" && [ $has_force == 1 ]; then + deepks_force=`python3 get_sum_numpy.py OUT.autotest/deepks_force.npy ` + echo "deepks_force $deepks_force" >> $1 + fi + if ! test -z "$has_stress" && [ $has_stress == 1 ]; then + deepks_stress=`python3 get_sum_numpy.py OUT.autotest/deepks_stress.npy ` + echo "deepks_stress $deepks_stress" >> $1 + fi + if ! test -z "$deepks_bandgap" && [ $deepks_bandgap == 1 ]; then + deepks_orbital=`python3 get_sum_numpy.py OUT.autotest/deepks_orbital.npy ` + echo "deepks_orbital $deepks_orbital" >> $1 + fi + if ! test -z "$deepks_v_delta" && [[ $deepks_v_delta -gt 0 ]]; then + deepks_hamiltonian=`python3 get_sum_numpy.py OUT.autotest/deepks_hamiltonian.npy ` + echo "deepks_hamiltonian $deepks_hamiltonian" >> $1 + deepks_overlap=`python3 get_sum_numpy.py OUT.autotest/deepks_overlap.npy ` + echo "deepks_overlap $deepks_overlap" >> $1 + fi +fi + +#-------------------------------------------- +# band gap information +#-------------------------------------------- +if ! test -z "$deepks_bandgap" && [ $deepks_bandgap == 1 ] && [ $deepks_out_labels == 1 ]; then + odelta=`python3 get_odelta.py` + echo "odelta $odelta" >>$1 + oprec=`python3 get_oprec.py` + echo "oprec $oprec" >> $1 +fi + + +#-------------------------------------------- +# check vdelta in deepks +#-------------------------------------------- +if ! test -z "$deepks_v_delta" && [ $deepks_v_delta == 1 ] && [ $deepks_out_labels == 1 ]; then + totalh=`python3 get_sum_numpy.py OUT.autotest/deepks_htot.npy ` + echo "totalh $totalh" >>$1 + totalvdelta=`python3 get_v_delta.py` + echo "totalvdelta $totalvdelta" >>$1 + totalvdp=`python3 get_sum_numpy.py OUT.autotest/deepks_vdpre.npy ` + echo "totalvdp $totalvdp" >> $1 +fi + +#-------------------------------------------- +# check vdelta in deepks +#-------------------------------------------- +if ! test -z "$deepks_v_delta" && [ $deepks_v_delta == 2 ] && [ $deepks_out_labels == 1 ]; then + totalh=`python3 get_sum_numpy.py OUT.autotest/deepks_htot.npy ` + echo "totalh $totalh" >>$1 + totalvdelta=`python3 get_v_delta.py` + echo "totalvdelta $totalvdelta" >>$1 + total_phialpha=`python3 get_sum_numpy.py OUT.autotest/deepks_phialpha.npy ` + echo "total_phialpha $total_phialpha" >> $1 + total_gevdm=`python3 get_sum_numpy.py OUT.autotest/deepks_gevdm.npy ` + echo "total_gevdm $total_gevdm" >> $1 +fi \ No newline at end of file diff --git a/tests/integrate/tools/catch_properties.sh b/tests/integrate/tools/catch_properties.sh index d534b4c0ca..179d713c07 100755 --- a/tests/integrate/tools/catch_properties.sh +++ b/tests/integrate/tools/catch_properties.sh @@ -45,10 +45,6 @@ has_xc=$(get_input_key_value "out_mat_xc" "INPUT") has_xc2=$(get_input_key_value "out_mat_xc2" "INPUT") has_eband_separate=$(get_input_key_value "out_eband_terms" "INPUT") has_r=$(get_input_key_value "out_mat_r" "INPUT") -deepks_out_labels=$(get_input_key_value "deepks_out_labels" "INPUT") -deepks_scf=$(get_input_key_value "deepks_scf" "INPUT") -deepks_bandgap=$(get_input_key_value "deepks_bandgap" "INPUT") -deepks_v_delta=$(get_input_key_value "deepks_v_delta" "INPUT") has_lowf=$(get_input_key_value "out_wfc_lcao" "INPUT") out_app_flag=$(get_input_key_value "out_app_flag" "INPUT") has_wfc_r=$(get_input_key_value "out_wfc_r" "INPUT") @@ -577,61 +573,8 @@ fi #-------------------------------------------- # deepks #-------------------------------------------- -if ! test -z "$deepks_out_labels" && [ $deepks_out_labels == 1 ]; then - sed '/n_des/d' OUT.autotest/deepks_desc.dat > des_tmp.txt - total_des=`sum_file des_tmp.txt 5` - rm des_tmp.txt - echo "totaldes $total_des" >>$1 - if ! test -z "$deepks_scf" && [ $deepks_scf == 1 ]; then - deepks_e_dm=`python3 get_dm_eig.py` - echo "deepks_e_dm $deepks_e_dm" >>$1 - fi - if ! test -z "$has_force" && [ $has_force == 1 ]; then - deepks_f_label=`python3 get_grad_vx.py` - echo "deepks_f_label $deepks_f_label" >>$1 - fi - if ! test -z "$has_stress" && [ $has_stress == 1 ]; then - deepks_s_label=`python3 get_grad_vepsl.py` - echo "deepks_s_label $deepks_s_label" >>$1 - fi -fi - -#-------------------------------------------- -# band gap information -#-------------------------------------------- -if ! test -z "$deepks_bandgap" && [ $deepks_bandgap == 1 ]; then - odelta=`python3 get_odelta.py` - echo "odelta $odelta" >>$1 - oprec=`python3 get_oprec.py` - echo "oprec $oprec" >> $1 -fi - - -#-------------------------------------------- -# check vdelta in deepks -#-------------------------------------------- -if ! test -z "$deepks_v_delta" && [ $deepks_v_delta == 1 ]; then - totalh=`python3 get_sum_numpy.py OUT.autotest/deepks_htot.npy ` - echo "totalh $totalh" >>$1 - totalvdelta=`python3 get_v_delta.py` - echo "totalvdelta $totalvdelta" >>$1 - totalvdp=`python3 get_sum_numpy.py OUT.autotest/deepks_vdpre.npy ` - echo "totalvdp $totalvdp" >> $1 -fi - -#-------------------------------------------- -# check vdelta in deepks -#-------------------------------------------- -if ! test -z "$deepks_v_delta" && [ $deepks_v_delta == 2 ]; then - totalh=`python3 get_sum_numpy.py OUT.autotest/deepks_htot.npy ` - echo "totalh $totalh" >>$1 - totalvdelta=`python3 get_v_delta.py` - echo "totalvdelta $totalvdelta" >>$1 - total_phialpha=`python3 get_sum_numpy.py OUT.autotest/deepks_phialpha.npy ` - echo "total_phialpha $total_phialpha" >> $1 - total_gevdm=`python3 get_sum_numpy.py OUT.autotest/deepks_gevdm.npy ` - echo "total_gevdm $total_gevdm" >> $1 -fi +script_dir=$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd) +bash ${script_dir}/catch_deepks_properties.sh $1 #-------------------------------------------- # check symmetry