Skip to content

Commit 43a67cb

Browse files
QianruipkuQianruipku
authored
Merge pull request #598 from Qianruipku/planewave
< feature >
2 parents 09d5339 + 97da4c0 commit 43a67cb

File tree

15 files changed

+861
-94
lines changed

15 files changed

+861
-94
lines changed

‎source/module_pw/pw_basis.cpp‎

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,6 @@
11
#include "pw_basis.h"
22
#include "../module_base/mymath.h"
33
#include <iostream>
4-
#include "../module_base/memory.h"
54
#include "../module_base/timer.h"
65
#include "../module_base/global_function.h"
76

‎source/module_pw/pw_basis.h‎

Lines changed: 6 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -100,14 +100,14 @@ class PW_Basis
100100

101101
void collect_local_pw();
102102

103-
void collect_tot_pw(
104-
double* gg_global,
105-
ModuleBase::Vector3<double> *gdirect_global,
106-
ModuleBase::Vector3<double> *gcar_global
107-
);
103+
// void collect_tot_pw(
104+
// double* gg_global,
105+
// ModuleBase::Vector3<double> *gdirect_global,
106+
// ModuleBase::Vector3<double> *gcar_global
107+
// );
108108

109109

110-
private:
110+
public:
111111
bool gamma_only; // only half g are used.
112112
double ggecut; //Energy cut off for g^2/2, unit in 1/lat0^2, ggecut=ecutwfc(Ry)*lat0^2/4pi^2
113113
double lat0; //unit length for lattice, unit in bohr
@@ -208,11 +208,6 @@ class PW_Basis
208208
// void gatherp_scatters_gamma(std::complex<double> *in, std::complex<double> *out); //gather planes and scatter sticks of all processors, used when gamma_only
209209
// void gathers_scatterp_gamma(std::complex<double> *in, std::complex<double> *out); //gather sticks of and scatter planes of all processors, used when gamma only
210210

211-
212-
213-
214-
215-
216211
};
217212

218213
}

‎source/module_pw/pw_basis_k.cpp‎

Lines changed: 104 additions & 46 deletions
Original file line numberDiff line numberDiff line change
@@ -1,52 +1,69 @@
11
#include "pw_basis_k.h"
2+
#include "../module_base/constants.h"
3+
namespace ModulePW
4+
{
5+
26
PW_Basis_K::PW_Basis_K()
37
{
48
nks = 1;
59
kvec_d = NULL;
6-
ngk = NULL;
7-
GR_index = NULL;
10+
kvec_c = NULL;
11+
npwk = NULL;
12+
igl2isz_k = NULL;
813
}
914
PW_Basis_K::~PW_Basis_K()
1015
{
1116
if(kvec_d != NULL) delete[] kvec_d;
12-
if(ngk != NULL) delete[] ngk;
13-
if(GR_index != NULL) delete[] GR_index;
17+
if(kvec_c != NULL) delete[] kvec_c;
18+
if(npwk != NULL) delete[] npwk;
19+
if(igl2isz_k != NULL) delete[] igl2isz_k;
1420
}
1521

1622
void PW_Basis_K:: initparameters(
1723
bool gamma_only_in,
18-
double ecut_in,
1924
double gk_ecut_in,
2025
int nks_in, //number of k points in this pool
2126
ModuleBase::Vector3<double> *kvec_d_in, // Direct coordinates of k points
2227
int poolnproc_in, // Number of processors in this pool
2328
int poolrank_in, // Rank in this pool
24-
int distribution_type_in,
29+
int distribution_type_in
2530
)
2631
{
27-
initparameters(gamma_only_in,ecut_in,poolnproc_in,poolrank_in,distribution_type_in);
2832
this->nks = nks_in;
29-
this->kvec_d = kvec_d_in;
30-
this->gk_ecut = gk_ecut_in;
33+
this->kvec_d = new ModuleBase::Vector3<double> [nks];
34+
this->kvec_c = new ModuleBase::Vector3<double> [nks];
35+
36+
double kmaxmod = 0;
3137
for(int ik = 0 ; ik < this->nks ; ++ik)
3238
{
33-
kvec_c =
39+
this->kvec_d[ik] = kvec_d_in[ik];
40+
this->kvec_c[ik] = this->kvec_d[ik] * this->G;
41+
double kmod = sqrt(this->kvec_c[ik] * this->kvec_c[ik]);
42+
if(kmod > kmaxmod) kmaxmod = kmod;
3443
}
44+
// MPI_Allreduce(MPI_IN_PLACE, &kmaxmod, 1, MPI_DOUBLE, MPI_MAX , MPI_COMM_WORLD);
45+
double tpiba2 = ModuleBase::TWO_PI * ModuleBase::TWO_PI / this->lat0 / this->lat0;
46+
this->gk_ecut = gk_ecut_in/tpiba2;
47+
this->ggecut = pow(sqrt(this->gk_ecut) + kmaxmod, 2);
48+
49+
this->gamma_only = gamma_only_in;
50+
if(kmaxmod > 0) this->gamma_only = false; //if it is not the gamma point, we do not use gamma_only
51+
if (this->gamma_only) this->ny = int(this->bigny / 2) + 1;
52+
else this->ny = bigny;
53+
this->nxy = this->nx * this->ny;
54+
this->nxyz = this->nxy * this->nz;
55+
56+
this->poolnproc = poolnproc_in;
57+
this->poolrank = poolrank_in;
58+
this->distribution_type = distribution_type_in;
3559
return;
3660
}
3761

3862
void PW_Basis_K::setupIndGk()
3963
{
40-
ModuleBase::TITLE("PW_Basis_K","setupIndGk");
41-
ModuleBase::timer::tick("PW_Basis_K","setupIndGk");
42-
43-
igk.create(nks, npw);
44-
ModuleBase::Memory::record("PW_Basis_K","igk",nks*npw,"int");
45-
46-
//=============================================
47-
// Done this in each cpu ( use nks,not nkstot )
48-
// notice : using cartesian coordinate
49-
//=============================================
64+
//count npwk
65+
this->npwk_max = 0;
66+
this->npwk = new int [this->nks];
5067
for (int ik = 0; ik < this->nks; ik++)
5168
{
5269
int ng = 0;
@@ -55,8 +72,7 @@ void PW_Basis_K::setupIndGk()
5572
const double gk2 = this->get_GPlusK_cartesian(ik, ig).norm2();
5673
if (gk2 <= this->gk_ecut)
5774
{
58-
this->igk(ik, ng) = ig;
59-
ng++;
75+
++ng;
6076
}
6177
}
6278
this->npwk[ik] = ng;
@@ -66,30 +82,72 @@ void PW_Basis_K::setupIndGk()
6682
}
6783
}
6884

