1+ // wenfei 2022-1-5
2+ // This file contains constructor and destructor of the class LCAO_deepks,
3+ // as well as subroutines for initializing and releasing relevant data structures
4+
5+ // Other than the constructor and the destructor, it contains 3 types of subroutines:
6+ // 1. subroutines that are related to calculating descriptors:
7+ // - init : allocates some arrays
8+ // - init_index : records the index (inl)
9+ // - allocate_nlm : allocates data structures (nlm_save) which is used to store <chi|alpha>
10+ // 2. subroutines that are related to calculating force label:
11+ // - init_gdmx : allocates gdmx; it is a private subroutine
12+ // - del_gdmx : releases gdmx
13+ // 3. subroutines that are related to V_delta:
14+ // - allocate_V_delta : allocates H_V_delta; if calculating force, it also calls
15+ // init_gdmx, as well as allocating F_delta
16+ // - allocate_V_deltaR : allcoates H_V_deltaR, for multi-k calculations
17+
18+ #ifdef __DEEPKS
19+
20+ #include " LCAO_deepks.h"
21+
22+ namespace GlobalC
23+ {
24+ LCAO_Deepks ld;
25+ }
26+
27+ // Constructor of the class
28+ LCAO_Deepks::LCAO_Deepks ()
29+ {
30+ alpha_index = new ModuleBase::IntArray[1 ];
31+ inl_index = new ModuleBase::IntArray[1 ];
32+ inl_l = new int [1 ];
33+ H_V_delta = new double [1 ];
34+ }
35+
36+ // Desctructor of the class
37+ LCAO_Deepks::~LCAO_Deepks ()
38+ {
39+ delete[] alpha_index;
40+ delete[] inl_index;
41+ delete[] inl_l;
42+ delete[] H_V_delta;
43+
44+ // =======1. "out_descriptor" part==========
45+ // delete pdm**
46+ for (int inl = 0 ;inl < this ->inlmax ;inl++)
47+ {
48+ delete[] pdm[inl];
49+ }
50+ delete[] pdm;
51+ // =======2. "deepks_scf" part==========
52+ if (GlobalV::deepks_scf)
53+ {
54+ // delete gedm**
55+ for (int inl = 0 ;inl < this ->inlmax ;inl++)
56+ {
57+ delete[] gedm[inl];
58+ }
59+ delete[] gedm;
60+ }
61+ }
62+
63+ void LCAO_Deepks::init (
64+ const LCAO_Orbitals &orb,
65+ const int nat,
66+ const int ntype,
67+ std::vector<int > na)
68+ {
69+ ModuleBase::TITLE (" LCAO_Deepks" " init"
70+
71+ GlobalV::ofs_running << " Initialize the descriptor index for DeePKS (lcao line)"
72+
73+ const int lm = orb.get_lmax_d ();
74+ const int nm = orb.get_nchimax_d ();
75+ const int tot_inl_per_atom = orb.Alpha [0 ].getTotal_nchi ();
76+
77+ assert (lm >= 0 );
78+ assert (nm >= 0 );
79+ assert (tot_inl_per_atom >= 0 );
80+
81+ const int tot_inl = tot_inl_per_atom * nat;
82+
83+ this ->lmaxd = lm;
84+ this ->nmaxd = nm;
85+ this ->inlmax = tot_inl;
86+ GlobalV::ofs_running << " lmax of descriptor = " this ->lmaxd << std::endl;
87+ GlobalV::ofs_running << " nmax of descriptor= "
88+ GlobalV::ofs_running << " total basis (all atoms) for descriptor= "
89+
90+ // init pdm**
91+ const int pdm_size = (this ->lmaxd * 2 + 1 ) * (this ->lmaxd * 2 + 1 );
92+ this ->pdm = new double * [this ->inlmax ];
93+ for (int inl = 0 ;inl < this ->inlmax ;inl++)
94+ {
95+ this ->pdm [inl] = new double [pdm_size];
96+ ModuleBase::GlobalFunc::ZEROS (this ->pdm [inl], pdm_size);
97+ }
98+
99+ // cal n(descriptor) per atom , related to Lmax, nchi(L) and m. (not total_nchi!)
100+ this ->des_per_atom =0 ; // mohan add 2021-04-21
101+ for (int l = 0 ; l <= this ->lmaxd ; l++)
102+ {
103+ this ->des_per_atom += orb.Alpha [0 ].getNchi (l) * (2 * l + 1 );
104+ }
105+
106+ this ->n_descriptor = nat * this ->des_per_atom ;
107+
108+ this ->init_index (ntype, nat, na, tot_inl, orb);
109+ this ->allocate_nlm (nat);
110+
111+ return ;
112+ }
113+
114+ void LCAO_Deepks::init_index (const int ntype, const int nat, std::vector<int > na, const int Total_nchi, const LCAO_Orbitals &orb)
115+ {
116+ delete[] this ->alpha_index ;
117+ this ->alpha_index = new ModuleBase::IntArray[ntype];
118+ delete[] this ->inl_index ;
119+ this ->inl_index = new ModuleBase::IntArray[ntype];
120+ delete[] this ->inl_l ;
121+ this ->inl_l = new int [this ->inlmax ];
122+ ModuleBase::GlobalFunc::ZEROS (this ->inl_l , this ->inlmax );
123+
124+ int inl = 0 ;
125+ int alpha = 0 ;
126+ for (int it = 0 ; it < ntype; it++)
127+ {
128+ this ->alpha_index [it].create (
129+ na[it],
130+ this ->lmaxd + 1 , // l starts from 0
131+ this ->nmaxd ,
132+ 2 * this ->lmaxd + 1 ); // m ==> 2*l+1
133+
134+ this ->inl_index [it].create (
135+ na[it],
136+ this ->lmaxd + 1 ,
137+ this ->nmaxd );
138+
139+ GlobalV::ofs_running << " Type " 1
140+ << " number_of_atoms "
141+
142+ for (int ia = 0 ; ia < na[it]; ia++)
143+ {
144+ // alpha
145+ for (int l = 0 ; l < this ->lmaxd + 1 ; l++)
146+ {
147+ for (int n = 0 ; n < orb.Alpha [0 ].getNchi (l); n++)
148+ {
149+ for (int m = 0 ; m < 2 * l + 1 ; m++)
150+ {
151+ this ->alpha_index [it](ia, l, n, m) = alpha;
152+ alpha++;
153+ }
154+ this ->inl_index [it](ia, l, n) = inl;
155+ this ->inl_l [inl] = l;
156+ inl++;
157+ }
158+ }
159+ }// end ia
160+ }// end it
161+ assert (this ->n_descriptor == alpha);
162+ assert (Total_nchi == inl);
163+ GlobalV::ofs_running << " descriptors_per_atom " this ->des_per_atom << std::endl;
164+ GlobalV::ofs_running << " total_descriptors " this ->n_descriptor << std::endl;
165+ return ;
166+ }
167+
168+ void LCAO_Deepks::allocate_nlm (const int nat)
169+ {
170+ if (GlobalV::GAMMA_ONLY_LOCAL)
171+ {
172+ this ->nlm_save .resize (nat);
173+ }
174+ else
175+ {
176+ this ->nlm_save_k .resize (nat);
177+ }
178+ }
179+
180+ void LCAO_Deepks::init_gdmx (const int nat)
181+ {
182+ this ->gdmx = new double ** [nat];
183+ this ->gdmy = new double ** [nat];
184+ this ->gdmz = new double ** [nat];
185+ for (int iat = 0 ;iat < nat;iat++)
186+ {
187+ this ->gdmx [iat] = new double * [inlmax];
188+ this ->gdmy [iat] = new double * [inlmax];
189+ this ->gdmz [iat] = new double * [inlmax];
190+ for (int inl = 0 ;inl < inlmax;inl++)
191+ {
192+ this ->gdmx [iat][inl] = new double [(2 * lmaxd + 1 ) * (2 * lmaxd + 1 )];
193+ this ->gdmy [iat][inl] = new double [(2 * lmaxd + 1 ) * (2 * lmaxd + 1 )];
194+ this ->gdmz [iat][inl] = new double [(2 * lmaxd + 1 ) * (2 * lmaxd + 1 )];
195+ ModuleBase::GlobalFunc::ZEROS (gdmx[iat][inl], (2 * lmaxd + 1 ) * (2 * lmaxd + 1 ));
196+ ModuleBase::GlobalFunc::ZEROS (gdmy[iat][inl], (2 * lmaxd + 1 ) * (2 * lmaxd + 1 ));
197+ ModuleBase::GlobalFunc::ZEROS (gdmz[iat][inl], (2 * lmaxd + 1 ) * (2 * lmaxd + 1 ));
198+ }
199+ }
200+ return ;
201+ }
202+
203+ void LCAO_Deepks::del_gdmx (const int nat)
204+ {
205+ for (int iat = 0 ;iat < nat;iat++)
206+ {
207+ for (int inl = 0 ;inl < inlmax;inl++)
208+ {
209+ delete[] this ->gdmx [iat][inl];
210+ delete[] this ->gdmy [iat][inl];
211+ delete[] this ->gdmz [iat][inl];
212+ }
213+ delete[] this ->gdmx [iat];
214+ delete[] this ->gdmy [iat];
215+ delete[] this ->gdmz [iat];
216+ }
217+ delete[] this ->gdmx ;
218+ delete[] this ->gdmy ;
219+ delete[] this ->gdmz ;
220+ return ;
221+ }
222+
223+ void LCAO_Deepks::allocate_V_delta (const int nat, const int nloc, const int nks)
224+ {
225+ ModuleBase::TITLE (" LCAO_Deepks" " allocate_V_delta"
226+
227+ // initialize the H matrix H_V_delta
228+ if (GlobalV::GAMMA_ONLY_LOCAL)
229+ {
230+ delete[] this ->H_V_delta ;
231+ this ->H_V_delta = new double [nloc];
232+ ModuleBase::GlobalFunc::ZEROS (this ->H_V_delta , nloc);
233+ }
234+ else
235+ {
236+ H_V_delta_k = new std::complex <double >* [nks];
237+ for (int ik=0 ;ik<nks;ik++)
238+ {
239+ this ->H_V_delta_k [ik] = new std::complex <double >[nloc];
240+ ModuleBase::GlobalFunc::ZEROS (this ->H_V_delta_k [ik], nloc);
241+ }
242+ }
243+
244+ // init gedm**
245+ const int pdm_size = (this ->lmaxd * 2 + 1 ) * (this ->lmaxd * 2 + 1 );
246+ this ->gedm = new double * [this ->inlmax ];
247+ for (int inl = 0 ;inl < this ->inlmax ;inl++)
248+ {
249+ this ->gedm [inl] = new double [pdm_size];
250+ ModuleBase::GlobalFunc::ZEROS (this ->gedm [inl], pdm_size);
251+ }
252+ if (GlobalV::FORCE)
253+ {
254+ // init F_delta
255+ F_delta.create (nat, 3 );
256+ this ->init_gdmx (nat);
257+ }
258+
259+ return ;
260+ }
261+
262+ void LCAO_Deepks::allocate_V_deltaR (const int nnr)
263+ {
264+ delete[] H_V_deltaR;
265+ H_V_deltaR = new double [nnr];
266+ ModuleBase::GlobalFunc::ZEROS (H_V_deltaR, nnr);
267+ }
268+
269+ #endif
0 commit comments