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