69-
// bool out_gk =0; //DIY! mohan 2011-10-03
70-
// if(out_gk)
71-
// {
72-
// std::stringstream ss;
73-
// ss << GlobalV::global_out_dir << "PW_GK" << GlobalV::MY_RANK+1 << ".dat";
74-
// std::ofstream ofs( ss.str().c_str() );
75-
// ofs << GlobalC::pw.ggpsi << " (ggpsi, Ry)" << std::endl;
76-
// ofs << GlobalC::pw.ggwfc << " (ggwfc, Ry)" << std::endl;
77-
// ofs << GlobalC::pw.ggwfc2 << " (ggwfc2, Ry)" << std::endl;
78-
// ModuleBase::Vector3<double> f;
79-
// for(int ik=0; ik < nks; ++ik)
80-
// {
81-
// ofs << ik+1 << " (Index of k)" << std::endl;
82-
// ofs << pwb.ngmw << " (Number of plane waves)" << std::endl;
83-
// for(int ig=0; ig < pwb.ngmw; ++ig)
84-
// {
85-
// f = pwb.get_GPlusK_cartesian(ik, ig);
86-
// ofs << f.x << " " << f.y << " " << f.z << " " << f.norm() << std::endl;
87-
// }
88-
// }
89-
// ofs.close();
90-
// }
91-
92-
ModuleBase::timer::tick("PW_Basis_K","setupIndGk");
85+
//get igl2isz_k
86+
this->igl2isz_k = new int [this->nks * this->npwk_max];
87+
for (int ik = 0; ik < this->nks; ik++)
88+
{
89+
int igl = 0;
90+
for (int ig = 0; ig < npw ; ig++)
91+
{
92+
const double gk2 = this->get_GPlusK_cartesian(ik, ig).norm2();
93+
if (gk2 <= this->gk_ecut)
94+
{
95+
this->igl2isz_k[ik*npwk_max + igl] = this->ig2isz[ig];
96+
++igl;
97+
}
98+
}
99+
}
100+
101+
delete[] this->ig2isz;
102+
this->ig2isz = NULL;
93103

