hsvac.F90 8.26 KB
Newer Older
Daniel Wortmann's avatar
Daniel Wortmann committed
1 2 3 4 5
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------
6
MODULE m_hsvac
7
  USE m_juDFT
8
CONTAINS
Daniel Wortmann's avatar
Daniel Wortmann committed
9 10 11 12
  !-----------------------------------------------------------
  !Calculate the vacuum contribution to the Hamiltonian and
  !Overlap matrix
  !-----------------------------------------------------------
13
  SUBROUTINE hsvac(&
Daniel Wortmann's avatar
Daniel Wortmann committed
14 15 16
       vacuum,stars,DIMENSION, atoms,mpi,jsp,input,v,evac,cell,&
       lapw,sym, noco,jij,hmat,smat)
 
17 18 19 20

    USE m_vacfun
    USE m_types
    IMPLICIT NONE
21 22 23 24 25 26 27 28 29 30
    TYPE(t_dimension),INTENT(IN)  :: DIMENSION
    TYPE(t_input),INTENT(IN)      :: input
    TYPE(t_vacuum),INTENT(IN)     :: vacuum
    TYPE(t_noco),INTENT(IN)       :: noco
    TYPE(t_jij),INTENT(IN)        :: jij
    TYPE(t_sym),INTENT(IN)        :: sym
    TYPE(t_stars),INTENT(IN)      :: stars
    TYPE(t_cell),INTENT(IN)       :: cell
    TYPE(t_atoms),INTENT(IN)      :: atoms
    TYPE(t_lapw),INTENT(IN)       :: lapw
Daniel Wortmann's avatar
Daniel Wortmann committed
31 32 33
    TYPE(t_mpi),INTENT(IN)        :: mpi
    TYPE(t_potden),INTENT(IN)     :: v
    CLASS(t_mat),INTENT(INOUT)     :: hmat(:,:),smat(:,:)
34 35
    !     ..
    !     .. Scalar Arguments ..
Daniel Wortmann's avatar
Daniel Wortmann committed
36
    INTEGER, INTENT (IN) :: jsp
37 38
    !     ..
    !     .. Array Arguments ..
39
    REAL,    INTENT (IN) :: evac(2,DIMENSION%jspd)
40 41 42
    !     ..
    !     .. Local Scalars ..
    COMPLEX hij,sij,apw_lo,c_1
Daniel Wortmann's avatar
Daniel Wortmann committed
43 44 45
    REAL d2,gz,sign,th,wronk
    INTEGER i,i2,ii,jj,ik,j,jk,k,jspin,ipot,ii0,i0
    INTEGER ivac,irec,imz,igvm2,igvm2i,s1,s2
46 47 48 49
    INTEGER jspin1,jspin2,jmax,jsp_start,jsp_end
    INTEGER i_start,nc,nc_0
    !     ..
    !     .. Local Arrays ..
Daniel Wortmann's avatar
Daniel Wortmann committed
50
    INTEGER:: nv2(DIMENSION%jspd)
51 52 53 54
    INTEGER kvac1(DIMENSION%nv2d,DIMENSION%jspd),kvac2(DIMENSION%nv2d,DIMENSION%jspd)
    INTEGER map2(DIMENSION%nvd,DIMENSION%jspd)
    COMPLEX tddv(DIMENSION%nv2d,DIMENSION%nv2d),tduv(DIMENSION%nv2d,DIMENSION%nv2d)
    COMPLEX tudv(DIMENSION%nv2d,DIMENSION%nv2d),tuuv(DIMENSION%nv2d,DIMENSION%nv2d)
55
    COMPLEX vxy_help(stars%ng2-1)
56 57 58 59
    COMPLEX a(DIMENSION%nvd,DIMENSION%jspd),b(DIMENSION%nvd,DIMENSION%jspd)
    REAL ddnv(DIMENSION%nv2d,DIMENSION%jspd),dudz(DIMENSION%nv2d,DIMENSION%jspd)
    REAL duz(DIMENSION%nv2d,DIMENSION%jspd), udz(DIMENSION%nv2d,DIMENSION%jspd)
    REAL uz(DIMENSION%nv2d,DIMENSION%jspd)
60 61 62 63
    ! l_J auxiliary potential array
    COMPLEX, ALLOCATABLE :: vxy1(:,:,:)
    !     ..

64
    d2 = SQRT(cell%omtil/cell%area)
65

66
    IF (jij%l_J) ALLOCATE (vxy1(vacuum%nmzxyd,stars%ng2-1,2))
