diff --git a/source/module_ri/Exx_LRI.hpp b/source/module_ri/Exx_LRI.hpp index 3cba82a0b5..cbe0f91027 100644 --- a/source/module_ri/Exx_LRI.hpp +++ b/source/module_ri/Exx_LRI.hpp @@ -58,10 +58,11 @@ void Exx_LRI::init(const MPI_Comm &mpi_comm_in, const K_Vectors &kv_in) const std::vector>> abfs_same_atom = Exx_Abfs::Construct_Orbs::abfs_same_atom( this->lcaos, this->info.kmesh_times, this->info.pca_threshold ); - if(this->info.files_abfs.empty()) + if(this->info.files_abfs.empty()) { this->abfs = abfs_same_atom; - else + } else { this->abfs = Exx_Abfs::IO::construct_abfs( abfs_same_atom, GlobalC::ORB, this->info.files_abfs, this->info.kmesh_times ); +} Exx_Abfs::Construct_Orbs::print_orbs_size(this->abfs, GlobalV::ofs_running); auto get_ccp_parameter = [this]() -> std::map @@ -86,8 +87,9 @@ void Exx_LRI::init(const MPI_Comm &mpi_comm_in, const K_Vectors &kv_in) this->abfs_ccp = Conv_Coulomb_Pot_K::cal_orbs_ccp(this->abfs, this->info.ccp_type, get_ccp_parameter(), this->info.ccp_rmesh_times); - for( size_t T=0; T!=this->abfs.size(); ++T ) + for( size_t T=0; T!=this->abfs.size(); ++T ) { GlobalC::exx_info.info_ri.abfs_Lmax = std::max( GlobalC::exx_info.info_ri.abfs_Lmax, static_cast(this->abfs[T].size())-1 ); +} this->cv.set_orbitals( this->lcaos, this->abfs, this->abfs_ccp, @@ -108,11 +110,13 @@ void Exx_LRI::cal_exx_ions() // this->m_abfslcaos_lcaos.init_radial_table(Rradial); std::vector atoms(GlobalC::ucell.nat); - for(int iat=0; iat atoms_pos; - for(int iat=0; iat latvec = {RI_Util::Vector3_to_array3(GlobalC::ucell.a1), RI_Util::Vector3_to_array3(GlobalC::ucell.a2), @@ -269,9 +273,11 @@ void Exx_LRI::cal_exx_force() for(int is=0; isexx_lri.cal_force({"","",std::to_string(is),"",""}); - for(std::size_t idim=0; idimexx_lri.force[idim]) + for(std::size_t idim=0; idimexx_lri.force[idim]) { this->force_exx(force_item.first, idim) += std::real(force_item.second); +} +} } const double SPIN_multiple = std::map{{1,2}, {2,1}, {4,1}}.at(GlobalV::NSPIN); // why? @@ -291,9 +297,11 @@ void Exx_LRI::cal_exx_stress() for(int is=0; isexx_lri.cal_stress({"","",std::to_string(is),"",""}); - for(std::size_t idim0=0; idim0stress_exx(idim0,idim1) += std::real(this->exx_lri.stress(idim0,idim1)); +} +} } const double SPIN_multiple = std::map{{1,2}, {2,1}, {4,1}}.at(GlobalV::NSPIN); // why? @@ -310,8 +318,9 @@ std::vector> Exx_LRI::get_abfs_nchis() const for (const auto& abfs_T : this->abfs) { std::vector abfs_nchi_T; - for (const auto& abfs_L : abfs_T) + for (const auto& abfs_L : abfs_T) { abfs_nchi_T.push_back(abfs_L.size()); +} abfs_nchis.push_back(abfs_nchi_T); } return abfs_nchis; diff --git a/source/module_ri/module_exx_symmetry/irreducible_sector.h b/source/module_ri/module_exx_symmetry/irreducible_sector.h index c12d04dc5e..1b2e18d626 100644 --- a/source/module_ri/module_exx_symmetry/irreducible_sector.h +++ b/source/module_ri/module_exx_symmetry/irreducible_sector.h @@ -22,18 +22,20 @@ namespace ModuleSymmetry { bool operator()(const Tap& lhs, const Tap& rhs) const { - if (lhs.first < rhs.first)return true; - else if (lhs.first > rhs.first)return false; - else return lhs.second < rhs.second; + if (lhs.first < rhs.first) {return true; + } else if (lhs.first > rhs.first) {return false; + } else { return lhs.second < rhs.second; +} } }; struct apR_less_func { bool operator()(const TapR& lhs, const TapR& rhs) const { - if (lhs.first < rhs.first)return true; - else if (lhs.first > rhs.first)return false; - else return lhs.second < rhs.second; + if (lhs.first < rhs.first) {return true; + } else if (lhs.first > rhs.first) {return false; + } else { return lhs.second < rhs.second; +} } }; struct len_less_func @@ -44,9 +46,10 @@ namespace ModuleSymmetry } bool operator()(const TC& lhs, const TC& rhs) const { - if (norm2(lhs) < norm2(rhs))return true; - else if (norm2(lhs) > norm2(rhs))return false; - else return lhs < rhs; + if (norm2(lhs) < norm2(rhs)) {return true; + } else if (norm2(lhs) > norm2(rhs)) {return false; + } else { return lhs < rhs; +} } }; diff --git a/source/module_ri/module_exx_symmetry/irreducible_sector_bvk.cpp b/source/module_ri/module_exx_symmetry/irreducible_sector_bvk.cpp index 4ff7344adb..2ed86cb802 100644 --- a/source/module_ri/module_exx_symmetry/irreducible_sector_bvk.cpp +++ b/source/module_ri/module_exx_symmetry/irreducible_sector_bvk.cpp @@ -45,7 +45,8 @@ namespace ModuleSymmetry { //the BvK supercell has the same symmetry as the original cell this->bvk_nsym_ = symm.nrotk; this->isymbvk_to_isym_.resize(symm.nrotk); - for (int isym = 0;isym < symm.nrotk;++isym) this->isymbvk_to_isym_[isym] = isym; + for (int isym = 0;isym < symm.nrotk;++isym) { this->isymbvk_to_isym_[isym] = isym; +} return; } @@ -59,7 +60,8 @@ namespace ModuleSymmetry for (int it = 0;it < st.ntype;++it) { bvk_na[it] = atoms[it].na * bvk_min_period[0] * bvk_min_period[1] * bvk_min_period[2]; - if (it > 0) bvk_istart[it] = bvk_istart[it - 1] + bvk_na[it - 1]; + if (it > 0) { bvk_istart[it] = bvk_istart[it - 1] + bvk_na[it - 1]; +} if (bvk_na[it] < bvk_na[bvk_itmin_type]) { bvk_itmin_type = it; @@ -76,10 +78,10 @@ namespace ModuleSymmetry s3 = a3 = lat.a3 * static_cast(bvk_min_period[2]); ModuleBase::Matrix3 bvk_min_lat = set_matrix3(s1, s2, s3); int at = 0; - for (int it = 0; it < st.ntype; ++it) - for (int c1 = 0;c1 < bvk_min_period[0];++c1) - for (int c2 = 0;c2 < bvk_min_period[1];++c2) - for (int c3 = 0;c3 < bvk_min_period[2];++c3) + for (int it = 0; it < st.ntype; ++it) { + for (int c1 = 0;c1 < bvk_min_period[0];++c1) { + for (int c2 = 0;c2 < bvk_min_period[1];++c2) { + for (int c3 = 0;c3 < bvk_min_period[2];++c3) { for (int ia = 0; ia < atoms[it].na; ++ia) { bvk_dpos[3 * at] = (static_cast (c1) + atoms[it].taud[ia].x) / static_cast(bvk_min_period[0]); @@ -92,6 +94,10 @@ namespace ModuleSymmetry } ++at; } +} +} +} +} // analyze bravis and generate optimized lattice for minimal BvK lattice double cel_const[6], pre_const[6]; @@ -157,4 +163,4 @@ namespace ModuleSymmetry // return in_plain; // } -} \ No newline at end of file +}; \ No newline at end of file diff --git a/source/module_ri/module_exx_symmetry/symmetry_rotation.cpp b/source/module_ri/module_exx_symmetry/symmetry_rotation.cpp index 6a78552a32..62661e799d 100644 --- a/source/module_ri/module_exx_symmetry/symmetry_rotation.cpp +++ b/source/module_ri/module_exx_symmetry/symmetry_rotation.cpp @@ -16,7 +16,8 @@ namespace ModuleSymmetry this->reduce_Cs_ = true; this->abfs_l_nchi_ = abfs_l_nchi; this->abfs_Lmax_ = 0; - for (auto& abfs_T : abfs_l_nchi) this->abfs_Lmax_ = std::max(this->abfs_Lmax_, static_cast(abfs_T.size()) - 1); + for (auto& abfs_T : abfs_l_nchi) { this->abfs_Lmax_ = std::max(this->abfs_Lmax_, static_cast(abfs_T.size()) - 1); +} } void Symmetry_rotation::cal_Ms(const K_Vectors& kv, //const std::vector>& kstars, @@ -34,7 +35,8 @@ namespace ModuleSymmetry } // 1. calculate the rotation matrix in real spherical harmonics representation for each symmetry operation: [T_l (isym)]_mm' std::vector gmatc(nsym_); - for (int i = 0;i < nsym_;++i) gmatc[i] = this->irs_.direct_to_cartesian(ucell.symm.gmatrix[i], ucell.latvec); + for (int i = 0;i < nsym_;++i) { gmatc[i] = this->irs_.direct_to_cartesian(ucell.symm.gmatrix[i], ucell.latvec); +} this->cal_rotmat_Slm(gmatc.data(), reduce_Cs_ ? std::max(this->abfs_Lmax_, ucell.lmax) : ucell.lmax); // 2. calculate the rotation matrix in AO-representation for each ibz_kpoint and symmetry operation: M(k, isym) @@ -44,9 +46,12 @@ namespace ModuleSymmetry kvec_res.x = fmod(kvec.x + 100.5 - 0.5 * symm_prec, 1) - 0.5 + 0.5 * symm_prec; kvec_res.y = fmod(kvec.y + 100.5 - 0.5 * symm_prec, 1) - 0.5 + 0.5 * symm_prec; kvec_res.z = fmod(kvec.z + 100.5 - 0.5 * symm_prec, 1) - 0.5 + 0.5 * symm_prec; - if (std::abs(kvec_res.x) < symm_prec) kvec_res.x = 0.0; - if (std::abs(kvec_res.y) < symm_prec) kvec_res.y = 0.0; - if (std::abs(kvec_res.z) < symm_prec) kvec_res.z = 0.0; + if (std::abs(kvec_res.x) < symm_prec) { kvec_res.x = 0.0; +} + if (std::abs(kvec_res.y) < symm_prec) { kvec_res.y = 0.0; +} + if (std::abs(kvec_res.z) < symm_prec) { kvec_res.z = 0.0; +} return kvec_res; }; int nks_ibz = kv.kstars.size(); // kv.nks = 2 * kv.nks_ibz when nspin=2 @@ -54,9 +59,11 @@ namespace ModuleSymmetry for (int ik_ibz = 0;ik_ibz < nks_ibz;++ik_ibz) { // const TCdouble& kvec_d_ibz = restrict_kpt((*kstars[ik_ibz].begin()).second * ucell.symm.kgmatrix[(*kstars[ik_ibz].begin()).first], ucell.symm.epsilon); - for (auto& isym_kvd : kv.kstars[ik_ibz]) - if (isym_kvd.first < nsym_) + for (auto& isym_kvd : kv.kstars[ik_ibz]) { + if (isym_kvd.first < nsym_) { this->Ms_[ik_ibz][isym_kvd.first] = this->contruct_2d_rot_mat_ao(ucell.symm, ucell.atoms, ucell.st, kv.kvec_d[ik_ibz], isym_kvd.first, pv); +} +} } // output Ms of isym=1 @@ -87,29 +94,35 @@ namespace ModuleSymmetry auto vec_conj = [](const std::vector>& z, const double scal = 1.0) -> std::vector> { std::vector> z_conj(z.size()); - for (int i = 0;i < z.size();++i) z_conj[i] = std::conj(z[i]) * scal; + for (int i = 0;i < z.size();++i) { z_conj[i] = std::conj(z[i]) * scal; +} return z_conj; }; std::vector>> dm_k_full; int nspin0 = GlobalV::NSPIN == 2 ? 2 : 1; dm_k_full.reserve(kv.get_nkstot_full() * nspin0); //nkstot_full didn't doubled by spin int nk = kv.get_nkstot() / nspin0; - for (int is = 0;is < nspin0;++is) - for (int ik_ibz = 0;ik_ibz < nk;++ik_ibz) - for (auto& isym_kvd : kv.kstars[ik_ibz]) + for (int is = 0;is < nspin0;++is) { + for (int ik_ibz = 0;ik_ibz < nk;++ik_ibz) { + for (auto& isym_kvd : kv.kstars[ik_ibz]) { if (isym_kvd.first == 0) { double factor = 1.0 / static_cast(kv.kstars[ik_ibz].size()); std::vector> dm_scaled(pv.get_local_size()); - for (int i = 0;i < pv.get_local_size();++i) dm_scaled[i] = factor * dm_k_ibz[ik_ibz + is * nk][i]; + for (int i = 0;i < pv.get_local_size();++i) { dm_scaled[i] = factor * dm_k_ibz[ik_ibz + is * nk][i]; +} dm_k_full.push_back(dm_scaled); } - else if (vec3_eq(isym_kvd.second, -kv.kvec_d[ik_ibz], this->eps_) && this->TRS_first_) + else if (vec3_eq(isym_kvd.second, -kv.kvec_d[ik_ibz], this->eps_) && this->TRS_first_) { dm_k_full.push_back(vec_conj(dm_k_ibz[ik_ibz + is * nk], 1.0 / static_cast(kv.kstars[ik_ibz].size()))); - else if (isym_kvd.first < nsym_) //space group operations + } else if (isym_kvd.first < nsym_) { //space group operations dm_k_full.push_back(this->rot_matrix_ao(dm_k_ibz[ik_ibz + is * nk], ik_ibz, kv.kstars[ik_ibz].size(), isym_kvd.first, pv)); - else // TRS*spacegroup operations + } else { // TRS*spacegroup operations dm_k_full.push_back(this->rot_matrix_ao(dm_k_ibz[ik_ibz + is * nk], ik_ibz, kv.kstars[ik_ibz].size(), isym_kvd.first - nsym_, pv, true)); +} +} +} +} // test for output @@ -149,14 +162,16 @@ namespace ModuleSymmetry { auto factorial = [](int n) -> int { int result = 1; - for (int i = 1;i <= n;++i) result *= i; + for (int i = 1;i <= n;++i) { result *= i; +} return result; }; double result = 0.0; - for (int i = std::max(0, m2 - m1);i <= std::min(l - m1, l + m2);++i) + for (int i = std::max(0, m2 - m1);i <= std::min(l - m1, l + m2);++i) { result += std::pow(-1, i) * std::sqrt(factorial(l + m1) * factorial(l - m1) * factorial(l + m2) * factorial(l - m2)) * std::pow(std::cos(beta / 2), 2 * l + m2 - m1 - 2 * i) * std::pow(-std::sin(beta / 2), m1 - m2 + 2 * i) / (factorial(i) * factorial(l - m1 - i) * factorial(l + m2 - i) * factorial(i - m2 + m1)); +} return result; } @@ -173,14 +188,19 @@ namespace ModuleSymmetry { if (m1 == m2) { - if (m1 == 0) return 1.0; - if (m1 > 0) return 1 / std::sqrt(2); - if (m1 < 0) return std::pow(-1, m1) * ModuleBase::IMAG_UNIT / std::sqrt(2); + if (m1 == 0) { return 1.0; +} + if (m1 > 0) { return 1 / std::sqrt(2); +} + if (m1 < 0) { return std::pow(-1, m1) * ModuleBase::IMAG_UNIT / std::sqrt(2); +} } else if (m1 == -m2) { - if (m1 > 0) return -ModuleBase::IMAG_UNIT / std::sqrt(2); - if (m1 < 0) return std::pow(-1, m1) / std::sqrt(2); + if (m1 > 0) { return -ModuleBase::IMAG_UNIT / std::sqrt(2); +} + if (m1 < 0) { return std::pow(-1, m1) / std::sqrt(2); +} } return 0.0; } @@ -197,19 +217,23 @@ namespace ModuleSymmetry { // use the 2-angle elements to get alpha and gamma alpha = std::atan2(gmatc.e32, gmatc.e31); - if (alpha < 0) alpha += 2 * ModuleBase::PI; + if (alpha < 0) { alpha += 2 * ModuleBase::PI; +} gamma = std::atan2(gmatc.e23, -gmatc.e13); - if (gamma < 0) gamma += 2 * ModuleBase::PI; + if (gamma < 0) { gamma += 2 * ModuleBase::PI; +} // use the larger one of 2-angle elements to calculate beta - if (std::fabs(gmatc.e32) > std::fabs(gmatc.e31)) + if (std::fabs(gmatc.e32) > std::fabs(gmatc.e31)) { beta = std::atan2(gmatc.e32 / std::sin(alpha), gmatc.e33); - else + } else { beta = std::atan2(gmatc.e31 / std::cos(alpha), gmatc.e33); +} } else {//sin(beta)=0, beta = 0 or pi, only (alpha+gamma) or (alpha-gamma) is important. now assign this to alpha. alpha = std::atan2(gmatc.e12, gmatc.e11); - if (alpha < 0) alpha += 2 * ModuleBase::PI; + if (alpha < 0) { alpha += 2 * ModuleBase::PI; +} // if beta=0, gmatc.e11=cos(alpha+gamma), gmatc.e21=sin(alpha+gamma) // if beta=pi, gmatc.e11=cos(pi+alpha-gamma), gmatc.e21=sin(pi+alpha-gamma) if (gmatc.e33 > 0) @@ -239,22 +263,29 @@ namespace ModuleSymmetry auto set_integer = [](RI::Tensor>& mat) -> void { double zero_thres = 1e-10; - for (int i = 0;i < mat.shape[0];++i) + for (int i = 0;i < mat.shape[0];++i) { for (int j = 0;j < mat.shape[1];++j) { - if (std::abs(mat(i, j).real() - std::round(mat(i, j).real())) < zero_thres) mat(i, j).real(std::round(mat(i, j).real())); - if (std::abs(mat(i, j).imag() - std::round(mat(i, j).imag())) < zero_thres) mat(i, j).imag(std::round(mat(i, j).imag())); + if (std::abs(mat(i, j).real() - std::round(mat(i, j).real())) < zero_thres) { mat(i, j).real(std::round(mat(i, j).real())); +} + if (std::abs(mat(i, j).imag() - std::round(mat(i, j).imag())) < zero_thres) { mat(i, j).imag(std::round(mat(i, j).imag())); +} } +} }; this->rotmat_Slm_.resize(nsym_); // c matrix is independent on isym std::vector>> c_mm(lmax + 1); - for (int l = 0;l <= lmax;++l) + for (int l = 0;l <= lmax;++l) { c_mm[l] = RI::Tensor>({ size_t(2 * l + 1), size_t(2 * l + 1) }); - for (int l = 0;l <= lmax;++l) - for (int m1 = -l;m1 <= l;++m1) - for (int m2 = -l;m2 <= l;++m2) +} + for (int l = 0;l <= lmax;++l) { + for (int m1 = -l;m1 <= l;++m1) { + for (int m2 = -l;m2 <= l;++m2) { c_mm[l](m2im(m1), m2im(m2)) = ovlp_Ylm_Slm(l, m1, m2); +} +} +} for (int isym = 0;isym < nsym_;++isym) { @@ -266,9 +297,11 @@ namespace ModuleSymmetry for (int l = 0;l <= lmax;++l) {// wigner D matrix RI::Tensor> D_mm({ size_t(2 * l + 1), size_t(2 * l + 1) }); - for (int m1 = -l;m1 <= l;++m1) - for (int m2 = -l;m2 <= l;++m2) + for (int m1 = -l;m1 <= l;++m1) { + for (int m2 = -l;m2 <= l;++m2) { D_mm(m2im(m1), m2im(m2)) = wigner_D(euler_angle, l, m1, m2, (gmatc[isym].Det() < 0)); +} +} this->rotmat_Slm_[isym][l] = c_mm[l].dagger() * D_mm * c_mm[l]; // set_integer(this->rotmat_Slm_[isym][l]); } @@ -314,25 +347,29 @@ namespace ModuleSymmetry void Symmetry_rotation::set_block_to_mat2d(const int starti, const int startj, const RI::Tensor>& block, std::vector>& obj_mat, const Parallel_2D& pv, const bool trans) const { // caution: ComplaxMatrix is row-major(col-continuous), but obj_mat is col-major(row-continuous) - for (int j = 0;j < block.shape[0];++j)//outside dimension - for (int i = 0;i < block.shape[1];++i) //inside dimension + for (int j = 0;j < block.shape[0];++j) {//outside dimension + for (int i = 0;i < block.shape[1];++i) { //inside dimension if (pv.in_this_processor(starti + i, startj + j)) { int index = pv.global2local_col(startj + j) * pv.get_row_size() + pv.global2local_row(starti + i); obj_mat[index] = trans ? block(i, j) : block(j, i); } +} +} } void Symmetry_rotation::set_block_to_mat2d(const int starti, const int startj, const RI::Tensor>& block, std::vector& obj_mat, const Parallel_2D& pv, const bool trans) const { // caution: ComplaxMatrix is row-major(col-continuous), but obj_mat is col-major(row-continuous) - for (int j = 0;j < block.shape[0];++j)//outside dimension - for (int i = 0;i < block.shape[1];++i) //inside dimension + for (int j = 0;j < block.shape[0];++j) {//outside dimension + for (int i = 0;i < block.shape[1];++i) { //inside dimension if (pv.in_this_processor(starti + i, startj + j)) { int index = pv.global2local_col(startj + j) * pv.get_row_size() + pv.global2local_row(starti + i); obj_mat[index] = trans ? block(i, j).real() : block(j, i).real(); } +} +} } // 2d-block parallized rotation matrix in AO-representation, denoted as M. @@ -428,14 +465,16 @@ namespace ModuleSymmetry { const ModuleBase::Vector3& R_index = adjs.box[ad]; if (ucell.cal_dtau(iat1, iat2, R_index).norm() * ucell.lat0 - < ucell.atoms[it1].Rcut + ucell.atoms[it2].Rcut) + < ucell.atoms[it1].Rcut + ucell.atoms[it2].Rcut) { Rs_set.insert({ R_index.x, R_index.y, R_index.z }); +} } } } // set to vector std::vector Rs(Rs_set.size()); - for (auto& R : Rs_set) Rs.push_back(R); + for (auto& R : Rs_set) { Rs.push_back(R); +} return Rs; } diff --git a/source/module_ri/module_exx_symmetry/symmetry_rotation_R.hpp b/source/module_ri/module_exx_symmetry/symmetry_rotation_R.hpp index fe54d89479..12d88d1a17 100644 --- a/source/module_ri/module_exx_symmetry/symmetry_rotation_R.hpp +++ b/source/module_ri/module_exx_symmetry/symmetry_rotation_R.hpp @@ -13,8 +13,9 @@ namespace ModuleSymmetry GlobalV::ofs_running << name << ":\n"; for (int i = 0;i < t.shape[0];++i) { - for (int j = 0;j < t.shape[1];++j) + for (int j = 0;j < t.shape[1];++j) { GlobalV::ofs_running << ((std::abs(t(i, j)) > threshold) ? t(i, j) : static_cast(0)) << " "; +} GlobalV::ofs_running << std::endl; } } @@ -27,8 +28,9 @@ namespace ModuleSymmetry GlobalV::ofs_running << "abf: " << a << '\n'; for (int i = 0;i < t.shape[1];++i) { - for (int j = 0;j < t.shape[2];++j) + for (int j = 0;j < t.shape[2];++j) { GlobalV::ofs_running << ((std::abs(t(a, i, j)) > threshold) ? t(a, i, j) : static_cast(0)) << " "; +} GlobalV::ofs_running << std::endl; } GlobalV::ofs_running << std::endl; @@ -115,9 +117,11 @@ namespace ModuleSymmetry inline void set_block(const int starti, const int startj, const RI::Tensor>& block, RI::Tensor& obj_tensor) { // no changing row/col order - for (int i = 0;i < block.shape[0];++i) - for (int j = 0;j < block.shape[1];++j) + for (int i = 0;i < block.shape[0];++i) { + for (int j = 0;j < block.shape[1];++j) { obj_tensor(starti + i, startj + j) = RI::Global_Func::convert(block(i, j)); +} +} } template @@ -180,7 +184,8 @@ namespace ModuleSymmetry RI::Tensor Symmetry_rotation::set_rotation_matrix_abf(const int& type, const int& isym)const { int nabfs = 0; - for (int l = 0;l < this->abfs_l_nchi_[type].size();++l)nabfs += this->abfs_l_nchi_[type][l] * (2 * l + 1); + for (int l = 0;l < this->abfs_l_nchi_[type].size();++l) {nabfs += this->abfs_l_nchi_[type][l] * (2 * l + 1); +} RI::Tensor T({ static_cast(nabfs), static_cast(nabfs) }); // check if zero int iw = 0; for (int L = 0;L < this->abfs_l_nchi_[type].size();++L) @@ -206,9 +211,10 @@ namespace ModuleSymmetry assert(C.shape.size() == 3); const int& slice_size = C.shape[1] * C.shape[2]; // step 1: multiply 2 AOs' rotation matrices - for (int iabf = 0;iabf < C.shape[0];++iabf) + for (int iabf = 0;iabf < C.shape[0];++iabf) { this->rotate_atompair_serial(Cout.ptr() + iabf * slice_size, C.ptr() + iabf * slice_size, a1.nw, a2.nw, isym, a1, a2, 'H'); +} // step 2: multiply the ABFs' rotation matrix from the left const RI::Tensor& Tabfs = this->set_rotation_matrix_abf(type1, isym); RI::Sym::T1_HR(Cout.ptr(), Cout.ptr(), Tabfs, slice_size);