94104
return;
95105
}
106+
void PW_Basis_K::setuptransform()
107+
{
108+
this->distribute_r();
109+
this->distribute_g();
110+
this->getstartgr();
111+
this->setupIndGk();
112+
this->ft.initfft(this->nx,this->bigny,this->nz,this->liy,this->riy,this->nst,this->nplane,this->poolnproc,this->gamma_only);
113+
this->ft.setupFFT();
114+
}
115+
116+
void PW_Basis_K::collect_local_pw()
117+
{
118+
if(gg != NULL) delete[] gg;
119+
if(gdirect != NULL) delete[] gdirect;
120+
if(gcar != NULL) delete[] gcar;
121+
this->gg = new double[this->npwk_max * this->nks];
122+
this->gdirect = new ModuleBase::Vector3<double>[this->npwk_max * this->nks];
123+
this->gcar = new ModuleBase::Vector3<double>[this->npwk_max * this->nks];
124+
125+
ModuleBase::Vector3<double> f;
126+
for(int ik = 0 ; ik < this->nks ; ++ik)
127+
{
128+
for(int igl = 0 ; igl < this-> npwk[ik] ; ++igl)
129+
{
130+
int isz = this->igl2isz_k[ik * npwk_max + igl];
131+
int iz = isz % this->nz;
132+
int is = isz / this->nz;
133+
int ixy = this->is2ixy[is];
134+
int ix = ixy / this->ny;
135+
int iy = ixy % this->ny;
136+
if (ix >= int(this->nx/2) + 1) ix -= this->nx;
137+
if (iy >= int(this->bigny/2) + 1) iy -= this->bigny;
138+
if (iz >= int(this->nz/2) + 1) iz -= this->nz;
139+
f.x = ix;
140+
f.y = iy;
141+
f.z = iz;
142+
this->gg[ik * npwk_max + igl] = f * (this->GGT * f);
143+
this->gdirect[ik * npwk_max + igl] = f;
144+
this->gcar[ik * npwk_max + igl] = f * this->G;
145+
}
146+
}
147+
}
148+
149+
150+
151+
152+
153+
}

‎source/module_pw/pw_basis_k.h‎

Lines changed: 23 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,8 @@
22
#define PWBASISK_H
33

44
#include "pw_basis.h"
5-
#include "../module_base/intarray.h"
5+
namespace ModulePW
6+
{
67

78
//
89
//Special pw_basis class.
@@ -15,38 +16,52 @@ class PW_Basis_K : public PW_Basis
1516
public:
1617
PW_Basis_K();
1718
~PW_Basis_K();
19+
1820
void initparameters(
1921
bool gamma_only_in,
2022
double ecut_in,
21-
double gk_ecut_in,
2223
int nk_in, //number of k points in this pool
2324
ModuleBase::Vector3<double> *kvec_d, // Direct coordinates of k points
2425
int poolnproc_in, // Number of processors in this pool
2526
int poolrank_in, // Rank in this pool
2627
int distribution_type_in
2728
);
28-
void setupIndGk(); //set up igk
29-
3029

3130

3231
public:
3332
int nks;//number of k points in this pool
3433
ModuleBase::Vector3<double> *kvec_d; // Direct coordinates of k points
3534
ModuleBase::Vector3<double> *kvec_c; // Cartesian coordinates of k points
36-
ModuleBase::IntArray igk; //[nks, npw_max] map igk_local to ig_local
3735
int *npwk; //[nks] number of plane waves of different k-points
3836
int npwk_max; //max npwk among all nks k-points
3937
double gk_ecut; //Energy cut off for (g+k)^2/2
4038

4139
public:
42-
void init_k();//initialize some data for current k-points
43-
//After inik_k()
44-
int *GR_index; //[npw_max] map igk_local to (is,iz) of current k, used after inik_k()
40+
void setuptransform();
41+
void setupIndGk();
42+
int *igl2isz_k; //[npwk_max*nks] map (ig,ik) to (is,iz)
43+
void collect_local_pw();
44+
45+
public:
46+
void real2recip(double * in, std::complex<double> * out, int ik); //in:(nplane,nx*ny) ; out(nz, ns)
47+
void real2recip(std::complex<double> * in, std::complex<double> * out, int ik); //in:(nplane,nx*ny) ; out(nz, ns)
48+
void recip2real(std::complex<double> * in, double *out, int ik); //in:(nz, ns) ; out(nplane,nx*ny)
49+
void recip2real(std::complex<double> * in, std::complex<double> * out, int ik); //in:(nz, ns) ; out(nplane,nx*ny)
50+
51+
#ifdef __MIX_PRECISION
52+
void real2recip(float * in, std::complex<float> * out, int ik); //in:(nplane,nx*ny) ; out(nz, ns)
53+
void real2recip(std::complex<float> * in, std::complex<float> * out, int ik); //in:(nplane,nx*ny) ; out(nz, ns)
54+
void recip2real(std::complex<float> * in, float *out, int ik); //in:(nz, ns) ; out(nplane,nx*ny)
55+
void recip2real(std::complex<float> * in, std::complex<float> * out, int ik); //in:(nz, ns) ; out(nplane,nx*ny)
56+
#endif
4557

4658
public:
4759
//operator
4860
ModuleBase::Vector3<double> get_GPlusK_cartesian(const int ik, const int ig) const;
4961
double get_GPlusK_cartesian_projection(const int ik, const int ig, const int axis) const;
5062
double get_SquareGPlusK_cartesian(const int ik, const int ig) const;
5163
};
64+
65+
}
5266
#endif //PlaneWave_K class
67+

