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