forked from abacusmodeling/abacus-develop
-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathcenter2_orb-orb22.cpp
More file actions
130 lines (110 loc) · 3.22 KB
/
center2_orb-orb22.cpp
File metadata and controls
130 lines (110 loc) · 3.22 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
//=========================================================
//AUTHOR : Peize Lin
//DATE : 2016-01-24
//=========================================================
#include "center2_orb-orb22.h"
Center2_Orb::Orb22::Orb22(
const Numerical_Orbital_Lm &nA1_in,
const Numerical_Orbital_Lm &nA2_in,
const Numerical_Orbital_Lm &nB1_in,
const Numerical_Orbital_Lm &nB2_in,
const ORB_table_phi &MOT_in,
const ORB_gaunt_table &MGT_in
)
:nA1(nA1_in),
nA2(nA2_in),
nB1(nB1_in),
nB2(nB2_in),
MOT(MOT_in),
MGT(MGT_in){}
void Center2_Orb::Orb22::init_radial_table()
{
const Numerical_Orbital_Lm & nB_short = (nB1.getNr()<=nB2.getNr()) ? nB1 : nB2;
std::vector<double> nB_tmp(nB_short.getNr());
for( size_t ir=0; ir!=nB_tmp.size(); ++ir)
{
nB_tmp[ir] = nB1.getPsi(ir) * nB2.getPsi(ir);
}
const int LB1 = nB1.getL();
const int LB2 = nB2.getL();
for( int LB = abs(LB1-LB2); LB<=LB1+LB2; ++LB)
{
if( (LB-std::abs(LB1-LB2))%2==1 ) // if LA+LB-LAB == odd, then Gaunt_Coefficients = 0
continue;
nB[LB].set_orbital_info(
nB_short.getLabel(),
nB_short.getType(),
LB,
1, // N?
nB_short.getNr(),
nB_short.getRab(),
nB_short.getRadial(),
Numerical_Orbital_Lm::Psi_Type::Psi,
ModuleBase::GlobalFunc::VECTOR_TO_PTR(nB_tmp),
nB_short.getNk(),
nB_short.getDk(),
nB_short.getDruniform(),
false,
true,
GlobalV::CAL_FORCE); // mohan add 2021-05-07
orb21s.insert( make_pair( LB, Center2_Orb::Orb21( nA1, nA2, nB[LB], MOT, MGT ) ) );
orb21s.at(LB).init_radial_table();
}
}
void Center2_Orb::Orb22::init_radial_table( const std::set<size_t> &radials )
{
const Numerical_Orbital_Lm & nB_short = (nB1.getNr()<=nB2.getNr()) ? nB1 : nB2;
std::vector<double> nB_tmp(nB_short.getNr());
for( size_t ir=0; ir!=nB_tmp.size(); ++ir)
{
nB_tmp[ir] = nB1.getPsi(ir) * nB2.getPsi(ir);
}
const int LB1 = nB1.getL();
const int LB2 = nB2.getL();
for( int LB = abs(LB1-LB2); LB<=LB1+LB2; ++LB)
{
if( (LB-std::abs(LB1-LB2))%2==1 ) // if LA+LB-LAB == odd, then Gaunt_Coefficients = 0
continue;
nB[LB].set_orbital_info(
nB_short.getLabel(),
nB_short.getType(),
LB,
1, // N?
nB_short.getNr(),
nB_short.getRab(),
nB_short.getRadial(),
Numerical_Orbital_Lm::Psi_Type::Psi,
ModuleBase::GlobalFunc::VECTOR_TO_PTR(nB_tmp),
nB_short.getNk(),
nB_short.getDk(),
nB_short.getDruniform(),
false,
true, GlobalV::CAL_FORCE);
orb21s.insert( make_pair( LB, Center2_Orb::Orb21( nA1, nA2, nB[LB], MOT, MGT ) ) );
orb21s.at(LB).init_radial_table(radials);
}
}
double Center2_Orb::Orb22::cal_overlap(
const ModuleBase::Vector3<double> &RA, const ModuleBase::Vector3<double> &RB,
const int &mA1, const int &mA2, const int &mB1, const int &mB2) const
{
const int LB1 = nB1.getL();
const int LB2 = nB2.getL();
double overlap = 0.0;
for( const auto& orb21 : orb21s )
{
const int LB = orb21.first;
for( int mB=0; mB!=2*LB+1; ++mB )
// const int mB=mB1+mB2;
{
const double Gaunt_real_B1_B2_B12 =
MGT.Gaunt_Coefficients (
MGT.get_lm_index(LB1,mB1),
MGT.get_lm_index(LB2,mB2),
MGT.get_lm_index(LB,mB));
if( 0==Gaunt_real_B1_B2_B12 ) continue;
overlap += Gaunt_real_B1_B2_B12 * orb21.second.cal_overlap(RA, RB, mA1, mA2, mB);
}
}
return overlap;
}