-
Notifications
You must be signed in to change notification settings - Fork 0
/
hf-bloch.cpp
207 lines (174 loc) · 5.95 KB
/
hf-bloch.cpp
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
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
#include <iostream>
#include <cmath>
#include <fstream>
#include <algorithm>
#include <complex>
#include <lapacke.h>
#include <Eigen/Dense>
using namespace std;
using namespace Eigen;
using namespace std::literals;
typedef complex<double> cd;
typedef vector < pair<double,VectorXcd> > spectrum;
int no_of_unitcell_pts=20;
double a = 1;
double dx = a/double(no_of_unitcell_pts);
int N = 20; //no of unit cells
int no_of_total_pts = no_of_unitcell_pts*N;
double init_pt= 0.0;
double final_pt = 10.0;
double A=1;
double epsilon = dx/2;
spectrum* arr = new spectrum [N];
VectorXd unitcell_point;
double V(double x){return -A*(1+cos(2*M_PI*x/a));}
bool compare(const pair<double, VectorXcd>&i, const pair<double, VectorXcd>&j) {return i.first < j.first;}
cd u(int m, int k, int i) { return arr[k].at(m).second(i%no_of_unitcell_pts);}
bool diagonalize(MatrixXcd A, VectorXcd& lambda, MatrixXcd& v)
{
int N = A.cols();
if (A.rows()!=N) return false;
v.resize(N,N);
lambda.resize(N);
int LDA = A.outerStride();
int LDV = v.outerStride();
int INFO = 0;
cd* w = const_cast<cd*>(lambda.data());
char Nchar = 'N';
char Vchar = 'V';
int LWORK = int(A.size())*10;
VectorXcd WORK(LWORK);
VectorXd RWORK(5*LDA);
zgeev_(&Nchar, &Vchar, &N, reinterpret_cast <__complex__ double*> (A.data()), &LDA, reinterpret_cast <__complex__ double*> (w), 0, &LDV, reinterpret_cast <__complex__ double*> (v.data()), &LDV, reinterpret_cast <__complex__ double*> (WORK.data()), &LWORK, RWORK.data(), &INFO);
return INFO==0;
}
double direct_core_integrand(int m, int k, int kappa, int i, int i_prime)
{
return 1/(4*M_PI)*norm(u(m,k,i))*norm(u(m,kappa,i_prime))/(abs(unitcell_point(i)-unitcell_point(i_prime))+epsilon);
}
double direct_integral(int m, int i, int k, int kappa)
{
double trapez_sum=0.0;
for(int i_prime=0; i_prime < unitcell_point.size(); i_prime++)
trapez_sum+= direct_core_integrand(m,k,kappa,i,i_prime);
return trapez_sum;
}
double vh(int m, int k, double i)
{
double vh = 0;
for(int kappa=0; kappa<N; kappa++)
{
if(k==kappa) continue;
vh += direct_integral(m,i,k,kappa);
}
return vh;
}
cd exchange_core_integrand(int m, int k, int kappa, int i, int i_prime)
{
cd num = conj(u(m,kappa,i))*u(m,kappa,i_prime)*conj(u(m,k,i_prime))*u(m,k,i)*exp(cd(0,(k-kappa)*(unitcell_point(i)-unitcell_point(i_prime))));
cd denom = (abs(unitcell_point(i)-unitcell_point(i_prime))+epsilon);
return 1/(4*M_PI)*num/denom;
}
cd exchange_integral(int m, int i, int k, int kappa)
{
cd trapez_sum=cd(0,0);
for(int i_prime=0; i_prime<unitcell_point.size(); i_prime++)
trapez_sum+= exchange_core_integrand(m,k,kappa,i,i_prime);
return trapez_sum;
}
double vhf(int m, int k, double i)
{
double vhf = 0;
for(int kappa=0; kappa<N; kappa++)
{
if(k==kappa) continue;
vhf += exchange_integral(m,i,k,kappa).real();
}
return vhf;
}
VectorXd sortascending(VectorXd v1)
{
vector<double> stdv1 (v1.data(),v1.data()+v1.size());
std::sort (stdv1.begin(), stdv1.end());
stdv1.resize(5);
Map<ArrayXd> sorted(stdv1.data(), stdv1.size());
return sorted;
}
int main()
{
unitcell_point.resize(no_of_unitcell_pts);
for(int i=0; i<unitcell_point.size(); i++) unitcell_point(i) = i*dx;
MatrixXcd H(no_of_unitcell_pts,no_of_unitcell_pts);
spectrum eigenspectrum;
VectorXcd v; MatrixXcd eigenvectors; VectorXd eigenvalues;
for(int n=0; n<N; n++)
{
double k = (2*M_PI*n)/(N*a)-M_PI/a;
for(int i=0; i<unitcell_point.size(); i++)
{
int j = (i==unitcell_point.size()-1)? 0 : i+1;
H(i,j)= cd(-1/(2*dx*dx), -k/(2*dx));
H(j,i)= cd(-1/(2*dx*dx), k/(2*dx));
H(i,i)= cd(1/(dx*dx)+ pow(k,2)/2+ V(unitcell_point(i)), 0);
}
diagonalize(H,v,eigenvectors);
for(int i=0; i<unitcell_point.size(); i++) eigenspectrum.push_back(make_pair(v[i].real(), eigenvectors.col(i)));
sort(eigenspectrum.begin(),eigenspectrum.end(),compare);
eigenspectrum.resize(N);
arr[n] = eigenspectrum;
eigenspectrum.clear();
}
spectrum* nonint_arr = new spectrum [N];
for(int i=0; i<N; i++) nonint_arr[i] = arr[i];
spectrum int_eigenspectrum;
char choice_for_result = 'n';
int master_loop = 1;
spectrum* int_arr = new spectrum [N];
for(int m=0; m<2; m++)
{
std::cout << "Band:" << m << "\n===================================\n" << '\n';
for( ; ; )
{
std::cout << "Loop:" << master_loop << "\n------------------------------------\n" << '\n';
for(int kindex=0; kindex<N; kindex++)
{
double k = (2*M_PI*kindex)/(N*a)-M_PI/a;
for(int i=0; i<unitcell_point.size(); i++)
{
int j = (i==unitcell_point.size()-1)? 0 : i+1;
H(i,j)= cd(-1/(2*dx*dx), -k/(2*dx));
H(j,i)= cd(-1/(2*dx*dx), k/(2*dx));
H(i,i)= cd(1/(dx*dx)+ pow(k,2)/2+ V(unitcell_point(i)), 0)+2*vh(m,kindex,i)-vhf(m,kindex,i);
}
diagonalize(H,v,eigenvectors);
for(int i=0; i<unitcell_point.size(); i++) int_eigenspectrum.push_back(make_pair(v[i].real(), eigenvectors.col(i)));
sort(int_eigenspectrum.begin(),int_eigenspectrum.end(),compare);
int_eigenspectrum.resize(N);
int_arr[kindex] = int_eigenspectrum;
int_eigenspectrum.clear();
}
VectorXd diff(N);
for(int i=0; i<N; i++) diff(i) = abs(int_arr[i].at(m).first-arr[i].at(m).first);
if(choice_for_result=='y')
{
for(int i=0; i<N; i++)
cout << i << " " << arr[i].at(m).first << " " << int_arr[i].at(m).first << " " << diff(i) << endl;
cout << endl;
}
for(int i=0; i<N; i++) arr[i] = int_arr[i];
cout << "Max diff= " << diff.maxCoeff() << endl << endl;
double tolerance = (m==0)?1e-3:0.003;
if(diff.maxCoeff()<tolerance) break;
master_loop++;
}
}
ofstream dataout("hf_data.txt");
for(int i=0; i<N; i++)
{
dataout << (2*M_PI*i)/(N*a)-M_PI/a << " ";
for(int m=0; m<2; m++)
dataout << nonint_arr[i].at(m).first << " " << arr[i].at(m).first << " ";
dataout << endl;
}
return 0;
}