-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathSurfacePhononGF.2.1
129 lines (115 loc) · 4.39 KB
/
SurfacePhononGF.2.1
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
#include "SurfacePhononGF.h"
#include "FuncUtils.h"
//#include "constants.h"
#include <fstream>
#include <iomanip>
using std::fstream;
void surfphGF::DR00(v4cd& d00r){
std::cout << string(45,'*')<<std::endl;
std::cout << "Solve surface phonon green's function iteratively." << std::endl;
std::cout << string(45,'*')<<std::endl;
int nk = mfcbylz.size();
int no = omgv.size();
// std::cout << "surfphGF: omgv: " << std::endl;
// for(auto i : omgv) std::cout << i << std::endl;
// std::cout << "no " << no << std::endl;
// std::cout << "nk " << nk << std::endl;
int nx = mfcbylz[0][0][0].size();
// std::cout << "nx: " << nx << std::endl;
for(int k = 0; k < nk; k++){
MatrixXcd K00(nx, nx);
MatrixXcd K01(nx, nx);
MatrixXcd K10(nx, nx);
for(int ia = 0; ia < nx; ia++){
for(int ja = 0; ja < nx; ja++){
K00(ia, ja) = mfcbylz[k][0][0][ia][ja];
K01(ia, ja) = mfcbylz[k][0][1][ia][ja];
K10(ia, ja) = mfcbylz[k][1][0][ia][ja];
}
}
for(int o = 0; o < no; o++){
MatrixXcd K0 = K00;
MatrixXcd K0s = K00;
MatrixXcd K1 = K01;
MatrixXcd K2 = K10;
complex<double> dtli = complex<double>(omgv[o], delta);
std::cout << "omega: " << omgv[o] << std::endl;
MatrixXcd ComplexIden = Eigen::MatrixXcd::Identity(nx, nx) * dtli * dtli;
DeciTech(K0, K0s, K1, K2, ComplexIden);
MatrixXcd d00rM = (ComplexIden - K0s).householderQr().solve(Eigen::MatrixXcd::Identity(nx, nx));
// d00rM = K01 * d00rM * K10;
for(int ia = 0; ia < nx; ia++){
for(int ja = 0; ja < nx; ja++){
d00r[k][o][ia][ja] = d00rM(ia, ja);
}
}
// if(no%(o+1) == 0)
// std::cout << std::fixed<<std::setprecision(0) << o*1.0/(no-1)*100 << " % completed." << std::endl;
}
}
}
// Periodic Pulay method
void surfphGF::DeciTech(MatrixXcd& K0, MatrixXcd& K0s,
MatrixXcd& K1, MatrixXcd& K2, MatrixXcd Iden)
{
const size_t nx = Iden.cols();
MatrixXcd d00 = (Iden - K0).householderQr().solve(Eigen::MatrixXcd::Identity(nx, nx));
MatrixXcd K0sNext = K1 * d00 * K2;
double meps = K0sNext.norm();
if(meps > 1.0e10){
std::cerr << "Surface Green function can not converge to a finite vaule." << std::endl;
abort();
}
if(meps < 1.0e-8) return;
// std::cout << "K0sNext: " << K0sNext << std::endl;
K0s = K0s + 0.05 * K0sNext;
K0 = K0 + 0.05 * (K2 * d00 * K1 + K1 * d00 * K2);
K1 = K1 + 0.05*(K1 * d00 * K1 - K1);
K2 = K2 + 0.05*(K2 * d00 * K2 - K2);
DeciTech(K0, K0s, K1, K2, Iden);
}
void surfphGF::DR00(MatrixXcd t0, MatrixXcd t0t, MatrixXcd Tphnext, MatrixXcd& Tph){
MatrixXcd Iden = MatrixXcd::Identity(t0.rows(), t0.cols());
// MatrixXcd t0next = 0.01*((Iden - t0*t0t - t0t*t0).fullPivLu().inverse()*t0*t0-t0);
// MatrixXcd t0tnext = 0.01*((Iden - t0*t0t - t0t*t0).fullPivLu().inverse()*t0t*t0t-t0t);
MatrixXcd t0next = (Iden - t0*t0t - t0t*t0).householderQr().solve(Iden) * t0 * t0;
MatrixXcd t0tnext = (Iden - t0*t0t - t0t*t0).householderQr().solve(Iden) * t0t * t0t;
MatrixXcd Tphtnext = Tphnext*t0tnext;
std::cout << "Tphtnext: "<< std::endl;
std::cout << Tphtnext << std::endl;
if(Tphnext.norm() < 1.0e-8) return;
else{
Tph += Tphtnext;
DR00(t0next, t0tnext, Tphtnext, Tph);
}
}
void surfphGF::writeDR00(string ssgf, double cellLz, int nlayers,
const vector<vector<double> >& kp2d,
const v4cd& d00r){
fstream ofs(ssgf, std::ios::out|std::ios::binary);
ofs.write((char*)&cellLz, sizeof(double));
ofs.write((char*)&nlayers, sizeof(int));
int kpsz1 = kp2d.size();
int kpsz2 = kp2d[0].size();
ofs.write((char*)&kpsz1, sizeof(int));
ofs.write((char*)&kpsz2, sizeof(int));
for(int i = 0; i < kpsz1; i++)
ofs.write((char*)&kp2d[i][0], kpsz2*sizeof(double));
int dsz1 = d00r.size();
int dsz2 = d00r[0].size();
int dsz3 = d00r[0][0].size();
int dsz4 = d00r[0][0][0].size();
// std::cout << "d00r size: "<<dsz1<<"*"<<dsz2<<"*"<<dsz3<<"*"<<dsz4<<"*"<<std::endl;
ofs.write((char*)&dsz1, sizeof(int));
ofs.write((char*)&dsz2, sizeof(int));
ofs.write((char*)&dsz3, sizeof(int));
ofs.write((char*)&dsz4, sizeof(int));
for(int i = 0; i < dsz1; i++)
for(int j = 0; j < dsz2; j++)
for(int k = 0; k < dsz3; k++)
ofs.write((char*)(&d00r[i][j][k][0]), dsz4*sizeof(complex<double>));
// double kp2da[kp2d.size()][kp2d[0].size()];
// mcopy(kp2d, (double*)kp2da);
// ofs.write(reinterpret_cast<char*>(kp2da), sizeof(kp2da));
ofs.close();
}