Skip to content

Commit 65db964

Browse files
dzzz2001dzzz2001
authored
Merge branch 'develop' into force_ew_opt
2 parents 390f9da + 94f71a4 commit 65db964

File tree

6 files changed

+50
-23
lines changed

6 files changed

+50
-23
lines changed

‎source/source_cell/k_vector_utils.cpp‎

Lines changed: 8 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -604,8 +604,8 @@ void kvec_ibz_kpoint(K_Vectors& kv,
604604
}
605605
return;
606606
};
607-
// for output in kpoints file
608-
int ibz_index[kv.get_nkstot()];
607+
// update map k -> irreducible k
608+
kv.ibz_index.assign( kv.get_nkstot_full(), -1); // -1 means not in ibz_kpoint list
609609
// search in all k-poins.
610610
for (int i = 0; i < kv.get_nkstot(); ++i)
611611
{
@@ -667,7 +667,7 @@ void kvec_ibz_kpoint(K_Vectors& kv,
667667
// nkstot_ibz indicate the index of ibz kpoint.
668668
kvec_d_ibz[nkstot_ibz] = kv.kvec_d[i];
669669
// output in kpoints file
670-
ibz_index[i] = nkstot_ibz;
670+
kv.ibz_index[i] = nkstot_ibz;
671671

672672
// the weight should be averged k-point weight.
673673
wk_ibz[nkstot_ibz] = weight;
@@ -688,7 +688,7 @@ void kvec_ibz_kpoint(K_Vectors& kv,
688688
double kmol_new = kv.kvec_d[i].norm2();
689689
double kmol_old = kvec_d_ibz[exist_number].norm2();
690690

691-
ibz_index[i] = exist_number;
691+
kv.ibz_index[i] = exist_number;
692692

693693
// std::cout << "\n kmol_new = " << kmol_new;
694694
// std::cout << "\n kmol_old = " << kmol_old;
@@ -765,10 +765,10 @@ void kvec_ibz_kpoint(K_Vectors& kv,
765765
kv.kvec_d[i].x,
766766
kv.kvec_d[i].y,
767767
kv.kvec_d[i].z,
768-
ibz_index[i] + 1,
769-
kvec_d_ibz[ibz_index[i]].x,
770-
kvec_d_ibz[ibz_index[i]].y,
771-
kvec_d_ibz[ibz_index[i]].z);
768+
kv.ibz_index[i] + 1,
769+
kvec_d_ibz[kv.ibz_index[i]].x,
770+
kvec_d_ibz[kv.ibz_index[i]].y,
771+
kvec_d_ibz[kv.ibz_index[i]].z);
772772
}
773773
ss << table << std::endl;
774774
skpt = ss.str();

‎source/source_cell/klist.cpp‎

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -160,6 +160,13 @@ void K_Vectors::set(const UnitCell& ucell,
160160
// set the k vectors for the up and down spin
161161
this->set_kup_and_kdw();
162162

163+
// initialize ibz_index
164+
this->ibz_index.resize(this->nkstot_full);
165+
for (int ik = 0; ik < this->nkstot_full; ik++)
166+
{
167+
this->ibz_index[ik] = ik;
168+
}
169+
163170
// get ik2iktot
164171
this->cal_ik_global();
165172

‎source/source_cell/klist.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -126,6 +126,7 @@ class K_Vectors
126126
}
127127

128128
std::vector<int> ik2iktot; ///<[nks] map ik to the global index of k points
129+
std::vector<int> ibz_index; ///< map k points (before symmetry reduction) to irreducible k-points
129130

130131
/**
131132
* @brief Updates the k-points to use the irreducible Brillouin zone (IBZ).

‎source/source_cell/pseudo.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,7 @@
11
#ifndef PSEUDO_H
22
#define PSEUDO_H
33

4+
#include <cstdint>
45
#include "source_base/global_function.h"
56
#include "source_io/output.h"
67

‎source/source_lcao/module_ri/RPA_LRI.h‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -49,7 +49,7 @@ template <typename T, typename Tdata> class RPA_LRI
4949
const psi::Psi<T>& psi,
5050
const elecstate::ElecState* pelec);
5151
void out_eigen_vector(const Parallel_Orbitals& parav, const psi::Psi<T>& psi);
52-
void out_struc(const ModuleBase::Matrix3& latvec, const ModuleBase::Matrix3& G);
52+
void out_struc(const UnitCell &ucell);
5353
void out_bands(const elecstate::ElecState *pelec);
5454

5555
void out_Cs(const UnitCell &ucell);

‎source/source_lcao/module_ri/RPA_LRI.hpp‎

Lines changed: 32 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -69,6 +69,7 @@ void RPA_LRI<T, Tdata>::cal_rpa_cv(const UnitCell& ucell)
6969
});
7070
std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>& Cs = std::get<0>(Cs_dCs);
7171
this->Cs_period = RI::RI_Tools::cal_period(Cs, period);
72+
this->Cs_period = exx_lri_rpa.exx_lri.post_2D.set_tensors_map2(this->Cs_period);
7273
}
7374

7475
template <typename T, typename Tdata>
@@ -112,6 +113,11 @@ void RPA_LRI<T, Tdata>::cal_postSCF_exx(const elecstate::DensityMatrix<T, Tdata>
112113
GlobalC::exx_info.info_global.ccp_type = Conv_Coulomb_Pot_K::Ccp_Type::Hf;
113114
GlobalC::exx_info.info_global.hybrid_alpha = 1;
114115
GlobalC::exx_info.info_ri.ccp_rmesh_times = PARAM.inp.rpa_ccp_rmesh_times;
116+
GlobalC::exx_info.info_global.coulomb_param[Conv_Coulomb_Pot_K::Coulomb_Type::Fock].resize(1);
117+
GlobalC::exx_info.info_global.coulomb_param[Conv_Coulomb_Pot_K::Coulomb_Type::Fock] = {{
118+
{"alpha", "1"},
119+
{"singularity_correction", "spencer"} }};
120+
115121

116122
exx_lri_rpa.init(mpi_comm_in, ucell, kv, orb);
117123
exx_lri_rpa.cal_exx_ions(ucell,PARAM.inp.out_ri_cv);
@@ -120,8 +126,7 @@ void RPA_LRI<T, Tdata>::cal_postSCF_exx(const elecstate::DensityMatrix<T, Tdata>
120126
exx_lri_rpa.cal_exx_elec(Ds, ucell,*dm.get_paraV_pointer(), &symrot);
121127
} else {
122128
exx_lri_rpa.cal_exx_elec(Ds, ucell,*dm.get_paraV_pointer());
123-
}
124-
// cout<<"postSCF_Eexx: "<<exx_lri_rpa.Eexx<<endl;
129+
}
125130
}
126131

127132
template <typename T, typename Tdata>
@@ -133,7 +138,7 @@ void RPA_LRI<T, Tdata>::out_for_RPA(const UnitCell& ucell,
133138
ModuleBase::TITLE("DFT_RPA_interface", "out_for_RPA");
134139
this->out_bands(pelec);
135140
this->out_eigen_vector(parav, psi);
136-
this->out_struc(ucell.latvec, ucell.G);
141+
this->out_struc(ucell);
137142

138143
this->cal_rpa_cv(ucell);
139144
std::cout << "rpa_pca_threshold: " << this->info.pca_threshold << std::endl;
@@ -217,37 +222,50 @@ void RPA_LRI<T, Tdata>::out_eigen_vector(const Parallel_Orbitals& parav, const p
217222
}
218223

219224
template <typename T, typename Tdata>
220-
void RPA_LRI<T, Tdata>::out_struc(const ModuleBase::Matrix3& latvec, const ModuleBase::Matrix3& G)
225+
void RPA_LRI<T, Tdata>::out_struc(const UnitCell &ucell)
221226
{
222227
if (GlobalV::MY_RANK != 0)
223228
{
224229
return;
225230
}
226231
ModuleBase::TITLE("DFT_RPA_interface", "out_struc");
227232
double TWOPI_Bohr2A = ModuleBase::TWO_PI * ModuleBase::BOHR_TO_A;
228-
const int nks_tot = PARAM.inp.nspin == 2 ? (int)p_kv->get_nks() / 2 : p_kv->get_nks();
229-
ModuleBase::Matrix3 lat = latvec / ModuleBase::BOHR_TO_A;
230-
ModuleBase::Matrix3 G_RPA = G * TWOPI_Bohr2A;
233+
const int nks_tot = PARAM.inp.nspin == 2 ? (int)p_kv->get_nkstot() / 2 : p_kv->get_nkstot();
234+
const int nks_tot_full = p_kv->get_nkstot_full();
235+
const int natom = ucell.nat;
236+
ModuleBase::Matrix3 lat = ucell.latvec / ModuleBase::BOHR_TO_A;
237+
ModuleBase::Matrix3 G_RPA = ucell.G * TWOPI_Bohr2A;
231238
std::stringstream ss;
232239
ss << "stru_out";
233240
std::ofstream ofs;
234241
ofs.open(ss.str().c_str(), std::ios::out);
235-
ofs << lat.e11 << std::setw(15) << lat.e12 << std::setw(15) << lat.e13 << std::endl;
242+
ofs << std::fixed << std::setprecision(9) << lat.e11 << std::setw(15) << lat.e12 << std::setw(15) << lat.e13 << std::endl;
236243
ofs << lat.e21 << std::setw(15) << lat.e22 << std::setw(15) << lat.e23 << std::endl;
237244
ofs << lat.e31 << std::setw(15) << lat.e32 << std::setw(15) << lat.e33 << std::endl;
238245

239246
ofs << G_RPA.e11 << std::setw(15) << G_RPA.e12 << std::setw(15) << G_RPA.e13 << std::endl;
240247
ofs << G_RPA.e21 << std::setw(15) << G_RPA.e22 << std::setw(15) << G_RPA.e23 << std::endl;
241248
ofs << G_RPA.e31 << std::setw(15) << G_RPA.e32 << std::setw(15) << G_RPA.e33 << std::endl;
242249

250+
ofs << natom << std::endl;
251+
for(int iat=0; iat < natom; ++iat) {
252+
int it = ucell.iat2it[iat];
253+
ModuleBase::Vector3<double> atom_pos = ucell.atoms[ it ].tau[ ucell.iat2ia[iat] ]/ ModuleBase::BOHR_TO_A;
254+
ofs << atom_pos.x << std::setw(15) << atom_pos.y << std::setw(15) << atom_pos.z << std::setw(15) << (it + 1) << std::endl;
255+
}
256+
243257
ofs << p_kv->nmp[0] << std::setw(6) << p_kv->nmp[1] << std::setw(6) << p_kv->nmp[2] << std::setw(6) << std::endl;
244258

245259
for (int ik = 0; ik != nks_tot; ik++)
246260
{
247-
ofs << std::setw(15) << std::fixed << std::setprecision(9) << p_kv->kvec_c[ik].x * TWOPI_Bohr2A << std::setw(15)
248-
<< std::fixed << std::setprecision(9) << p_kv->kvec_c[ik].y * TWOPI_Bohr2A << std::setw(15) << std::fixed
249-
<< std::setprecision(9) << p_kv->kvec_c[ik].z * TWOPI_Bohr2A << std::endl;
261+
ofs << std::setw(15) << p_kv->kvec_c[ik].x * TWOPI_Bohr2A << std::setw(15)
262+
<< p_kv->kvec_c[ik].y * TWOPI_Bohr2A << std::setw(15) << p_kv->kvec_c[ik].z * TWOPI_Bohr2A << std::endl;
263+
}
264+
265+
for (int ik = 0; ik != nks_tot_full; ++ik){
266+
ofs << (p_kv->ibz_index[ik] + 1) << std::endl;
250267
}
268+
251269
ofs.close();
252270
return;
253271
}
@@ -381,11 +399,11 @@ void RPA_LRI<T, Tdata>::out_coulomb_k(const UnitCell &ucell)
381399
size_t nu_num = exx_lri_rpa.exx_objs[Conv_Coulomb_Pot_K::Coulomb_Method::Center2].cv.get_index_abfs_size(ucell.iat2it[iJ]);
382400
ofs << all_mu << " " << mu_shift[I] + 1 << " " << mu_shift[I] + mu_num << " " << mu_shift[iJ] + 1
383401
<< " " << mu_shift[iJ] + nu_num << std::endl;
384-
ofs << ik + 1 << " " << p_kv->wk[ik] / 2.0 * PARAM.inp.nspin << std::endl;
402+
ofs << ik + 1 << " " << std::fixed << std::setprecision(12) << p_kv->wk[ik] / 2.0 * PARAM.inp.nspin << std::endl;
385403
for (int i = 0; i != vq_J.data->size(); i++)
386404
{
387-
ofs << std::setw(21) << std::fixed << std::setprecision(12) << (*vq_J.data)[i].real()
388-
<< std::setw(21) << std::fixed << std::setprecision(12) << (*vq_J.data)[i].imag() << std::endl;
405+
ofs << std::setw(21) << (*vq_J.data)[i].real()
406+
<< std::setw(21) << (*vq_J.data)[i].imag() << std::endl;
389407
}
390408
}
391409
}

0 commit comments

Comments
 (0)