LCOV - code coverage report
Current view: top level - eigen - od_hsvac.F90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 0 163 0.0 %
Date: 2019-09-08 04:53:50 Functions: 0 1 0.0 %

          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_od_hsvac
       8             :   USE m_juDFT
       9             : CONTAINS
      10           0 :   SUBROUTINE od_hsvac(&
      11           0 :        vacuum,stars,DIMENSION, oneD,atoms, jsp,input,vxy,vz,evac,cell,&
      12           0 :        bkpt,lapw, MM,vM,m_cyl,n2d_1, n_size,n_rank,sym,noco,nv2,l_real,hamOvlp)
      13             : 
      14             :     !     subroutine for calculating the hamiltonian and overlap matrices in
      15             :     !     the vacuum in the case of 1-dimensional calculations
      16             :     !     Y. Mokrousov June 2002             
      17             : 
      18             :     USE m_cylbes
      19             :     USE m_dcylbs
      20             :     USE m_od_vacfun
      21             :     USE m_types
      22             :     IMPLICIT NONE
      23             :     TYPE(t_dimension),INTENT(IN)  :: DIMENSION
      24             :     TYPE(t_oneD),INTENT(IN)       :: oneD
      25             :     TYPE(t_input),INTENT(IN)      :: input
      26             :     TYPE(t_vacuum),INTENT(IN)     :: vacuum
      27             :     TYPE(t_noco),INTENT(IN)       :: noco
      28             :     TYPE(t_sym),INTENT(IN)        :: sym
      29             :     TYPE(t_stars),INTENT(IN)      :: stars
      30             :     TYPE(t_cell),INTENT(IN)       :: cell
      31             :     TYPE(t_atoms),INTENT(IN)      :: atoms
      32             :     TYPE(t_lapw),INTENT(IN)       :: lapw
      33             :     TYPE(t_hamOvlp),INTENT(INOUT) :: hamOvlp
      34             :     !     ..
      35             :     !     .. Scalar Arguments ..
      36             :     INTEGER, INTENT (IN) :: vM
      37             :     INTEGER, INTENT (IN) :: MM 
      38             :     INTEGER, INTENT (IN) :: jsp ,n_size,n_rank,n2d_1 
      39             :     INTEGER, INTENT (IN) :: m_cyl
      40             :     !     ..
      41             :     !     .. Array Arguments ..
      42             :     COMPLEX, INTENT (INOUT) :: vxy(vacuum%nmzxyd,n2d_1-1,2)
      43             :     INTEGER, INTENT (OUT):: nv2(input%jspins)
      44             :     REAL,    INTENT (INOUT) :: vz(vacuum%nmzd,2,4)
      45             :     REAL,    INTENT (IN) :: evac(2,input%jspins)
      46             :     REAL,    INTENT (IN) :: bkpt(3)
      47             : 
      48             :     LOGICAL, INTENT(IN)  :: l_real
      49             : 
      50             :     !     ..
      51             :     !     .. Local Scalars ..
      52             :     COMPLEX hij,sij,apw_lo,exp1,exp2,exp3,am,bm,ic
      53             :     REAL    d2,wronk,gr,gphi,qq,x,y
      54             :     INTEGER i,i2,ii,ik,j,jk,k,jspin,ipot,npot,ii0 ,l,i3,imz,m
      55             :     INTEGER jspin1,jspin2,jmax,irec2,irec3,ivac,ind1,gi
      56             :     INTEGER i_start,nc,nc_0,rotax,chiral,zi,m1,z,indm,indl
      57             :     !     ..
      58             :     !     .. Local Arrays ..
      59             : 
      60           0 :     INTEGER, ALLOCATABLE :: nvp(:,:),ind(:,:,:)
      61           0 :     INTEGER, ALLOCATABLE :: kvac3(:,:),map1(:,:)
      62           0 :     COMPLEX, ALLOCATABLE :: tddv(:,:,:,:)
      63           0 :     COMPLEX, ALLOCATABLE :: tduv(:,:,:,:)
      64           0 :     COMPLEX, ALLOCATABLE :: tudv(:,:,:,:)
      65           0 :     COMPLEX, ALLOCATABLE :: tuuv(:,:,:,:)
      66           0 :     COMPLEX, ALLOCATABLE ::  a(:,:,:),b(:,:,:)
      67           0 :     COMPLEX, ALLOCATABLE :: ai(:,:,:),bi(:,:,:)
      68           0 :     REAL, ALLOCATABLE :: bess(:),dbss(:),bess1(:)
      69           0 :     REAL, ALLOCATABLE :: ddnv(:,:,:),dudz(:,:,:)
      70           0 :     REAL, ALLOCATABLE :: duz(:,:,:)
      71           0 :     REAL, ALLOCATABLE :: udz(:,:,:),uz(:,:,:)
      72             :     !     ..
      73           0 :     ic  = CMPLX(0.,1.)
      74           0 :     d2 = SQRT(cell%omtil/cell%area)
      75             : 
      76             :  
      77             :     ALLOCATE (&
      78             :          ai(-vM:vM,DIMENSION%nv2d,DIMENSION%nvd),bi(-vM:vM,DIMENSION%nv2d,DIMENSION%nvd),&
      79             :          nvp(DIMENSION%nv2d,input%jspins),ind(stars%ng2,DIMENSION%nv2d,input%jspins),&
      80             :          kvac3(DIMENSION%nv2d,input%jspins),map1(DIMENSION%nvd,input%jspins),&
      81             :          tddv(-vM:vM,-vM:vM,DIMENSION%nv2d,DIMENSION%nv2d),&
      82             :          tduv(-vM:vM,-vM:vM,DIMENSION%nv2d,DIMENSION%nv2d),&
      83             :          tudv(-vM:vM,-vM:vM,DIMENSION%nv2d,DIMENSION%nv2d),&
      84             :          tuuv(-vM:vM,-vM:vM,DIMENSION%nv2d,DIMENSION%nv2d),&
      85             :          a(-vM:vM,DIMENSION%nvd,input%jspins),b(-vM:vM,DIMENSION%nvd,input%jspins),&
      86             :          bess(-vM:vM),dbss(-vM:vM),bess1(-vM:vM),&
      87             :          ddnv(-vM:vM,DIMENSION%nv2d,input%jspins),dudz(-vM:vM,DIMENSION%nv2d,input%jspins),&
      88             :          duz(-vM:vM,DIMENSION%nv2d,input%jspins),&
      89           0 :          udz(-vM:vM,DIMENSION%nv2d,input%jspins),uz(-vM:vM,DIMENSION%nv2d,input%jspins) )
      90             : 
      91             :     !--->     set up mapping function from 3d-->1d lapws
      92             :     !--->            creating arrays ind and nvp
      93             : 
      94           0 :     DO jspin = 1,input%jspins
      95             : 
      96           0 :        nv2(jspin) = 0
      97           0 :        k_loop:DO  k = 1,lapw%nv(jspin)
      98           0 :           DO  j = 1,nv2(jspin)
      99           0 :              IF (lapw%k3(k,jspin).EQ.kvac3(j,jspin)) THEN
     100           0 :                 map1(k,jspin) = j
     101           0 :                 CYCLE k_loop
     102             :              END IF
     103             :           ENDDO
     104           0 :           nv2(jspin) = nv2(jspin) + 1
     105           0 :           IF (nv2(jspin)>DIMENSION%nv2d)  CALL juDFT_error("dimension%nv2d",calledby ="od_hsvac")
     106           0 :           kvac3(nv2(jspin),jspin) = lapw%k3(k,jspin)
     107           0 :           map1(k,jspin) = nv2(jspin)
     108             :        END DO k_loop
     109             : 
     110           0 :        DO ik = 1,DIMENSION%nv2d
     111           0 :           nvp(ik,jspin) = 0
     112           0 :           DO i = 1,stars%ng2
     113           0 :              ind(i,ik,jspin) = 0
     114             :           END DO
     115             :        END DO
     116             : 
     117           0 :        DO k = 1,lapw%nv(jspin)
     118           0 :           ik = map1(k,jspin)
     119           0 :           nvp(ik,jspin) = nvp(ik,jspin) + 1
     120           0 :           ind(nvp(ik,jspin),ik,jspin) = k
     121             :        END DO
     122             : 
     123             :     ENDDO
     124             : 
     125           0 :     npot = 1      
     126           0 :     ivac = 1
     127             : 
     128           0 :     IF (noco%l_noco) THEN
     129             :        !--->         load the non-warping part of the potential
     130           0 :        READ (25)((vz(imz,ivac,ipot),imz=1,vacuum%nmzd),ipot=1,4)
     131           0 :        npot = 3
     132             :     ENDIF
     133             : 
     134           0 :     DO ipot = 1,npot
     135             : 
     136           0 :        IF (noco%l_noco) THEN
     137           0 :           READ (25)((vxy(imz,k,ivac), imz=1,vacuum%nmzxy),k=1,n2d_1-1)
     138             :           !--->  l_J we want to average the diagonal elements of the pot. matrix
     139             :        ENDIF ! loco
     140             : 
     141             :        !     get the wavefunctions and set up the tuuv, etc matrices
     142             : 
     143             :        CALL od_vacfun(&
     144             :             m_cyl,cell,vacuum,DIMENSION,stars,&
     145             :             jsp,input,noco,ipot,oneD,n2d_1, ivac,evac(1,1),bkpt,MM,vM,&
     146           0 :             vxy(1,1,ivac),vz,kvac3,nv2, tuuv,tddv,tudv,tduv,uz,duz,udz,dudz,ddnv)
     147             : 
     148           0 :        IF (noco%l_noco) THEN
     149             : 
     150           0 :           DO jspin = 1,input%jspins
     151             : 
     152           0 :              DO k = 1,lapw%nv(jspin)
     153           0 :                 irec3 = stars%ig(lapw%k1(k,jspin),lapw%k2(k,jspin),lapw%k3(k,jspin))
     154           0 :                 IF (irec3.NE.0) THEN
     155           0 :                    irec2 = stars%ig2(irec3)
     156           0 :                    gr = stars%sk2(irec2)
     157           0 :                    gphi = stars%phi2(irec2)
     158           0 :                    i2 = map1(k,jspin)
     159           0 :                    qq = gr*cell%z1
     160           0 :                    CALL cylbes(vM,qq,bess)
     161           0 :                    CALL dcylbs(vM,qq,bess,dbss)
     162           0 :                    DO m = -vM,vM
     163           0 :                       wronk = uz(m,i2,jspin)*dudz(m,i2,jspin) - udz(m,i2,jspin)*duz(m,i2,jspin)
     164             :                       a(m,k,jspin)=EXP(-CMPLX(0.0,m*gphi))*(ic**m)*&
     165           0 :                            CMPLX(dudz(m,i2,jspin)*bess(m)- udz(m,i2,jspin)*gr*dbss(m),0.0) /(d2*wronk)
     166             :                       b(m,k,jspin)=EXP(-CMPLX(0.0,m*gphi))*(ic**m)*&
     167           0 :                            CMPLX(-duz(m,i2,jspin)*bess(m)+ uz(m,i2,jspin)*gr*dbss(m),0.0) /(d2*wronk)
     168             :                    END DO
     169             :                 END IF
     170             :              ENDDO
     171             : 
     172             :           ENDDO  ! jspin
     173             : 
     174             :        ELSE 
     175             : 
     176           0 :           DO k = 1,lapw%nv(jsp)
     177           0 :              irec3 = stars%ig(lapw%k1(k,jsp),lapw%k2(k,jsp),lapw%k3(k,jsp))
     178           0 :              IF (irec3.NE.0) THEN
     179           0 :                 irec2 = stars%ig2(irec3)
     180           0 :                 gr = stars%sk2(irec2)
     181           0 :                 gphi = stars%phi2(irec2)
     182           0 :                 i2 = map1(k,jsp)
     183           0 :                 qq = gr*cell%z1
     184           0 :                 CALL cylbes(vM,qq,bess) 
     185           0 :                 CALL dcylbs(vM,qq,bess,dbss)
     186           0 :                 DO m = -vM,vM
     187           0 :                    wronk = uz(m,i2,jsp)*dudz(m,i2,jsp) - udz(m,i2,jsp)*duz(m,i2,jsp) 
     188             :                    a(m,k,1)=EXP(-CMPLX(0.0,m*gphi))*(ic**m)*&
     189           0 :                         CMPLX(dudz(m,i2,jsp)*bess(m)- udz(m,i2,jsp)*gr*dbss(m),0.0) /(d2*wronk)
     190             : 
     191             :                    b(m,k,1)=EXP(-CMPLX(0.0,m*gphi))*(ic**m)*&
     192           0 :                         CMPLX(-duz(m,i2,jsp)*bess(m)+ uz(m,i2,jsp)*gr*dbss(m),0.0) /(d2*wronk)
     193             : 
     194             :                 END DO
     195             :              END IF
     196             :           ENDDO
     197             : 
     198             :        ENDIF ! loco
     199             :        !     update hamiltonian and overlap matrices
     200             : 
     201           0 :        IF (ipot.EQ.1 .OR. ipot.EQ.2) THEN
     202           0 :           jspin = ipot
     203             :           !+gb||
     204           0 :           IF (ipot.EQ.1) THEN
     205           0 :              nc = 0
     206           0 :              i_start = n_rank
     207             :           ELSE
     208           0 :              nc = nc + atoms%nlotot
     209           0 :              nc_0 = nc
     210           0 :              i_start = MOD(MOD(n_rank - (lapw%nv(1)+atoms%nlotot),n_size) + n_size,n_size)
     211             :           ENDIF
     212             : 
     213           0 :           DO  i = i_start+1,lapw%nv(jspin),n_size
     214           0 :              ik = map1(i,jspin)
     215           0 :              nc = nc + 1
     216           0 :              IF (ipot.EQ.1) THEN
     217           0 :                 jspin = 1
     218           0 :                 ii0 = nc*(nc-1)/2*n_size- (nc-1)*(n_size-n_rank-1)
     219           0 :              ELSEIF (ipot.EQ.2) THEN
     220           0 :                 jspin = 2
     221             :                 ii0=nc*(nc-1)/2*n_size-(nc-1)*(n_size-n_rank-1)+&
     222           0 :                      lapw%nv(1)+atoms%nlotot
     223             :              ENDIF
     224           0 :              jspin1 = jsp
     225           0 :              IF (noco%l_noco) jspin1 = jspin
     226           0 :              DO j = 1,i - 1
     227           0 :                 ii = ii0 + j
     228             :                 !     overlap: only  (g-g') parallel=0        
     229           0 :                 IF (map1(j,jspin).EQ.ik) THEN
     230           0 :                    sij = (0.0,0.0)
     231           0 :                    DO m = -vM,vM
     232             :                       sij = sij + CONJG(a(m,i,jspin))*a(m,j,jspin) &
     233           0 :                            +CONJG(b(m,i,jspin))*b(m,j,jspin) *ddnv(m,ik,jspin1)
     234             :                    END DO
     235           0 :                    IF (l_real) THEN 
     236           0 :                       hamOvlp%b_r(ii) = hamOvlp%b_r(ii) + REAL(sij)
     237             :                    ELSE
     238           0 :                       hamOvlp%b_c(ii) = hamOvlp%b_c(ii) + sij
     239             :                    ENDIF
     240             :                 END IF
     241             :              ENDDO
     242           0 :              ii = ii0 + i
     243           0 :              sij = (0.0,0.0)
     244           0 :              DO m = -vM,vM
     245             :                 sij = sij + CONJG(a(m,i,jspin))*a(m,i,jspin)+ &
     246           0 :                      CONJG(b(m,i,jspin))*b(m,i,jspin)*ddnv(m,ik,jspin1)
     247             :              END DO
     248             : 
     249           0 :              IF (l_real) THEN
     250           0 :                 hamOvlp%b_r(ii) = hamOvlp%b_r(ii) + REAL(sij)
     251             :              ELSE 
     252           0 :                 hamOvlp%b_c(ii) = hamOvlp%b_c(ii) + sij
     253             :              ENDIF
     254             :           ENDDO
     255             :        ENDIF ! ipot.eq.1.or.2
     256             :        !   hamiltonian update 
     257             :        !   for the noncylindr. contributions we use the cutoff of m_cyl        
     258           0 :        IF (ipot.EQ.1) THEN
     259           0 :           jspin1 = 1
     260           0 :           jspin2 = 1
     261           0 :           nc = 0
     262           0 :           i_start = n_rank
     263           0 :        ELSEIF (ipot.EQ.2) THEN
     264           0 :           jspin1 = 2
     265           0 :           jspin2 = 2
     266           0 :           nc = nc_0
     267           0 :           i_start = MOD(MOD(n_rank - (lapw%nv(1)+atoms%nlotot),n_size) + n_size,n_size)
     268           0 :        ELSEIF (ipot.EQ.3) THEN
     269           0 :           jspin1 = 2
     270           0 :           jspin2 = 1
     271           0 :           nc = nc_0
     272           0 :           i_start = MOD(MOD(n_rank - (lapw%nv(1)+atoms%nlotot),n_size) + n_size,n_size)
     273             :        ENDIF
     274             : 
     275           0 :        ai(:,:,:) = CMPLX(0.,0.)
     276           0 :        bi(:,:,:) = CMPLX(0.,0.)
     277             : 
     278           0 :        DO ik = 1,nv2(jspin1)
     279           0 :           DO jk = 1,nv2(jspin2)
     280           0 :              i3 = kvac3(ik,jspin1) - kvac3(jk,jspin2) 
     281           0 :              DO l = -vM,vM
     282           0 :                 DO m = -vM,vM
     283           0 :                    IF (l.EQ.m .OR. (iabs(m).LE.m_cyl .AND. iabs(l).LE.m_cyl)) THEN
     284           0 :                       ind1 = oneD%ig1(i3,m-l)
     285           0 :                       IF (ind1.NE.0) THEN
     286           0 :                          DO gi = 1,nvp(ik,jspin1)
     287           0 :                             i = ind(gi,ik,jspin1)
     288             :                             ai(l,jk,i) = ai(l,jk,i) +&
     289           0 :                                  CONJG(a(m,i,jspin1))*tuuv(m,l,ik,jk) + CONJG(b(m,i,jspin1))*tduv(m,l,ik,jk)
     290             :                             bi(l,jk,i) = bi(l,jk,i) +&
     291           0 :                                  CONJG(a(m,i,jspin1))*tudv(m,l,ik,jk) + CONJG(b(m,i,jspin1))*tddv(m,l,ik,jk)
     292             : 
     293             :                          END DO
     294             :                       END IF
     295             :                    END IF   ! noncyl. contributions
     296             :                 END DO
     297             :              END DO
     298             :           END DO
     299             :        END DO
     300             : 
     301           0 :        DO i = i_start+1, lapw%nv(jspin1), n_size
     302           0 :           ik = map1(i,jspin1)
     303           0 :           nc = nc + 1
     304           0 :           IF (ipot.EQ.1) THEN
     305           0 :              ii0 = nc*(nc-1)/2*n_size- (nc-1)*(n_size-n_rank-1)
     306           0 :              jmax = i
     307           0 :           ELSEIF (ipot.EQ.2) THEN
     308           0 :              ii0 = nc*(nc-1)/2*n_size- (nc-1)*(n_size-n_rank-1) + lapw%nv(1)+atoms%nlotot
     309           0 :              jmax = i
     310           0 :           ELSEIF (ipot.EQ.3) THEN
     311           0 :              ii0 = nc*(nc-1)/2*n_size- (nc-1)*(n_size-n_rank-1)
     312           0 :              jmax = lapw%nv(jspin2)
     313             :           ENDIF
     314           0 :           DO j = 1,jmax
     315           0 :              ii = ii0 + j
     316           0 :              jk = map1(j,jspin2)
     317           0 :              hij = CMPLX(0.,0.)
     318           0 :              DO l = -vM,vM
     319           0 :                 hij = hij + ai(l,jk,i)*a(l,j,jspin2) + bi(l,jk,i)*b(l,j,jspin2)
     320             :              END DO
     321           0 :              IF (l_real) THEN
     322           0 :                 hamOvlp%a_r(ii) = hamOvlp%a_r(ii) + REAL(hij)
     323             :              ELSE 
     324           0 :                 hamOvlp%a_c(ii) = hamOvlp%a_c(ii) + hij
     325             :              ENDIF
     326             :           END DO
     327             :        END DO
     328             : 
     329             :     ENDDO !ipot
     330             : 
     331           0 :     DEALLOCATE (ai,bi,nvp,ind,kvac3,map1, tddv,tduv,tudv,tuuv,a,b,bess,dbss,bess1, ddnv,dudz,duz,udz,uz )
     332             : 
     333           0 :     RETURN
     334             :   END SUBROUTINE od_hsvac
     335             : END MODULE m_od_hsvac

Generated by: LCOV version 1.13