ABACUS develop
Atomic-orbital Based Ab-initio Computation at UStc
Loading...
Searching...
No Matches
Matrix_Orbs22.hpp
Go to the documentation of this file.
1//=======================
2// AUTHOR : Peize Lin
3// DATE : 2023-02-23
4//=======================
5
6#ifndef MATRIX_ORB22_HPP
7#define MATRIX_ORB22_HPP
8
9#include "Matrix_Orbs22.h"
10#include "RI_Util.h"
11
12template<typename Tdata>
14 const size_t TA,
15 const size_t TB,
22 const Matrix_Order &matrix_order) const
23{
24 const double lat0 = *this->lat0;
25 RI::Tensor<Tdata> m;
26 const size_t sizeA1 = index_A1[TA].count_size;
27 const size_t sizeA2 = index_A2[TA].count_size;
28 const size_t sizeB1 = index_B1[TB].count_size;
29 const size_t sizeB2 = index_B2[TB].count_size;
30 switch(matrix_order)
31 {
32 case Matrix_Order::A1A2B1B2: m = RI::Tensor<Tdata>({sizeA1, sizeA2, sizeB1, sizeB2}); break;
33 case Matrix_Order::A1A2B2B1: m = RI::Tensor<Tdata>({sizeA1, sizeA2, sizeB2, sizeB1}); break;
34 case Matrix_Order::A1B1A2B2: m = RI::Tensor<Tdata>({sizeA1, sizeB1, sizeA2, sizeB2}); break;
35 case Matrix_Order::A1B1B2A2: m = RI::Tensor<Tdata>({sizeA1, sizeB1, sizeB2, sizeA2}); break;
36 case Matrix_Order::A1B2A2B1: m = RI::Tensor<Tdata>({sizeA1, sizeB2, sizeA2, sizeB1}); break;
37 case Matrix_Order::A1B2B1A2: m = RI::Tensor<Tdata>({sizeA1, sizeB2, sizeB1, sizeA2}); break;
38 case Matrix_Order::A2A1B1B2: m = RI::Tensor<Tdata>({sizeA2, sizeA1, sizeB1, sizeB2}); break;
39 case Matrix_Order::A2A1B2B1: m = RI::Tensor<Tdata>({sizeA2, sizeA1, sizeB2, sizeB1}); break;
40 case Matrix_Order::A2B1A1B2: m = RI::Tensor<Tdata>({sizeA2, sizeB1, sizeA1, sizeB2}); break;
41 case Matrix_Order::A2B1B2A1: m = RI::Tensor<Tdata>({sizeA2, sizeB1, sizeB2, sizeA1}); break;
42 case Matrix_Order::A2B2A1B1: m = RI::Tensor<Tdata>({sizeA2, sizeB2, sizeA1, sizeB1}); break;
43 case Matrix_Order::A2B2B1A1: m = RI::Tensor<Tdata>({sizeA2, sizeB2, sizeB1, sizeA1}); break;
44 case Matrix_Order::B1A1A2B2: m = RI::Tensor<Tdata>({sizeB1, sizeA1, sizeA2, sizeB2}); break;
45 case Matrix_Order::B1A1B2A2: m = RI::Tensor<Tdata>({sizeB1, sizeA1, sizeB2, sizeA2}); break;
46 case Matrix_Order::B1A2A1B2: m = RI::Tensor<Tdata>({sizeB1, sizeA2, sizeA1, sizeB2}); break;
47 case Matrix_Order::B1A2B2A1: m = RI::Tensor<Tdata>({sizeB1, sizeA2, sizeB2, sizeA1}); break;
48 case Matrix_Order::B1B2A1A2: m = RI::Tensor<Tdata>({sizeB1, sizeB2, sizeA1, sizeA2}); break;
49 case Matrix_Order::B1B2A2A1: m = RI::Tensor<Tdata>({sizeB1, sizeB2, sizeA2, sizeA1}); break;
50 case Matrix_Order::B2A1A2B1: m = RI::Tensor<Tdata>({sizeB2, sizeA1, sizeA2, sizeB1}); break;
51 case Matrix_Order::B2A1B1A2: m = RI::Tensor<Tdata>({sizeB2, sizeA1, sizeB1, sizeA2}); break;
52 case Matrix_Order::B2A2A1B1: m = RI::Tensor<Tdata>({sizeB2, sizeA2, sizeA1, sizeB1}); break;
53 case Matrix_Order::B2A2B1A1: m = RI::Tensor<Tdata>({sizeB2, sizeA2, sizeB1, sizeA1}); break;
54 case Matrix_Order::B2B1A1A2: m = RI::Tensor<Tdata>({sizeB2, sizeB1, sizeA1, sizeA2}); break;
55 case Matrix_Order::B2B1A2A1: m = RI::Tensor<Tdata>({sizeB2, sizeB1, sizeA2, sizeA1}); break;
56 default: throw std::invalid_argument(std::string(__FILE__)+" line "+std::to_string(__LINE__));
57 }
58
59 for( const auto &co3 : center2_orb22_s.at(TA).at(TB) )
60 {
61 const int LA1 = co3.first;
62 for( const auto &co4 : co3.second )
63 {
64 const size_t NA1 = co4.first;
65 for( size_t MA1=0; MA1!=2*LA1+1; ++MA1 )
66 {
67 for( const auto &co5 : co4.second )
68 {
69 const int LA2 = co5.first;
70 for( const auto &co6 : co5.second )
71 {
72 const size_t NA2 = co6.first;
73 for( size_t MA2=0; MA2!=2*LA2+1; ++MA2 )
74 {
75 for( const auto &co7 : co6.second )
76 {
77 const int LB1 = co7.first;
78 for( const auto &co8 : co7.second )
79 {
80 const size_t NB1 = co8.first;
81 for( size_t MB1=0; MB1!=2*LB1+1; ++MB1 )
82 {
83 for( const auto &co9 : co8.second )
84 {
85 const int LB2 = co9.first;
86 for( const auto &co10 : co9.second )
87 {
88 const size_t NB2 = co10.first;
89 for( size_t MB2=0; MB2!=2*LB2+1; ++MB2 )
90 {
91 const Tdata overlap = co10.second.cal_overlap( tauA*lat0, tauB*lat0, MA1, MA2, MB1, MB2 );
92 const size_t iA1 = index_A1[TA][LA1][NA1][MA1];
93 const size_t iA2 = index_A2[TA][LA2][NA2][MA2];
94 const size_t iB1 = index_B1[TB][LB1][NB1][MB1];
95 const size_t iB2 = index_B2[TB][LB2][NB2][MB2];
96 switch(matrix_order)
97 {
98 case Matrix_Order::A1A2B1B2: m(iA1,iA2,iB1,iB2) = overlap; break;
99 case Matrix_Order::A1A2B2B1: m(iA1,iA2,iB2,iB1) = overlap; break;
100 case Matrix_Order::A1B1A2B2: m(iA1,iB1,iA2,iB2) = overlap; break;
101 case Matrix_Order::A1B1B2A2: m(iA1,iB1,iB2,iA2) = overlap; break;
102 case Matrix_Order::A1B2A2B1: m(iA1,iB2,iA2,iB1) = overlap; break;
103 case Matrix_Order::A1B2B1A2: m(iA1,iB2,iB1,iA2) = overlap; break;
104 case Matrix_Order::A2A1B1B2: m(iA2,iA1,iB1,iB2) = overlap; break;
105 case Matrix_Order::A2A1B2B1: m(iA2,iA1,iB2,iB1) = overlap; break;
106 case Matrix_Order::A2B1A1B2: m(iA2,iB1,iA1,iB2) = overlap; break;
107 case Matrix_Order::A2B1B2A1: m(iA2,iB1,iB2,iA1) = overlap; break;
108 case Matrix_Order::A2B2A1B1: m(iA2,iB2,iA1,iB1) = overlap; break;
109 case Matrix_Order::A2B2B1A1: m(iA2,iB2,iB1,iA1) = overlap; break;
110 case Matrix_Order::B1A1A2B2: m(iB1,iA1,iA2,iB2) = overlap; break;
111 case Matrix_Order::B1A1B2A2: m(iB1,iA1,iB2,iA2) = overlap; break;
112 case Matrix_Order::B1A2A1B2: m(iB1,iA2,iA1,iB2) = overlap; break;
113 case Matrix_Order::B1A2B2A1: m(iB1,iA2,iB2,iA1) = overlap; break;
114 case Matrix_Order::B1B2A1A2: m(iB1,iB2,iA1,iA2) = overlap; break;
115 case Matrix_Order::B1B2A2A1: m(iB1,iB2,iA2,iA1) = overlap; break;
116 case Matrix_Order::B2A1A2B1: m(iB2,iA1,iA2,iB1) = overlap; break;
117 case Matrix_Order::B2A1B1A2: m(iB2,iA1,iB1,iA2) = overlap; break;
118 case Matrix_Order::B2A2A1B1: m(iB2,iA2,iA1,iB1) = overlap; break;
119 case Matrix_Order::B2A2B1A1: m(iB2,iA2,iB1,iA1) = overlap; break;
120 case Matrix_Order::B2B1A1A2: m(iB2,iB1,iA1,iA2) = overlap; break;
121 case Matrix_Order::B2B1A2A1: m(iB2,iB1,iA2,iA1) = overlap; break;
122 default: throw std::invalid_argument(std::string(__FILE__)+" line "+std::to_string(__LINE__));
123 }
124 }
125 }
126 }
127 }
128 }
129 }
130 }
131 }
132 }
133 }
134 }
135 }
136 return m;
137}
138
139template<typename Tdata>
140std::array<RI::Tensor<Tdata>,3> Matrix_Orbs22::cal_grad_overlap_matrix(
141 const size_t TA,
142 const size_t TB,
143 const ModuleBase::Vector3<double> &tauA,
144 const ModuleBase::Vector3<double> &tauB,
149 const Matrix_Order &matrix_order) const
150{
151 std::array<RI::Tensor<Tdata>,3> m;
152 const size_t sizeA1 = index_A1[TA].count_size;
153 const size_t sizeA2 = index_A2[TA].count_size;
154 const size_t sizeB1 = index_B1[TB].count_size;
155 const size_t sizeB2 = index_B2[TB].count_size;
156 for(int i=0; i<m.size(); ++i)
157 {
158 switch(matrix_order)
159 {
160 case Matrix_Order::A1A2B1B2: m[i] = RI::Tensor<Tdata>({sizeA1, sizeA2, sizeB1, sizeB2}); break;
161 case Matrix_Order::A1A2B2B1: m[i] = RI::Tensor<Tdata>({sizeA1, sizeA2, sizeB2, sizeB1}); break;
162 case Matrix_Order::A1B1A2B2: m[i] = RI::Tensor<Tdata>({sizeA1, sizeB1, sizeA2, sizeB2}); break;
163 case Matrix_Order::A1B1B2A2: m[i] = RI::Tensor<Tdata>({sizeA1, sizeB1, sizeB2, sizeA2}); break;
164 case Matrix_Order::A1B2A2B1: m[i] = RI::Tensor<Tdata>({sizeA1, sizeB2, sizeA2, sizeB1}); break;
165 case Matrix_Order::A1B2B1A2: m[i] = RI::Tensor<Tdata>({sizeA1, sizeB2, sizeB1, sizeA2}); break;
166 case Matrix_Order::A2A1B1B2: m[i] = RI::Tensor<Tdata>({sizeA2, sizeA1, sizeB1, sizeB2}); break;
167 case Matrix_Order::A2A1B2B1: m[i] = RI::Tensor<Tdata>({sizeA2, sizeA1, sizeB2, sizeB1}); break;
168 case Matrix_Order::A2B1A1B2: m[i] = RI::Tensor<Tdata>({sizeA2, sizeB1, sizeA1, sizeB2}); break;
169 case Matrix_Order::A2B1B2A1: m[i] = RI::Tensor<Tdata>({sizeA2, sizeB1, sizeB2, sizeA1}); break;
170 case Matrix_Order::A2B2A1B1: m[i] = RI::Tensor<Tdata>({sizeA2, sizeB2, sizeA1, sizeB1}); break;
171 case Matrix_Order::A2B2B1A1: m[i] = RI::Tensor<Tdata>({sizeA2, sizeB2, sizeB1, sizeA1}); break;
172 case Matrix_Order::B1A1A2B2: m[i] = RI::Tensor<Tdata>({sizeB1, sizeA1, sizeA2, sizeB2}); break;
173 case Matrix_Order::B1A1B2A2: m[i] = RI::Tensor<Tdata>({sizeB1, sizeA1, sizeB2, sizeA2}); break;
174 case Matrix_Order::B1A2A1B2: m[i] = RI::Tensor<Tdata>({sizeB1, sizeA2, sizeA1, sizeB2}); break;
175 case Matrix_Order::B1A2B2A1: m[i] = RI::Tensor<Tdata>({sizeB1, sizeA2, sizeB2, sizeA1}); break;
176 case Matrix_Order::B1B2A1A2: m[i] = RI::Tensor<Tdata>({sizeB1, sizeB2, sizeA1, sizeA2}); break;
177 case Matrix_Order::B1B2A2A1: m[i] = RI::Tensor<Tdata>({sizeB1, sizeB2, sizeA2, sizeA1}); break;
178 case Matrix_Order::B2A1A2B1: m[i] = RI::Tensor<Tdata>({sizeB2, sizeA1, sizeA2, sizeB1}); break;
179 case Matrix_Order::B2A1B1A2: m[i] = RI::Tensor<Tdata>({sizeB2, sizeA1, sizeB1, sizeA2}); break;
180 case Matrix_Order::B2A2A1B1: m[i] = RI::Tensor<Tdata>({sizeB2, sizeA2, sizeA1, sizeB1}); break;
181 case Matrix_Order::B2A2B1A1: m[i] = RI::Tensor<Tdata>({sizeB2, sizeA2, sizeB1, sizeA1}); break;
182 case Matrix_Order::B2B1A1A2: m[i] = RI::Tensor<Tdata>({sizeB2, sizeB1, sizeA1, sizeA2}); break;
183 case Matrix_Order::B2B1A2A1: m[i] = RI::Tensor<Tdata>({sizeB2, sizeB1, sizeA2, sizeA1}); break;
184 default: throw std::invalid_argument(std::string(__FILE__)+" line "+std::to_string(__LINE__));
185 }
186 }
187 const double lat0 = *this->lat0;
188 for( const auto &co3 : center2_orb22_s.at(TA).at(TB) )
189 {
190 const int LA1 = co3.first;
191 for( const auto &co4 : co3.second )
192 {
193 const size_t NA1 = co4.first;
194 for( size_t MA1=0; MA1!=2*LA1+1; ++MA1 )
195 {
196 for( const auto &co5 : co4.second )
197 {
198 const int LA2 = co5.first;
199 for( const auto &co6 : co5.second )
200 {
201 const size_t NA2 = co6.first;
202 for( size_t MA2=0; MA2!=2*LA2+1; ++MA2 )
203 {
204 for( const auto &co7 : co6.second )
205 {
206 const int LB1 = co7.first;
207 for( const auto &co8 : co7.second )
208 {
209 const size_t NB1 = co8.first;
210 for( size_t MB1=0; MB1!=2*LB1+1; ++MB1 )
211 {
212 for( const auto &co9 : co8.second )
213 {
214 const int LB2 = co9.first;
215 for( const auto &co10 : co9.second )
216 {
217 const size_t NB2 = co10.first;
218 for( size_t MB2=0; MB2!=2*LB2+1; ++MB2 )
219 {
220 const std::array<double,3> grad_overlap = RI_Util::Vector3_to_array3(co10.second.cal_grad_overlap( tauA*lat0, tauB*lat0, MA1, MA2, MB1, MB2 ));
221 const size_t iA1 = index_A1[TA][LA1][NA1][MA1];
222 const size_t iA2 = index_A2[TA][LA2][NA2][MA2];
223 const size_t iB1 = index_B1[TB][LB1][NB1][MB1];
224 const size_t iB2 = index_B2[TB][LB2][NB2][MB2];
225 for(size_t i=0; i<m.size(); ++i)
226 {
227 switch(matrix_order)
228 {
229 case Matrix_Order::A1A2B1B2: m[i](iA1,iA2,iB1,iB2) = grad_overlap; break;
230 case Matrix_Order::A1A2B2B1: m[i](iA1,iA2,iB2,iB1) = grad_overlap; break;
231 case Matrix_Order::A1B1A2B2: m[i](iA1,iB1,iA2,iB2) = grad_overlap; break;
232 case Matrix_Order::A1B1B2A2: m[i](iA1,iB1,iB2,iA2) = grad_overlap; break;
233 case Matrix_Order::A1B2A2B1: m[i](iA1,iB2,iA2,iB1) = grad_overlap; break;
234 case Matrix_Order::A1B2B1A2: m[i](iA1,iB2,iB1,iA2) = grad_overlap; break;
235 case Matrix_Order::A2A1B1B2: m[i](iA2,iA1,iB1,iB2) = grad_overlap; break;
236 case Matrix_Order::A2A1B2B1: m[i](iA2,iA1,iB2,iB1) = grad_overlap; break;
237 case Matrix_Order::A2B1A1B2: m[i](iA2,iB1,iA1,iB2) = grad_overlap; break;
238 case Matrix_Order::A2B1B2A1: m[i](iA2,iB1,iB2,iA1) = grad_overlap; break;
239 case Matrix_Order::A2B2A1B1: m[i](iA2,iB2,iA1,iB1) = grad_overlap; break;
240 case Matrix_Order::A2B2B1A1: m[i](iA2,iB2,iB1,iA1) = grad_overlap; break;
241 case Matrix_Order::B1A1A2B2: m[i](iB1,iA1,iA2,iB2) = grad_overlap; break;
242 case Matrix_Order::B1A1B2A2: m[i](iB1,iA1,iB2,iA2) = grad_overlap; break;
243 case Matrix_Order::B1A2A1B2: m[i](iB1,iA2,iA1,iB2) = grad_overlap; break;
244 case Matrix_Order::B1A2B2A1: m[i](iB1,iA2,iB2,iA1) = grad_overlap; break;
245 case Matrix_Order::B1B2A1A2: m[i](iB1,iB2,iA1,iA2) = grad_overlap; break;
246 case Matrix_Order::B1B2A2A1: m[i](iB1,iB2,iA2,iA1) = grad_overlap; break;
247 case Matrix_Order::B2A1A2B1: m[i](iB2,iA1,iA2,iB1) = grad_overlap; break;
248 case Matrix_Order::B2A1B1A2: m[i](iB2,iA1,iB1,iA2) = grad_overlap; break;
249 case Matrix_Order::B2A2A1B1: m[i](iB2,iA2,iA1,iB1) = grad_overlap; break;
250 case Matrix_Order::B2A2B1A1: m[i](iB2,iA2,iB1,iA1) = grad_overlap; break;
251 case Matrix_Order::B2B1A1A2: m[i](iB2,iB1,iA1,iA2) = grad_overlap; break;
252 case Matrix_Order::B2B1A2A1: m[i](iB2,iB1,iA2,iA1) = grad_overlap; break;
253 default: throw std::invalid_argument(std::string(__FILE__)+" line "+std::to_string(__LINE__));
254 }
255 }
256 }
257 }
258 }
259 }
260 }
261 }
262 }
263 }
264 }
265 }
266 }
267 }
268 return m;
269}
270
271// 4-parameter interface, for opt_orb
272template<typename Tdata>
273std::map < size_t, std::map<size_t, std::map<size_t, std::map<size_t, RI::Tensor<Tdata>>>>> Matrix_Orbs22::cal_overlap_matrix_all(
274 const UnitCell &ucell,
278 const ModuleBase::Element_Basis_Index::IndexLNM &index_B2 ) const
279{
280 std::map<size_t,std::map<size_t,std::map<size_t,std::map<size_t, RI::Tensor<Tdata>>>>> matrixes;
281
282 for( const auto &co1 : center2_orb22_s )
283 {
284 const size_t TA = co1.first;
285 for( size_t IA=0; IA!=ucell.atoms[TA].na; ++IA )
286 {
287 const ModuleBase::Vector3<double> &tauA( ucell.atoms[TA].tau[IA] );
288
289 for( const auto &co2 : co1.second )
290 {
291 const size_t TB = co2.first;
292 for( size_t IB=0; IB!=ucell.atoms[TB].na; ++IB )
293 {
294 const ModuleBase::Vector3<double> &tauB( ucell.atoms[TB].tau[IB] );
295
296 matrixes[TA][IA][TB][IB] = cal_overlap_matrix<Tdata>(
297 TA,
298 TB,
299 ucell.atoms[TA].tau[IA],
300 ucell.atoms[TB].tau[IB],
301 index_A1,
302 index_A2,
303 index_B1,
304 index_B2,
306 }
307 }
308 }
309 }
310 return matrixes;
311}
312#endif
const std::complex< double > i
Definition cal_pLpR.cpp:46
int na
Definition atom_spec.h:27
std::vector< ModuleBase::Vector3< double > > tau
Definition atom_spec.h:35
std::array< RI::Tensor< Tdata >, 3 > cal_grad_overlap_matrix(const size_t TA, const size_t TB, const ModuleBase::Vector3< double > &tauA, const ModuleBase::Vector3< double > &tauB, const ModuleBase::Element_Basis_Index::IndexLNM &index_A1, const ModuleBase::Element_Basis_Index::IndexLNM &index_A2, const ModuleBase::Element_Basis_Index::IndexLNM &index_B1, const ModuleBase::Element_Basis_Index::IndexLNM &index_B2, const Matrix_Order &matrix_order) const
Definition Matrix_Orbs22.hpp:140
std::map< size_t, std::map< size_t, std::map< size_t, std::map< size_t, RI::Tensor< Tdata > > > > > cal_overlap_matrix_all(const UnitCell &ucell, const ModuleBase::Element_Basis_Index::IndexLNM &index_A1, const ModuleBase::Element_Basis_Index::IndexLNM &index_A2, const ModuleBase::Element_Basis_Index::IndexLNM &index_B1, const ModuleBase::Element_Basis_Index::IndexLNM &index_B2) const
Definition Matrix_Orbs22.hpp:273
std::map< size_t, std::map< size_t, std::map< int, std::map< size_t, std::map< int, std::map< size_t, std::map< int, std::map< size_t, std::map< int, std::map< size_t, Center2_Orb::Orb22 > > > > > > > > > > center2_orb22_s
Definition Matrix_Orbs22.h:112
Matrix_Order
Definition Matrix_Orbs22.h:37
double * lat0
Definition Matrix_Orbs22.h:99
RI::Tensor< Tdata > cal_overlap_matrix(const size_t TA, const size_t TB, const ModuleBase::Vector3< double > &tauA, const ModuleBase::Vector3< double > &tauB, const ModuleBase::Element_Basis_Index::IndexLNM &index_A1, const ModuleBase::Element_Basis_Index::IndexLNM &index_A2, const ModuleBase::Element_Basis_Index::IndexLNM &index_B1, const ModuleBase::Element_Basis_Index::IndexLNM &index_B2, const Matrix_Order &matrix_order) const
Definition Matrix_Orbs22.hpp:13
3 elements vector
Definition vector3.h:24
Definition unitcell.h:15
Atom * atoms
Definition unitcell.h:45
std::vector< Index_T > IndexLNM
Definition element_basis_index.h:42
std::array< Tcell, 3 > Vector3_to_array3(const ModuleBase::Vector3< Tcell > &v)
Definition RI_Util.h:32