Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
3 : ! This file is part of FLEUR and available as free software under the conditions
4 : ! of the MIT license as expressed in the LICENSE file in more detail.
5 : !--------------------------------------------------------------------------------
6 :
7 : module m_wann_projmethod
8 : use m_juDFT
9 : contains
10 0 : subroutine wann_projmethod(
11 : > fullnkpts,
12 : > l_projmethod,l_bestproj,
13 0 : > ikpt,nwfs,nslibd,amn,eig,
14 0 : < psiw,hwfr)
15 : c**********************************************************************
16 : c Construct Wannier functions from the projections of the
17 : c Bloch functions onto the trial orbitals. These projections
18 : c are provided as input to this subroutine in the amn-matrix.
19 : c The Wannier functions are orthonormalized. The orthonormalization
20 : c requires the construction of square root of the overlap
21 : c matrix (omn). Two different algorithms for the calculation of the
22 : c square root are available (projmethod/bestproj).
23 : c
24 : c YM && FF
25 : c**********************************************************************
26 :
27 : use m_types
28 : USE m_constants
29 :
30 : implicit none
31 :
32 : integer, intent(in) :: fullnkpts
33 : logical, intent(in) :: l_projmethod
34 : logical, intent(in) :: l_bestproj
35 : integer, intent(in) :: ikpt
36 : integer, intent(in) :: nwfs
37 : integer, intent(in) :: nslibd
38 : complex, intent(in) :: amn(:,:,:)
39 : real, intent(in) :: eig(:)
40 :
41 : complex, intent(inout) :: psiw(:,:,:)
42 : complex, intent(inout) :: hwfr(:,:)
43 :
44 0 : real :: ei(nwfs)
45 0 : complex, allocatable :: work(:)
46 0 : integer, allocatable :: iwork(:)
47 0 : real, allocatable :: rwork(:)
48 : integer :: info,lwork,lrwork,liwork,lee
49 : integer :: nwf,nwfp,i,n,lda,j
50 0 : complex :: omn(nwfs,nwfs),vec(nwfs,nwfs)
51 0 : complex :: smn(nwfs,nwfs)
52 0 : complex, allocatable :: amn_copy(:,:)
53 : character :: jobu*1,jobv*1
54 : character :: jobz*1,uplo*1
55 0 : complex, allocatable :: sunit(:,:)
56 0 : real, allocatable :: sdiag(:)
57 0 : complex, allocatable :: umn(:,:),vmn(:,:)
58 : integer :: ldu,ldv,m
59 0 : complex :: hwfk(nwfs,nwfs)
60 :
61 0 : call timestart("wann_projmethod")
62 0 : if(l_projmethod) then
63 0 : omn=cmplx(0.0,0.0)
64 : c.. now we determine the Omn matrix, PRB 71,125119, Eq.(13)
65 0 : do nwf = 1,nwfs
66 0 : do nwfp = 1,nwfs
67 0 : do i = 1,nslibd
68 : omn(nwf,nwfp) = omn(nwf,nwfp) +
69 0 : + conjg(amn(i,nwf,ikpt))*amn(i,nwfp,ikpt)
70 : enddo
71 : enddo
72 : enddo
73 :
74 : c lwork = 2*nwfs + nwfs*nwfs
75 : c lrwork = 1 + 4*nwfs + 5*nwfs*nwfs
76 : c liwork = 2 + 6*nwfs
77 0 : lee = log( dble(nwfs) )/log(2.d0) + 1
78 0 : lwork = 1 + 5*nwfs + 2*nwfs*lee + 3*(nwfs**2)
79 0 : lrwork = 1 + 4*nwfs + 2*nwfs*lee + 3*(nwfs**2)
80 0 : liwork = 2 + 5*nwfs +1
81 :
82 0 : allocate (work(lwork),rwork(lrwork),iwork(liwork))
83 :
84 0 : jobz = 'V' ; uplo = 'L' ; n = nwfs ; lda = nwfs
85 :
86 0 : vec=omn
87 :
88 : call zheevd(jobz,uplo,n,vec,lda,ei,work,lwork,
89 0 : & rwork,lrwork,iwork,liwork,info)
90 :
91 0 : if (info.lt.0) write (oUnit,*)
92 0 : & 'ith argument had an illegal value ',info
93 0 : IF (info>0) CALL juDFT_error("not converged diagonalization"
94 0 : + ,calledby ="wann_projmethod")
95 :
96 0 : do i = 1,nwfs
97 0 : if (ei(i).le.(-1.e-6)) then
98 0 : print*,'eigenvalue is less than zero'
99 0 : elseif (abs(ei(i)).le.1.e-16) then
100 0 : print*,'eigenvalue is very close to zero:',ei(i)
101 : endif
102 : enddo
103 :
104 : c.. constructing the square root of the O matrix, Eq.(14)
105 :
106 0 : smn(:,:) = cmplx(0.,0.)
107 :
108 0 : do i = 1,nwfs
109 0 : do j = 1,nwfs
110 0 : do n = 1,nwfs
111 : smn(i,j) = smn(i,j) +
112 0 : + conjg(vec(i,n))*(1./sqrt(ei(n)))*vec(j,n)
113 : enddo
114 : enddo
115 : enddo
116 :
117 0 : deallocate (work,rwork,iwork)
118 : c.. now constructing the overlaps of the wavefunctions
119 : c.. with the wannier functions in reciprocal space, Eq.(16)
120 :
121 0 : psiw(:,:,ikpt) = cmplx(0.,0.)
122 0 : do i = 1,nslibd
123 0 : do j = 1,nwfs
124 0 : do n = 1,nwfs
125 : psiw(i,j,ikpt) = psiw(i,j,ikpt) +
126 0 : + smn(j,n)*amn(i,n,ikpt)
127 : enddo
128 : enddo
129 : enddo
130 : endif !wann%l_projmethod
131 :
132 0 : if(l_bestproj) then
133 :
134 : c.. psiw matrix has the meaning of the Umn rotational matrix
135 : c.. applied to the wavefunctions at each point in the BZone
136 : c.. here it is calculated differently,
137 : c.. according to the PRB 65,035109 (Eq.23)
138 : c..
139 : c.. the singe value decomposition of the Amn matrix is found:
140 : c.. A = U*SS*V
141 : c.. Then the psiw matrix, psiw = Amn*Smn is given by:
142 : c.. psiw = U*1*V
143 :
144 0 : allocate (amn_copy(nslibd,nwfs),sunit(nslibd,nwfs))
145 0 : allocate (umn(nslibd,nslibd),vmn(nwfs,nwfs))
146 0 : allocate (sdiag(max(nslibd,nwfs)))
147 :
148 : lwork = max( 3*min(nslibd,nwfs) + max(nslibd,nwfs),
149 0 : & 5*min(nslibd,nwfs) )
150 :
151 0 : allocate ( work(lwork),rwork(max(1,5*min(nslibd,nwfs))) )
152 : !
153 : ! Compute the eigenvalues and eigenvectors.
154 : !
155 0 : jobu = 'A' ; jobv = 'A'
156 0 : lda = nslibd ; ldu = nslibd ; ldv = nwfs
157 : !
158 : ! The input matrix is destroyed by the routine. Since we need to keep
159 : ! it around, we only pass a copy to the routine.
160 : !
161 0 : amn_copy(1:nslibd,1:nwfs) = amn(1:nslibd,1:nwfs,ikpt)
162 :
163 : call zgesvd(jobu,jobv,nslibd,nwfs,amn_copy,lda,
164 0 : > sdiag,umn,ldu,vmn,ldv,work,lwork,rwork,info)
165 :
166 0 : if (info /= 0) then
167 0 : write (*,*) 'LAPACK routine ZGESVD returned a nonzero'
168 0 : write (*,*) 'value of the error flag, INFO =',info
169 : end if
170 : !
171 : ! Make the MxN identical Sunit matrix
172 : !
173 0 : sunit(1:nslibd,1:nwfs) = cmplx(0.,0.)
174 0 : do i = 1,min(nslibd,nwfs)
175 0 : sunit(i,i) = cmplx(1.,0.)
176 : end do
177 :
178 :
179 : ! and finally the psiw matrix
180 :
181 0 : psiw(:,:,ikpt) = cmplx(0.,0.)
182 0 : do i = 1,nslibd
183 0 : do j = 1,nwfs
184 0 : do n = 1,nslibd
185 0 : do m = 1,nwfs
186 : psiw(i,j,ikpt) = psiw(i,j,ikpt) +
187 0 : + umn(i,n)*sunit(n,m)*vmn(m,j)
188 : enddo
189 : enddo
190 : enddo
191 : enddo
192 :
193 0 : deallocate (work,rwork,amn_copy,sunit,vmn,umn,sdiag)
194 :
195 0 : write (*,*) 'The Psiw matrix was calculated via SVD'
196 :
197 : endif !wann%l_bestproj
198 :
199 :
200 :
201 : c..constructing the WF-hamiltonian in reciprocal space Eq.(23)
202 :
203 0 : hwfk(:,:) = cmplx(0.,0.)
204 :
205 0 : do i = 1,nwfs
206 0 : do j = 1,nwfs
207 0 : do n = 1,nslibd
208 : hwfk(i,j) = hwfk(i,j) +
209 0 : + eig(n)*psiw(n,i,ikpt)*conjg(psiw(n,j,ikpt))
210 :
211 : enddo
212 : enddo
213 : enddo
214 :
215 : c.. adding up the k-point reciprocal hamiltonians Eq.(20)
216 : c.. to get the hopping elements inside the same unit cell
217 :
218 0 : hwfr=hwfr+hwfk/fullnkpts
219 :
220 :
221 : c.. now we diagonalize the reciprocal hamiltonian for the
222 : c.. purpose of testing
223 :
224 0 : do nwf = 1,nwfs
225 0 : do nwfp = 1,nwfs
226 0 : vec(nwf,nwfp) = hwfk(nwf,nwfp)
227 : enddo
228 : enddo
229 :
230 0 : lee = log( dble(nwfs) )/log(2.d0) + 1
231 0 : lwork = 1 + 5*nwfs + 2*nwfs*lee + 3*(nwfs**2)
232 0 : lrwork = 1 + 4*nwfs + 2*nwfs*lee + 3*(nwfs**2)
233 0 : liwork = 2 + 5*nwfs +1
234 :
235 0 : allocate (work(lwork),rwork(lrwork),iwork(liwork))
236 :
237 0 : jobz = 'V' ; uplo = 'L' ; n = nwfs ; lda = nwfs
238 :
239 : call zheevd(jobz,uplo,n,vec,lda,ei,work,lwork,
240 0 : & rwork,lrwork,iwork,liwork,info)
241 :
242 0 : if (info.lt.0) write (oUnit,*)
243 0 : & 'ith argument had an illegal value ',info
244 0 : IF (info>0) CALL juDFT_error("not converged diagonalization"
245 0 : + ,calledby ="wann_projmethod")
246 :
247 0 : do i = 1,nwfs
248 0 : write(oUnit,*) ei(i)
249 : enddo
250 :
251 0 : deallocate (work,rwork,iwork)
252 0 : call timestop("wann_projmethod")
253 :
254 0 : end subroutine wann_projmethod
255 : end module m_wann_projmethod
|