‎source/module_pw/pw_operation.cpp‎

Lines changed: 47 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -1,18 +1,45 @@
11
#include "pw_basis_k.h"
2-
#include "pw_basis.h"
32
#include <assert.h>
3+
namespace ModulePW
4+
{
45

56
ModuleBase::Vector3<double> PW_Basis_K:: get_GPlusK_cartesian(const int ik, const int ig) const {
67
assert(ig>=0 && ig<this->npw && ik>=0 && ik<this->nks);
7-
ModuleBase::Vector3<double> g_temp_ = this->kvec_c[ik] + this->gcar[ig];
8+
int isz = this->ig2isz[ig];
9+
int iz = isz % this->nz;
10+
int is = isz / this->nz;
11+
int ix = this->is2ixy[is] / this->ny;
12+
int iy = this->is2ixy[is] % this->ny;
13+
if (ix >= int(this->nx/2) + 1) ix -= this->nx;
14+
if (iy >= int(this->bigny/2) + 1) iy -= this->bigny;
15+
if (iz >= int(this->nz/2) + 1) iz -= this->nz;
16+
ModuleBase::Vector3<double> f;
17+
f.x = ix;
18+
f.y = iy;
19+
f.z = iz;
20+
f = f * this->G;
21+
ModuleBase::Vector3<double> g_temp_ = this->kvec_c[ik] + f;
822
return g_temp_;
923
}
1024

1125

1226
double PW_Basis_K::get_GPlusK_cartesian_projection(const int ik, const int ig, const int axis) const
1327
{
14-
assert(ig >= 0 && ig < this->npw && ik >= 0 && ik<this->nks && axis >= 0 && axis <= 2);
15-
ModuleBase::Vector3<double> g_temp_ = this->kvec_c[ik] + this->gcar[ig];
28+
assert(ig>=0 && ig<this->npw && ik>=0 && ik<this->nks);
29+
int isz = this->ig2isz[ig];
30+
int iz = isz % this->nz;
31+
int is = isz / this->nz;
32+
int ix = this->is2ixy[is] / this->ny;
33+
int iy = this->is2ixy[is] % this->ny;
34+
if (ix >= int(this->nx/2) + 1) ix -= this->nx;
35+
if (iy >= int(this->bigny/2) + 1) iy -= this->bigny;
36+
if (iz >= int(this->nz/2) + 1) iz -= this->nz;
37+
ModuleBase::Vector3<double> f;
38+
f.x = ix;
39+
f.y = iy;
40+
f.z = iz;
41+
f = f * this->G;
42+
ModuleBase::Vector3<double> g_temp_ = this->kvec_c[ik] + f;
1643
if (axis == 0)
1744
{
1845
return g_temp_.x;
@@ -31,6 +58,21 @@ double PW_Basis_K::get_GPlusK_cartesian_projection(const int ik, const int ig, c
3158
double PW_Basis_K::get_SquareGPlusK_cartesian(const int ik, const int ig) const
3259
{
3360
assert(ig>=0 && ig<this->npw && ik>=0 && ik<this->nks);
34-
ModuleBase::Vector3<double> g_temp_ = this->kvec_c[ik] + this->gcar[ig];
61+
int isz = this->ig2isz[ig];
62+
int iz = isz % this->nz;
63+
int is = isz / this->nz;
64+
int ix = this->is2ixy[is] / this->ny;
65+
int iy = this->is2ixy[is] % this->ny;
66+
if (ix >= int(this->nx/2) + 1) ix -= this->nx;
67+
if (iy >= int(this->bigny/2) + 1) iy -= this->bigny;
68+
if (iz >= int(this->nz/2) + 1) iz -= this->nz;
69+
ModuleBase::Vector3<double> f;
70+
f.x = ix;
71+
f.y = iy;
72+
f.z = iz;
73+
f = f * this->G;
74+
ModuleBase::Vector3<double> g_temp_ = this->kvec_c[ik] + f;
3575
return (g_temp_ * g_temp_);
76+
}
77+
3678
}

0 commit comments

Comments
 (0)