1!
2! Copyright (C) 2003 A. Smogunov
3! This file is distributed under the terms of the
4! GNU General Public License. See the file `License'
5! in the root directory of the present distribution,
6! or http://www.gnu.org/copyleft/gpl.txt .
7!
8! Generalized to spinor wavefunctions and spin-orbit Oct. 2004 (ADC).
9!
10!
11subroutine compbs(lleft, nocros, norb, nchan, kval, kfun,  &
12                  kfund, kint, kcoef, ikk, ien)
13!
14! Using the basis functions obtained by scatter_forw it computes
15! the complex band structure (CBS) of the lead.
16! Some variables needed for wave-function matching in transmission
17! calculation are constructed and saved.
18!
19  USE constants, ONLY : tpi
20  USE noncollin_module, ONLY : noncolin, npol
21  USE spin_orb,  ONLY : lspinorb
22  USE lsda_mod,  ONLY : nspin
23  USE cond
24  USE cell_base, ONLY : alat, at, omega
25  USE ions_base, ONLY : nat, ityp
26
27  implicit none
28  integer ::   &
29     nocros,   & ! number of orbitals crossing the boundary
30     noins,    & ! number of interior orbitals
31     norb,     & ! total number of orbitals
32     lleft       ! 1/0 if it is left/right tip
33  integer :: ik, ikk, i, j, ig, n, iorb, iorb1, &
34             iorb2, aorb, borb, nchan,     &
35             ij, is, js, ichan, ien
36  REAL(DP), PARAMETER :: eps=1.d-8
37  REAL(DP) :: raux, ddot
38  REAL(DP), ALLOCATABLE :: zpseu(:,:,:), zps(:,:)
39  COMPLEX(DP), PARAMETER :: cim=(0.d0,1.d0)
40  COMPLEX(DP) :: x1,          &
41            kval(2*(n2d+npol*nocros)), kfun(n2d,2*(n2d+npol*nocros)),  &
42            kint(nocros*npol,2*(n2d+npol*nocros)),                     &
43            kcoef(nocros*npol,2*(n2d+npol*nocros)),                    &
44            kfund(n2d,2*(n2d+npol*nocros))
45  COMPLEX(DP), ALLOCATABLE :: amat(:,:), bmat(:,:), vec(:,:), &
46                              zpseu_nc(:,:,:,:), zps_nc(:,:), &
47                              aux(:,:), veceig(:,:), korb(:,:)
48  COMPLEX(DP), PARAMETER :: one=(1.d0,0.d0), zero=(0.d0,0.d0)
49
50
51  call start_clock('compbs')
52
53  noins = norb-2*nocros
54  IF (norb>0) THEN
55     if(lleft.eq.1) then
56       if (noncolin) then
57         allocate( zpseu_nc(2,norb,norb,nspin) )
58         zpseu_nc = zpseul_nc
59       else
60         allocate( zpseu(2,norb,norb) )
61         zpseu = zpseul
62       endif
63     else
64       if (noncolin) then
65          allocate( zpseu_nc(2,norb,norb,nspin) )
66          zpseu_nc = zpseur_nc
67       else
68          allocate( zpseu(2,norb,norb) )
69          zpseu = zpseur
70       endif
71     endif
72     if (noncolin) then
73        allocate( zps_nc( norb*npol, norb*npol ) )
74     else
75        allocate( zps( norb, norb ) )
76     endif
77  END IF
78
79  allocate( amat( (2*n2d+npol*norb), (2*n2d+npol*norb) ) )
80  allocate( bmat( (2*n2d+npol*norb), (2*n2d+npol*norb) ) )
81  allocate( vec( (2*n2d+npol*norb), 2*(n2d+npol*nocros) ) )
82  allocate( aux( n2d, 2*n2d+npol*norb))
83  IF (lorb) allocate( korb(npol*(nocros+noins),2*(n2d+npol*nocros)) )
84
85  amat=(0.d0,0.d0)
86  bmat=(0.d0,0.d0)
87
88!
89! zps=zpseu-e*qq for US-PP and zps=zpseu for norm-conserv. PP
90!
91  do iorb=1, norb
92    do iorb1=1, norb
93      if (noncolin) then
94        ij=0
95        do is=1,npol
96          do js=1,npol
97            ij=ij+1
98            zps_nc(npol*(iorb-1)+is, npol*(iorb1-1)+js)=         &
99                zpseu_nc(1,iorb,iorb1,ij)
100            if (lspinorb) then
101              zps_nc(npol*(iorb-1)+is,npol*(iorb1-1)+js)=        &
102                    zps_nc(npol*(iorb-1)+is,npol*(iorb1-1)+js) &
103                            -eryd*zpseu_nc(2,iorb,iorb1,ij)
104            else
105              if (is.eq.js)                                      &
106                  zps_nc(npol*(iorb-1)+is,npol*(iorb1-1)+js)=    &
107                      zps_nc(npol*(iorb-1)+is,npol*(iorb1-1)+js) &
108                             -eryd*zpseu_nc(2,iorb,iorb1,ij)
109            endif
110          enddo
111        enddo
112      else
113        zps(iorb,iorb1)=zpseu(1,iorb,iorb1)-eryd*zpseu(2,iorb,iorb1)
114      endif
115    enddo
116  enddo
117
118!
119! Forming the matrices A and B for generalized eigenvalue problem
120
121!
122!   1
123!
124  do n=1, 2*n2d
125    do ig=1, n2d
126      amat(ig, n)=fun1(ig, n)
127      amat(ig+n2d,n)=fund1(ig, n)
128      bmat(ig,n)=fun0(ig, n)
129      bmat(ig+n2d,n)=fund0(ig, n)
130     enddo
131  enddo
132
133!
134!   2
135!
136  do iorb=1, norb*npol
137    do ig=1, n2d
138      amat(ig, 2*n2d+iorb)=funl1(ig, iorb)
139      amat(n2d+ig, 2*n2d+iorb)=fundl1(ig, iorb)
140      bmat(ig, 2*n2d+iorb)=funl0(ig, iorb)
141      bmat(n2d+ig, 2*n2d+iorb)=fundl0(ig, iorb)
142    enddo
143  enddo
144!
145!  3
146!
147  do iorb=1, norb*npol
148    aorb=iorb
149    borb=iorb
150    if (iorb.le.npol*nocros) aorb=iorb+npol*(noins+nocros)
151    if (iorb.gt.npol*nocros) borb=iorb-npol*(noins+nocros)
152    do n=1, 2*n2d
153      do iorb1=1, norb*npol
154        if (noncolin) then
155           amat(2*n2d+iorb,n)=amat(2*n2d+iorb,n)+                  &
156                             zps_nc(aorb, iorb1)*intw1(iorb1,n)
157           if (borb.gt.0) bmat(2*n2d+iorb,n)=                      &
158              bmat(2*n2d+iorb,n)-zps_nc(borb,iorb1)*intw1(iorb1,n)
159        else
160           amat(2*n2d+iorb,n)=amat(2*n2d+iorb,n)+                  &
161                             zps(aorb,iorb1)*intw1(iorb1,n)
162           if (borb.gt.0) bmat(2*n2d+iorb,n)=                      &
163              bmat(2*n2d+iorb,n)-zps(borb,iorb1)*intw1(iorb1,n)
164        endif
165      enddo
166    enddo
167  enddo
168!
169!  4
170!
171  do iorb=1, nocros*npol
172    do iorb1=1, norb*npol
173      do iorb2=1, norb*npol
174        if (noncolin) then
175           bmat(2*n2d+iorb,2*n2d+iorb1)=bmat(2*n2d+iorb,2*n2d+iorb1) &
176                 -zps_nc(iorb,iorb2)*intw2(iorb2, iorb1)
177        else
178           bmat(2*n2d+iorb,2*n2d+iorb1)=bmat(2*n2d+iorb,2*n2d+iorb1) &
179                 -zps(iorb,iorb2)*intw2(iorb2, iorb1)
180        endif
181      enddo
182      bmat(2*n2d+iorb+npol*(noins+nocros),2*n2d+iorb1)=                  &
183                              bmat(2*n2d+iorb,2*n2d+iorb1)
184    enddo
185    bmat(2*n2d+iorb,2*n2d+iorb)=                                  &
186                         bmat(2*n2d+iorb,2*n2d+iorb)+(1.d0,0.d0)
187  enddo
188!
189!   5
190!
191  do iorb=1, norb*npol
192    aorb=iorb
193    if (iorb.le.npol*nocros) aorb=iorb+npol*(noins+nocros)
194    do iorb1=1, norb*npol
195      do iorb2=1, norb*npol
196        if (noncolin) then
197           amat(2*n2d+iorb,2*n2d+iorb1)=amat(2*n2d+iorb,2*n2d+iorb1)+  &
198                   zps_nc(aorb,iorb2)*intw2(iorb2, iorb1)
199        else
200           amat(2*n2d+iorb,2*n2d+iorb1)=amat(2*n2d+iorb,2*n2d+iorb1)+  &
201                   zps(aorb,iorb2)*intw2(iorb2, iorb1)
202        endif
203      enddo
204    enddo
205    if (aorb.eq.iorb) amat(2*n2d+iorb,2*n2d+iorb)=                  &
206                        amat(2*n2d+iorb,2*n2d+iorb)-(1.d0,0.d0)
207  enddo
208
209!
210! To reduce matrices and solve GEP A X = c B X; X = {a_n, a_\alpha}
211!
212
213  call compbs_2(npol*nocros, npol*norb, n2d, 2*(n2d+npol*nocros),   &
214                amat, bmat, vec, kval)
215
216!
217! To normalize (over XY plane) all the states
218!
219   call zgemm('n', 'n', n2d, 2*(n2d+npol*nocros), 2*n2d+npol*norb,   &
220               one, amat, 2*n2d+npol*norb, vec, 2*n2d+npol*norb,     &
221               zero, kfun, n2d)
222   do ig=1,n2d
223      do ik=1, 2*n2d+npol*norb
224         aux(ig,ik)=amat(n2d+ig,ik)
225      enddo
226   enddo
227   call zgemm('n', 'n', n2d, 2*(n2d+npol*nocros), 2*n2d+npol*norb,   &
228               one, aux, n2d, vec, 2*n2d+npol*norb, zero, kfund, n2d)
229
230  do ik=1, 2*(n2d+npol*nocros)
231    raux=ddot(2*n2d,kfun(1,ik),1,kfun(1,ik),1)*sarea
232    raux=1.d0/sqrt(raux)
233    call dscal(2*(2*n2d+npol*norb),raux,vec(1,ik),1)
234    call dscal(2*n2d,raux,kfun(1,ik),1)
235    call dscal(2*n2d,raux,kfund(1,ik),1)
236  enddo
237
238!
239! To find k-vector and the current of Bloch states
240!
241
242  call kbloch (2*(n2d+npol*nocros), kval)
243
244  call jbloch(2*(n2d+npol*nocros), n2d, norbf, norb, nocros,  &
245              kfun, kfund, vec, kval, intw1, intw2, nchan, npol)
246!
247! To save band structure result
248!
249  kfun=(0.d0,0.d0)
250  kfund=(0.d0,0.d0)
251  kint=(0.d0,0.d0)
252!
253! To account for the case of the right lead
254!
255  if (lleft.eq.0) then
256    do i=1, 2*n2d
257      do j=1, 2*n2d+npol*norb
258        amat(i,j)=bmat(i,j)
259      enddo
260    enddo
261    do i=2*n2d+1, 2*n2d+npol*nocros
262      do j=1, 2*n2d+npol*norb
263        amat(i,j)=-bmat(i+npol*(nocros+noins),j)
264      enddo
265    enddo
266  endif
267!
268! psi_k and psi'_k on the scattering region boundary
269!
270   call zgemm('n', 'n', n2d, 2*(n2d+npol*nocros), 2*n2d+npol*norb,&
271               one, amat, 2*n2d+npol*norb, vec, 2*n2d+npol*norb,  &
272               zero, kfun, n2d)
273   do ig=1,n2d
274      do ik=1, 2*n2d+npol*norb
275         aux(ig,ik)=amat(n2d+ig,ik)
276      enddo
277   enddo
278   call zgemm('n', 'n', n2d, 2*(n2d+npol*nocros), 2*n2d+npol*norb,&
279               one, aux, n2d, vec, 2*n2d+npol*norb, zero, kfund, n2d)
280
281!
282! kint(iorb, ik)=\sum_{iorb1} D_{iorb,iorb1}
283!            \int_{cell} W_iorb1^* psi_ik )
284! for the orbitals crossing the boundary
285!
286  do ik=1,  2*(n2d+npol*nocros)
287    do iorb=1, nocros*npol
288      do j=1, 2*n2d+npol*norb
289        kint(iorb,ik)=kint(iorb,ik)+amat(2*n2d+iorb,j)*vec(j,ik)
290      enddo
291    enddo
292  enddo
293!
294! a_iorb = kcoef(iorb,ik) = \sum_{iorb1} D_{iorb,iorb1}
295!           \int_{all space} W_iorb1^* psi_ik )
296! for the orbitals crossing the boundary
297!
298  do ik=1, 2*(n2d+npol*nocros)
299    do iorb=1, nocros*npol
300      if (lleft.eq.0) then
301        kcoef(iorb,ik)=vec(2*n2d+iorb,ik)
302      else
303        kcoef(iorb,ik)=vec(2*n2d+npol*(nocros+noins)+iorb,ik)
304      endif
305    enddo
306  enddo
307
308!
309! to set up B.S. for the right lead in the case of identical tips
310!
311  if(lleft.eq.1.and.ikind.eq.1) then
312    nchanr=nchan
313    call dcopy(2*(n2d+npol*nocros), kval, 1, kvalr, 1)
314    kfunr=(0.d0,0.d0)
315    kfundr=(0.d0,0.d0)
316    kintr=(0.d0,0.d0)
317
318    do i=1, 2*n2d
319      do j=1, 2*n2d+npol*norb
320        amat(i,j)=bmat(i,j)
321      enddo
322    enddo
323    do i=2*n2d+1, 2*n2d+npol*nocros
324      do j=1, 2*n2d+npol*norb
325        amat(i,j)=-bmat(i+npol*(nocros+noins),j)
326      enddo
327    enddo
328
329    do ik=1, 2*(n2d+npol*nocros)
330      do ig=1, n2d
331        do j=1,  2*n2d+npol*norb
332          kfunr(ig,ik)= kfunr(ig,ik)+amat(ig,j)*vec(j,ik)
333          kfundr(ig,ik)= kfundr(ig,ik)+amat(n2d+ig,j)*vec(j,ik)
334        enddo
335      enddo
336    enddo
337    do ik=1,  2*(n2d+npol*nocros)
338      do iorb=1, nocros*npol
339        do j=1, 2*n2d+npol*norb
340          kintr(iorb,ik)=kintr(iorb,ik)+amat(2*n2d+iorb,j)*vec(j,ik)
341        enddo
342      enddo
343    enddo
344    do ik=1, 2*(n2d+npol*nocros)
345      do iorb=1, nocros*npol
346        kcoefr(iorb,ik)=vec(2*n2d+iorb,ik)
347      enddo
348    enddo
349
350  endif
351
352!--
353!  integrals of Bloch states with boundary orbitals for left/right leads
354!
355 if (lorb) then
356   korb = 0.d0
357   do ik = 1, 2*(n2d+npol*nocros)
358
359     do iorb = 1, npol*nocros
360      iorb1 = iorb + npol*(nocros+noins)
361      do ig = 1, 2*n2d
362        korb(iorb,ik) = korb(iorb,ik)+       &
363                      intw1(iorb1,ig)*vec(ig,ik)
364      enddo
365      do ig = 1, norb*npol
366        korb(iorb,ik) = korb(iorb,ik)+       &
367               intw2(iorb1,ig)*vec(2*n2d+ig,ik)
368      enddo
369     enddo
370
371     do iorb = 1, npol*nocros
372      x1 = 0.d0
373      do ig = 1, 2*n2d
374        x1 = x1 + intw1(iorb,ig)*vec(ig,ik)
375      enddo
376      do ig = 1, norb*npol
377        x1 = x1 + intw2(iorb,ig)*vec(2*n2d+ig,ik)
378      enddo
379      korb(iorb,ik) = korb(iorb,ik) + x1* &
380                   exp(kval(ik)*(0.d0,1.d0)*tpi)
381     enddo
382   enddo
383
384   if (ikind.ne.2.or.lleft.ne.0) korbl(:,:) = korb(:,:)
385   if (ikind.ne.2.or.lleft.ne.1) then
386     do ik = 1, 2*(n2d+npol*nocros)
387       x1 = exp(-kval(ik)*(0.d0,1.d0)*tpi)
388       do iorb = 1, npol*nocros
389        korbr(iorb,ik) = korb(iorb,ik) * x1
390       enddo
391     enddo
392   endif
393
394 endif
395!--
396
397!--
398! Computes and writes the propagating Bloch states
399!
400  if (lorb.and.ikind.eq.0.and.nchan /= 0) then
401    allocate( veceig(nchan, nchan) )
402    deallocate( aux )
403    allocate( aux(4*n2d+npol*(norb+2*nocros),nchan) )
404
405!-- right moving states
406    veceig = 0.d0
407    aux = 0.d0
408    do ichan = 1, nchan
409      do ig = 1, 2*n2d + npol*norb
410        aux(ig, ichan) = vec(ig,ichan)
411      enddo
412    enddo
413    CALL scat_states_plot(ikk,ien,norb,nocros,nchan,aux,veceig,.true.)
414
415!-- left moving states
416    veceig = 0.d0
417    aux = 0.d0
418    do ichan = 1, nchan
419      do ig = 1, 2*n2d + npol*norb
420        aux(ig, ichan) = vec(ig,n2d+npol*nocros+ichan)
421      enddo
422    enddo
423    CALL scat_states_plot(ikk,ien,norb,nocros,nchan,aux,veceig,.false.)
424
425    deallocate( veceig )
426  endif
427!--
428
429  deallocate(amat)
430  deallocate(bmat)
431  deallocate(vec)
432  deallocate(aux)
433  IF (norb>0) THEN
434     if (noncolin) then
435        deallocate(zpseu_nc)
436        deallocate(zps_nc)
437     else
438        deallocate(zpseu)
439        deallocate(zps)
440     endif
441  ENDIF
442  if (lorb) deallocate(korb)
443  call stop_clock('compbs')
444
445  return
446end subroutine compbs
447