67 68 69 70 71 72 73 74 75 76 77 78

    !--->    set up mapping function from 3d-->2d lapws

    DO jspin = 1,input%jspins
       nv2(jspin) = 0
       k_loop:DO  k = 1,lapw%nv(jspin)
          DO  j = 1,nv2(jspin)
             IF (lapw%k1(k,jspin).EQ.kvac1(j,jspin)&
                  .AND. lapw%k2(k,jspin).EQ.kvac2(j,jspin)) THEN
                map2(k,jspin) = j
                CYCLE k_loop
             END IF
79
          ENDDO
80
          nv2(jspin) = nv2(jspin) + 1
81
          IF (nv2(jspin)>DIMENSION%nv2d)  CALL juDFT_error("hsvac:dimension%nv2d",calledby ="hsvac")
82 83 84
          kvac1(nv2(jspin),jspin) = lapw%k1(k,jspin)
          kvac2(nv2(jspin),jspin) = lapw%k2(k,jspin)
          map2(k,jspin) = nv2(jspin)
85
       ENDDO k_loop
86 87 88
    ENDDO
    !--->    loop over the two vacuua (1: upper; 2: lower)
    DO ivac = 1,2
Daniel Wortmann's avatar
Daniel Wortmann committed
89 90 91 92 93 94 95 96 97
       sign = 3. - 2.*ivac !+/- 1
       DO jspin1=MERGE(1,jsp,noco%l_noco),MERGE(2,jsp,noco%l_noco) !loop over global spin
          s1=MIN(SIZE(hmat,1),jspin1) !in colinear case s1=1
          DO jspin2=MERGE(1,jsp,noco%l_noco),MERGE(2,jsp,noco%l_noco) !loop over global spin
             s2=MIN(SIZE(hmat,1),jspin2) !in colinear case s2=1
             ipot=3
             IF (jspin1==jspin2) ipot=jspin1
          !--->       get the wavefunctions and set up the tuuv, etc matrices          
             jspin=jsp
98
          CALL vacfun(&
99
               vacuum,DIMENSION,stars,&
100
               jsp,input,noco,ipot,&
Daniel Wortmann's avatar
Daniel Wortmann committed
101
               sym, cell,ivac,evac(1,1),lapw%bkpt,v%vacxy(:,:,ivac,ipot),v%vacz(:,:,:),kvac1,kvac2,nv2,&
102 103 104 105 106 107 108 109 110 111
               tuuv,tddv,tudv,tduv,uz,duz,udz,dudz,ddnv,wronk)
          !
          !--->       generate a and b coeffficients
          !
          IF (noco%l_noco) THEN
             DO jspin = 1,input%jspins
                DO k = 1,lapw%nv(jspin)
                   gz = sign*cell%bmat(3,3)*lapw%k3(k,jspin)
                   i2 = map2(k,jspin)
                   th = gz*cell%z1
Daniel Wortmann's avatar
Daniel Wortmann committed
112
                   c_1 = CMPLX( COS(th), SIN(th) )/ (d2*wronk)
113 114
                   a(k,jspin) = - c_1 * CMPLX(dudz(i2,jspin), gz*udz(i2,jspin) )
                   b(k,jspin) =   c_1 * CMPLX(duz(i2,jspin), gz* uz(i2,jspin) )
115 116 117 118
                ENDDO
             ENDDO
          ELSE
             DO k = 1,lapw%nv(jsp)
Daniel Wortmann's avatar
Daniel Wortmann committed
119
                gz = sign*cell%bmat(3,3)*lapw%gvec(3,k,jsp)
120 121
                i2 = map2(k,jsp)
                th = gz*cell%z1
Daniel Wortmann's avatar
Daniel Wortmann committed
122 123 124
                c_1 = CMPLX( COS(th), SIN(th) )/ (d2*wronk)
                a(k,jsp) = - c_1 * CMPLX(dudz(i2,jsp), gz*udz(i2,jsp) )
                b(k,jsp) =   c_1 * CMPLX(duz(i2,jsp), gz* uz(i2,jsp) )
125 126 127
             ENDDO
          ENDIF
          !--->       update hamiltonian and overlap matrices
Daniel Wortmann's avatar
Daniel Wortmann committed
128 129 130
          IF (jspin1==jspin2) THEN
             DO  i = mpi%n_rank+1,lapw%nv(jspin1),mpi%n_size
                i0=(i-1)/mpi%n_size+1 !local column index
131
                ik = map2(i,jspin)
Daniel Wortmann's avatar
Daniel Wortmann committed
132
                DO j = 1,i - 1 !TODO check noco case
