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