-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathcbs.f90
221 lines (182 loc) · 6.29 KB
/
cbs.f90
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
208
209
210
211
212
213
214
215
216
217
218
219
220
!main of the program
program cbs
use readers, only: read_hsx_file, hsx_t, read_OrbIndx_file
use FDF
use calculators
implicit none
logical :: Ksetting, SymSetting
double precision :: EMIN, EMAX,p1,p2,dist
type(hsx_t) :: hw
character(len=10) :: pre
integer :: bins,j,l,maxsuperc(3),dir,norb,nspin,ooo
integer,allocatable :: isc(:,:)
double precision :: k(3),a1(3),a2(3),a3(3),p(2),ap(3),b1(3),b2(3)
!Reading cbs.in, mandatory call it cbs.in
CALL FDF_INIT("cbs.in","out")
Ksetting = FDF_BOOLEAN('VaringKsetting',.FALSE.)
SymSetting = FDF_BOOLEAN('SymSetting',.FALSE.)
EMIN = FDF_PHYSICAL('CBS.Emin',-1.0d10,'eV')
EMAX = FDF_PHYSICAL('CBS.Emax',1.0d10,'eV')
bins = FDF_INTEGER('CBS.bins',100)
dir = FDF_INTEGER('CBS.dir',1)
pre = FDF_STRING('SystemLabel','siesta')
p(1) = FDF_GET('CBS.2DBZk1',0d0)
p(2) = FDF_GET('CBS.2DBZk2',0d0)
!Reading .HSX file
call read_hsx_file(hw,trim(adjustl(pre))//".HSX")
!Reading .ORB_INDX file
allocate(isc(hw%no_s,3))
call read_OrbIndx_file(isc,trim(adjustl(pre))//".ORB_INDX")
!Reading .XV file
!In bohr
open(189,file=trim(adjustl(pre))//".XV",status="old")
read(189,*) a1(1),a1(2),a1(3)
read(189,*) a2(1),a2(2),a2(3)
read(189,*) a3(1),a3(2),a3(3)
close(189)
!Summary on the direction you calculate the complex k
!calculate the 2D reciprocal vectors!
!from ap you also get the equation of the plane containing two lattice vector
!and then the distance from the third vector
if (dir==1) then
write(*,*) "You are calculating the complex band structure in the direction perpendicular"
write(*,*) "to the plane generated by lattice vectors a2 and a3 (in next lines in bohr)"
write(*,*) "a2=",a2
write(*,*) "a3=",a3
call get_2D_rec_vect(a2,a3,b1,b2,ap)
write(*,*) "In cartesian units this direction is"
write(*,*) ap
dist=abs(ap(1)*a1(1)+ap(2)*a1(2)+ap(3)*a1(3))/sqrt(ap(1)**2+ap(2)**2+ap(3)**2)
write(*,*) "Height of the layer in the direction under scuriny (bohr):"
write(*,*) dist
write(*,*) "K of the CBS is normalized dividing by the number above"
else if (dir==2) then
write(*,*) "You are calculating the complex band structure in the direction perpendicular"
write(*,*) "to the plane generated by lattice vectors a3 and a1 (in next lines in bohr)"
write(*,*) "a3=",a3
write(*,*) "a1=",a1
call get_2D_rec_vect(a3,a1,b1,b2,ap)
write(*,*) "In cartesian units this direction is"
write(*,*) ap
dist=abs(ap(1)*a2(1)+ap(2)*a2(2)+ap(3)*a2(3))/sqrt(ap(1)**2+ap(2)**2+ap(3)**2)
write(*,*) "Height of the layer in the direction under scuriny (bohr):"
write(*,*) dist
write(*,*) "K of the CBS is normalized dividing by the number above"
elseif (dir==3) then
write(*,*) "You are calculating the complex band structure in the direction perpendicular"
write(*,*) "to the plane generated by lattice vectors a1 and a2 (in next lines in bohr)"
write(*,*) "a1=",a1
write(*,*) "a2=",a2
call get_2D_rec_vect(a1,a2,b1,b2,ap)
write(*,*) "In cartesian units this direction is"
write(*,*) ap
dist=abs(ap(1)*a3(1)+ap(2)*a3(2)+ap(3)*a3(3))/sqrt(ap(1)**2+ap(2)**2+ap(3)**2)
write(*,*) "Height of the layer in the direction under scuriny (bohr):"
write(*,*) dist
write(*,*) "K of the CBS is normalized dividing by the number above"
else
STOP "Error: you messed up with direction, allowed values are 1, 2 or 3"
endif
!Norbital and Nspin
norb=hw%no_u
nspin=hw%nspin
!Understanding how many replicas of the elementary cell
!are in the supercell in each direction
maxsuperc(:)=1
do j=1,3
do l=1,hw%no_s
if (isc(l,j).gt.maxsuperc(j)) then
maxsuperc(j)=isc(l,j)
endif
enddo
enddo
!Last summary on the calculation properties
write(*,*) "N cell replicas in choosen direc",maxsuperc(dir), "norb",hw%no_u, "nspin", hw%nspin
if (Ksetting.and.SymSetting) then
STOP "Error: VaringKsetting OR SymSetting is allowed, can't do both at same time"
endif
if (Ksetting) then
write(*,*) "b1 and b2 are the 2D reciprocal vector:"
write(*,*) "b1=",b1(:)
write(*,*) "b2=",b2(:)
write(*,*) "Varing k parallel for E =", emax
call CBSvariableK(hw,dir,isc,emax,dist,b1,b2)
else if (SymSetting) then
k(:)=p(1)*b1(:)+p(2)*b2(:)
write(*,*) "Analysing projections at"
write(*,*) "E =", emax
write(*,*) "K =", k
call CBSsymmetries(hw,dir,isc,k,emax,dist)
else
!Define the k parallel where to calculate the CBS and print out a summary
k(:)=p(1)*b1(:)+p(2)*b2(:)
write(*,*) "The k parallel you choose is"
write(*,*) p(1),"*b1 + ",p(2),"*b2"
write(*,*) "where b1 and b2 are reciprocal lattice vectors of the 2D BZ"
write(*,*) "b1 and b2 in the next lines:"
write(*,*) "b1=",b1(:)
write(*,*) "b2=",b2(:)
write(*,*) "The k parallel you choose in cartesian units (bohr) is:"
write(*,*) k(:)
call CBSfixedK(hw,dir,isc,k,emin,emax,bins,dist)
endif
deallocate(isc)
end program
subroutine get_2D_rec_vect(aa,aabis,b1,b2,ap)
!The theory behind:
!I first calculate the vector perpendicular to the aa and aabis,
!I call it ap, second I normalise it.
!Third I solve the system aa*b1=2pi,aabis*b1=0,ap*b1=0
!which is the definition of reciprocal vectors and to impose b1
!is in the same plane of aa and aabis
!Same for b2
implicit none
double precision, intent(in) :: aa(3),aabis(3)
double precision, intent(out) :: b1(3),b2(3),ap(3)
integer :: i,info
double precision :: Matrix(3,3), OtherMat(3,3),OtherMat2(3,3),IPIV(3,3)
double precision :: PI,apmod,det,bisdet,trisdet
ap(:)=0.d0
ap(1)=aa(2)*aabis(3)-aa(3)*aabis(2)
ap(2)=aabis(1)*aa(3)-aa(1)*aabis(3)
ap(3)=aa(1)*aabis(2)-aa(2)*aabis(1)
apmod=sqrt(ap(1)**2+ap(2)**2+ap(3)**2)
ap(:)=ap(:)/apmod
Matrix(1,:)=ap(:)
Matrix(2,:)=aa(:)
Matrix(3,:)=aabis(:)
det=0.d0
bisdet=0.d0
trisdet=0.d0
call a33det(Matrix, det)
b1(:)=0d0
b2(:)=0d0
PI=4.D0*DATAN(1.D0)
do i=1,3
OtherMat=Matrix
OtherMat(1,i)=0
OtherMat(2,i)=0
OtherMat(3,i)=2*PI
OtherMat2=Matrix
OtherMat2(1,i)=0
OtherMat2(2,i)=2*PI
OtherMat2(3,i)=0
call a33det(OtherMat, bisdet)
call a33det(OtherMat2, trisdet)
b1(i)=bisdet/det
b2(i)=trisdet/det
enddo
return
end subroutine
subroutine a33det(A,DET)
implicit none
double precision, intent(in) :: A(3,3)
double precision, intent(out) :: DET
DET = A(1,1)*A(2,2)*A(3,3) &
- A(1,1)*A(2,3)*A(3,2) &
- A(1,2)*A(2,1)*A(3,3) &
+ A(1,2)*A(2,3)*A(3,1) &
+ A(1,3)*A(2,1)*A(3,2) &
- A(1,3)*A(2,2)*A(3,1)
return
end subroutine a33det