133 134
                   !--->             overlap: only  (g-g') parallel=0       '
                   IF (map2(j,jspin).EQ.ik) THEN
135 136
                      sij = CONJG(a(i,jspin))*a(j,jspin) + &
                           CONJG(b(i,jspin))*b(j,jspin)*ddnv(ik,jspin1)
137
                      !+APW_LO
138 139 140 141 142 143
                      IF (input%l_useapw) THEN
                         apw_lo = CONJG(a(i,jspin)*  uz(ik,jspin1) + b(i,jspin)* udz(ik,jspin1) ) &
                              * (a(j,jspin)* duz(ik,jspin1) + b(j,jspin)*dudz(ik,jspin1) )&
                              +      (a(j,jspin)*  uz(ik,jspin1) + b(j,jspin)* udz(ik,jspin1) ) &
                              * CONJG(a(i,jspin)* duz(ik,jspin1) + b(i,jspin)*dudz(ik,jspin1) )
                         !            IF (i.lt.10) write (3,'(2i4,2f20.10)') i,j,apw_lo
Daniel Wortmann's avatar
Daniel Wortmann committed
144 145
                         IF (hmat(1,1)%l_real) THEN
                            hmat(s1,s2)%data_r(j,i0) = hmat(s1,s2)%data_r(j,i0) + 0.25 * REAL(apw_lo) 
146
                         ELSE 
Daniel Wortmann's avatar
Daniel Wortmann committed
147
                            hmat(s1,s2)%data_c(j,i0) = hmat(s1,s2)%data_c(j,i0) + 0.25 * apw_lo
148
                         ENDIF
149
                      ENDIF
Daniel Wortmann's avatar
Daniel Wortmann committed
150 151 152
                      !Overlapp Matrix
                      IF (hmat(1,1)%l_real) THEN
                         smat(s1,s2)%data_r(j,i0) = smat(s1,s2)%data_r(j,i0) + REAL(sij)
153
                      ELSE 
Daniel Wortmann's avatar
Daniel Wortmann committed
154
                         smat(s1,s2)%data_c(j,i0) = smat(s1,s2)%data_c(j,i0) + sij
155
                      ENDIF
156 157
                   END IF
                ENDDO
Daniel Wortmann's avatar
Daniel Wortmann committed
158
                !Diagonal term of Overlapp matrix, Hamiltonian later
159
                sij = CONJG(a(i,jspin))*a(i,jspin) + CONJG(b(i,jspin))*b(i,jspin)*ddnv(ik,jspin1)
Daniel Wortmann's avatar
Daniel Wortmann committed
160 161
                IF (hmat(1,1)%l_real) THEN
                   smat(s1,s2)%data_r(j,i0) = smat(s1,s2)%data_r(j,i0) + REAL(sij)
162
                ELSE
Daniel Wortmann's avatar
Daniel Wortmann committed
163
                   smat(s1,s2)%data_c(j,i0) = smat(s1,s2)%data_c(j,i0) + sij
164
                ENDIF
165
             ENDDO
166 167 168
          ENDIF

          !--->    hamiltonian update
Daniel Wortmann's avatar
Daniel Wortmann committed
169 170
          DO  i = mpi%n_rank+1,lapw%nv(jspin2),mpi%n_size
             i0=(i-1)/mpi%n_size+1 !local column index
171
             ik = map2(i,jspin1)
Daniel Wortmann's avatar
Daniel Wortmann committed
172
             DO j = 1,i
173
                jk = map2(j,jspin2)
174 175
                hij = CONJG(a(i,jspin1))* (tuuv(ik,jk)*a(j,jspin2) +tudv(ik,jk)*b(j,jspin2))&
                     + CONJG(b(i,jspin1))* (tddv(ik,jk)*b(j,jspin2) +tduv(ik,jk)*a(j,jspin2))
Daniel Wortmann's avatar
Daniel Wortmann committed
176 177
                IF (hmat(1,1)%l_real) THEN
                   hmat(s1,s2)%data_r(j,i0) = hmat(s1,s2)%data_r(j,i0) + REAL(hij)
178
                ELSE
Daniel Wortmann's avatar
Daniel Wortmann committed
179
                   hmat(s1,s2)%data_c(j,i0) = hmat(s1,s2)%data_c(j,i0) + hij
180
                ENDIF
181 182 183 184 185 186 187 188
             ENDDO
          ENDDO

          !--->    end of loop over different parts of the potential matrix
       ENDDO

       !---> end of loop over vacua
    ENDDO
Daniel Wortmann's avatar
Daniel Wortmann committed
189
 ENDDO
190 191 192

  END SUBROUTINE hsvac
END MODULE m_hsvac