ABACUS develop
Atomic-orbital Based Ab-initio Computation at UStc
Loading...
Searching...
No Matches
Matrix_Orbs21.hpp
Go to the documentation of this file.
1//=======================
2// AUTHOR : Peize Lin
3// DATE : 2022-08-17
4//=======================
5
6#ifndef MATRIX_ORB21_HPP
7#define MATRIX_ORB21_HPP
8
9#include "Matrix_Orbs21.h"
10#include "RI_Util.h"
11
12template<typename Tdata>
13RI::Tensor<Tdata> Matrix_Orbs21::cal_overlap_matrix(
14 const size_t TA,
15 const size_t TB,
16 const ModuleBase::Vector3<double> &tauA,
17 const ModuleBase::Vector3<double> &tauB,
18 const ModuleBase::Element_Basis_Index::IndexLNM &index_A1,
19 const ModuleBase::Element_Basis_Index::IndexLNM &index_A2,
20 const ModuleBase::Element_Basis_Index::IndexLNM &index_B,
21 const Matrix_Order &matrix_order) const
22{
23 RI::Tensor<Tdata> m;
24 const double lat0 = *this->lat0;
25 const size_t sizeA1 = index_A1[TA].count_size;
26 const size_t sizeA2 = index_A2[TA].count_size;
27 const size_t sizeB = index_B[TB].count_size;
28 switch(matrix_order)
29 {
30 case Matrix_Order::A1A2B: m = RI::Tensor<Tdata>({sizeA1, sizeA2, sizeB}); break;
31 case Matrix_Order::A1BA2: m = RI::Tensor<Tdata>({sizeA1, sizeB, sizeA2}); break;
32 case Matrix_Order::BA1A2: m = RI::Tensor<Tdata>({sizeB, sizeA1, sizeA2}); break;
33 case Matrix_Order::BA2A1: m = RI::Tensor<Tdata>({sizeB, sizeA2, sizeA1}); break;
34 case Matrix_Order::A2A1B: m = RI::Tensor<Tdata>({sizeA2, sizeA1, sizeB}); break;
35 case Matrix_Order::A2BA1: m = RI::Tensor<Tdata>({sizeA2, sizeB, sizeA1}); break;
36 default: throw std::invalid_argument(std::string(__FILE__)+" line "+std::to_string(__LINE__));
37 }
38
39 for( const auto &co3 : center2_orb21_s.at(TA).at(TB) )
40 {
41 const int LA1 = co3.first;
42 for( const auto &co4 : co3.second )
43 {
44 const size_t NA1 = co4.first;
45 for( size_t MA1=0; MA1!=2*LA1+1; ++MA1 )
46 {
47 for( const auto &co5 : co4.second )
48 {
49 const int LA2 = co5.first;
50 for( const auto &co6 : co5.second )
51 {
52 const size_t NA2 = co6.first;
53 for( size_t MA2=0; MA2!=2*LA2+1; ++MA2 )
54 {
55 for( const auto &co7 : co6.second )
56 {
57 const int LB = co7.first;
58 for( const auto &co8 : co7.second )
59 {
60 const size_t NB = co8.first;
61 for( size_t MB=0; MB!=2*LB+1; ++MB )
62 {
63 const Tdata overlap = co8.second.cal_overlap( tauA*lat0, tauB*lat0, MA1, MA2, MB );
64 const size_t iA1 = index_A1[TA][LA1][NA1][MA1];
65 const size_t iA2 = index_A2[TA][LA2][NA2][MA2];
66 const size_t iB = index_B[TB][LB][NB][MB];
67 switch(matrix_order)
68 {
69 case Matrix_Order::A1A2B: m(iA1,iA2,iB) = overlap; break;
70 case Matrix_Order::A1BA2: m(iA1,iB,iA2) = overlap; break;
71 case Matrix_Order::A2A1B: m(iA2,iA1,iB) = overlap; break;
72 case Matrix_Order::A2BA1: m(iA2,iB,iA1) = overlap; break;
73 case Matrix_Order::BA1A2: m(iB,iA1,iA2) = overlap; break;
74 case Matrix_Order::BA2A1: m(iB,iA2,iA1) = overlap; break;
75 default: throw std::invalid_argument(std::string(__FILE__)+" line "+std::to_string(__LINE__));
76 }
77 }
78 }
79 }
80 }
81 }
82 }
83 }
84 }
85 }
86 return m;
87}
88
89template<typename Tdata>
90std::array<RI::Tensor<Tdata>,3> Matrix_Orbs21::cal_grad_overlap_matrix(
91 const size_t TA,
92 const size_t TB,
93 const ModuleBase::Vector3<double> &tauA,
94 const ModuleBase::Vector3<double> &tauB,
95 const ModuleBase::Element_Basis_Index::IndexLNM &index_A1,
96 const ModuleBase::Element_Basis_Index::IndexLNM &index_A2,
97 const ModuleBase::Element_Basis_Index::IndexLNM &index_B,
98 const Matrix_Order &matrix_order) const
99{
100 std::array<RI::Tensor<Tdata>,3> m;
101 const double lat0 = *this->lat0;
102 const size_t sizeA1 = index_A1[TA].count_size;
103 const size_t sizeA2 = index_A2[TA].count_size;
104 const size_t sizeB = index_B[TB].count_size;
105 for(int i=0; i<m.size(); ++i)
106 {
107 switch(matrix_order)
108 {
109 case Matrix_Order::A1A2B: m[i] = RI::Tensor<Tdata>({sizeA1, sizeA2, sizeB}); break;
110 case Matrix_Order::A1BA2: m[i] = RI::Tensor<Tdata>({sizeA1, sizeB, sizeA2}); break;
111 case Matrix_Order::BA1A2: m[i] = RI::Tensor<Tdata>({sizeB, sizeA1, sizeA2}); break;
112 case Matrix_Order::BA2A1: m[i] = RI::Tensor<Tdata>({sizeB, sizeA2, sizeA1}); break;
113 case Matrix_Order::A2A1B: m[i] = RI::Tensor<Tdata>({sizeA2, sizeA1, sizeB}); break;
114 case Matrix_Order::A2BA1: m[i] = RI::Tensor<Tdata>({sizeA2, sizeB, sizeA1}); break;
115 default: throw std::invalid_argument(std::string(__FILE__)+" line "+std::to_string(__LINE__));
116 }
117 }
118
119 for( const auto &co3 : center2_orb21_s.at(TA).at(TB) )
120 {
121 const int LA1 = co3.first;
122 for( const auto &co4 : co3.second )
123 {
124 const size_t NA1 = co4.first;
125 for( size_t MA1=0; MA1!=2*LA1+1; ++MA1 )
126 {
127 for( const auto &co5 : co4.second )
128 {
129 const int LA2 = co5.first;
130 for( const auto &co6 : co5.second )
131 {
132 const size_t NA2 = co6.first;
133 for( size_t MA2=0; MA2!=2*LA2+1; ++MA2 )
134 {
135 for( const auto &co7 : co6.second )
136 {
137 const int LB = co7.first;
138 for( const auto &co8 : co7.second )
139 {
140 const size_t NB = co8.first;
141 for( size_t MB=0; MB!=2*LB+1; ++MB )
142 {
143 const std::array<double,3> grad_overlap = RI_Util::Vector3_to_array3(co8.second.cal_grad_overlap( tauA*lat0, tauB*lat0, MA1, MA2, MB ));
144 const size_t iA1 = index_A1[TA][LA1][NA1][MA1];
145 const size_t iA2 = index_A2[TA][LA2][NA2][MA2];
146 const size_t iB = index_B[TB][LB][NB][MB];
147 for(size_t i=0; i<m.size(); ++i)
148 {
149 switch(matrix_order)
150 {
151 case Matrix_Order::A1A2B: m[i](iA1,iA2,iB) = grad_overlap[i]; break;
152 case Matrix_Order::A1BA2: m[i](iA1,iB,iA2) = grad_overlap[i]; break;
153 case Matrix_Order::A2A1B: m[i](iA2,iA1,iB) = grad_overlap[i]; break;
154 case Matrix_Order::A2BA1: m[i](iA2,iB,iA1) = grad_overlap[i]; break;
155 case Matrix_Order::BA1A2: m[i](iB,iA1,iA2) = grad_overlap[i]; break;
156 case Matrix_Order::BA2A1: m[i](iB,iA2,iA1) = grad_overlap[i]; break;
157 default: throw std::invalid_argument(std::string(__FILE__)+" line "+std::to_string(__LINE__));
158 }
159 }
160 }
161 }
162 }
163 }
164 }
165 }
166 }
167 }
168 }
169 return m;
170}
171
172template <typename Tdata>
173std::map<size_t, std::map<size_t, std::map<size_t, std::map<size_t, std::vector<RI::Tensor<Tdata>>>>>> Matrix_Orbs21::cal_overlap_matrix_all(
174 const UnitCell &ucell,
175 const ModuleBase::Element_Basis_Index::IndexLNM &index_A1,
176 const ModuleBase::Element_Basis_Index::IndexLNM &index_A2,
177 const ModuleBase::Element_Basis_Index::IndexLNM &index_B) const
178{
179 ModuleBase::TITLE("Matrix_Orbs21","cal_overlap_matrix");
180
181 std::map<size_t,std::map<size_t,std::map<size_t,std::map<size_t,std::vector<RI::Tensor<Tdata>>>>>> matrixes;
182
183 for( const auto &co1 : center2_orb21_s )
184 {
185 const size_t TA = co1.first;
186 for( size_t IA=0; IA!=ucell.atoms[TA].na; ++IA )
187 {
188 const ModuleBase::Vector3<double> &tauA( ucell.atoms[TA].tau[IA] );
189
190 for( const auto &co2 : co1.second )
191 {
192 const size_t TB = co2.first;
193 for( size_t IB=0; IB!=ucell.atoms[TB].na; ++IB )
194 {
195 const ModuleBase::Vector3<double> &tauB( ucell.atoms[TB].tau[IB] );
196
197 const RI::Tensor<Tdata> &&m = cal_overlap_matrix<Tdata>( TA, TB, tauA, tauB, index_A1, index_A2, index_B, Matrix_Order::A2BA1 );
198 matrixes[TA][IA][TB][IB].resize(2);
199 matrixes[TA][IA][TB][IB][0] = std::move(m);
200 const RI::Tensor<Tdata> &&n = cal_overlap_matrix<Tdata>( TA, TB, tauA, tauB, index_A1, index_A2, index_B, Matrix_Order::BA2A1 );
201 matrixes[TB][IB][TA][IA].resize(2);
202 matrixes[TB][IB][TA][IA][1] = std::move(n);
203 }
204 }
205 }
206 }
207 // matrixes[T][I][T][I][0] = matrixes[T][I][T][I][1], so delete repeat
208 for (auto m1 : matrixes)
209 {
210 const size_t T = m1.first;
211 for( auto m2 : m1.second )
212 {
213 const size_t I = m2.first;
214 matrixes[T][I][T][I].resize(1);
215 }
216 }
217
218 return matrixes;
219}
220#endif
Matrix_Orbs21.h
RI_Util.h
i
const std::complex< double > i
Definition cal_pLpR.cpp:46
Atom::na
int na
Definition atom_spec.h:27
Atom::tau
std::vector< ModuleBase::Vector3< double > > tau
Definition atom_spec.h:35
Matrix_Orbs21::center2_orb21_s
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, Center2_Orb::Orb21 > > > > > > > > center2_orb21_s
Definition Matrix_Orbs21.h:86
Matrix_Orbs21::lat0
double * lat0
Definition Matrix_Orbs21.h:76
Matrix_Orbs21::cal_grad_overlap_matrix
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_B, const Matrix_Order &matrix_order) const
Definition Matrix_Orbs21.hpp:90
Matrix_Orbs21::cal_overlap_matrix_all
std::map< size_t, std::map< size_t, std::map< size_t, std::map< size_t, std::vector< 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_B) const
Definition Matrix_Orbs21.hpp:173
Matrix_Orbs21::cal_overlap_matrix
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_B, const Matrix_Order &matrix_order) const
Definition Matrix_Orbs21.hpp:13
Matrix_Orbs21::Matrix_Order
Matrix_Order
Definition Matrix_Orbs21.h:35
Matrix_Orbs21::Matrix_Order::BA2A1
@ BA2A1
Matrix_Orbs21::Matrix_Order::A2A1B
@ A2A1B
Matrix_Orbs21::Matrix_Order::A2BA1
@ A2BA1
Matrix_Orbs21::Matrix_Order::A1BA2
@ A1BA2
Matrix_Orbs21::Matrix_Order::A1A2B
@ A1A2B
Matrix_Orbs21::Matrix_Order::BA1A2
@ BA1A2
ModuleBase::Vector3
3 elements vector
Definition vector3.h:24
UnitCell
Definition unitcell.h:15
UnitCell::atoms
Atom * atoms
Definition unitcell.h:45
T
#define T
Definition exp.cpp:237
ModuleBase::Element_Basis_Index::IndexLNM
std::vector< Index_T > IndexLNM
Definition element_basis_index.h:42
ModuleBase::TITLE
void TITLE(const std::string &class_name, const std::string &function_name, const bool disable)
Definition tool_title.cpp:18
RI_Util::Vector3_to_array3
std::array< Tcell, 3 > Vector3_to_array3(const ModuleBase::Vector3< Tcell > &v)
Definition RI_Util.h:32
Generated by doxygen 1.9.8