hsvac.F90 7.41 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(&
14
       vacuum,stars,DIMENSION, mpi,jsp,input,v,evac,cell,&
15
       lapw,sym, noco,hmat,smat)
Daniel Wortmann's avatar
Daniel Wortmann committed
16
 
17 18 19 20

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

60

61
    d2 = SQRT(cell%omtil/cell%area)
62 63 64 65 66 67 68

    !--->    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)
69
             IF (all(lapw%gvec(1:2,k,jspin)==kvac(1:2,j,jspin))) THEN
70 71 72
                map2(k,jspin) = j
                CYCLE k_loop
             END IF
73
          ENDDO
74
          nv2(jspin) = nv2(jspin) + 1
75
          IF (nv2(jspin)>DIMENSION%nv2d)  CALL juDFT_error("hsvac:dimension%nv2d",calledby ="hsvac")
76
          kvac(1:2,nv2(jspin),jspin) = lapw%gvec(1:2,k,jspin)
77
          map2(k,jspin) = nv2(jspin)
78
       ENDDO k_loop
79 80 81
    ENDDO
    !--->    loop over the two vacuua (1: upper; 2: lower)
    DO ivac = 1,2
Daniel Wortmann's avatar
Daniel Wortmann committed
82 83 84 85 86 87
       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
          !--->       get the wavefunctions and set up the tuuv, etc matrices          
88 89 90 91 92
             CALL vacfun(&
                  vacuum,stars,&
                  input,noco,jspin1,jspin2,&
                  sym, cell,ivac,evac,lapw%bkpt,v%vacxy,v%vacz,kvac,nv2,&
                  tuuv,tddv,tudv,tduv,uz,duz,udz,dudz,ddnv,wronk)
93 94 95
          !
          !--->       generate a and b coeffficients
          !
96
             DO jspin = MIN(jspin1,jspin2),MAX(jspin1,jspin2)
97 98 99 100
                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
101
                   c_1 = CMPLX( COS(th), SIN(th) )/ (d2*wronk)
102 103
                   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) )
104 105 106
                ENDDO
             ENDDO
          !--->       update hamiltonian and overlap matrices
Daniel Wortmann's avatar
Daniel Wortmann committed
107
          IF (jspin1==jspin2) THEN
108
             DO  i = mpi%n_rank+1,lapw%nv(jspin2),mpi%n_size
Daniel Wortmann's avatar
Daniel Wortmann committed
109
                i0=(i-1)/mpi%n_size+1 !local column index
110
                ik = map2(i,jspin2)
Daniel Wortmann's avatar
Daniel Wortmann committed
111
                DO j = 1,i - 1 !TODO check noco case
112
                   !--->             overlap: only  (g-g') parallel=0       '
113
                   IF (map2(j,jspin1).EQ.ik) THEN
114 115
                      sij = CONJG(a(i,jspin2))*a(j,jspin2) + &
                           CONJG(b(i,jspin2))*b(j,jspin2)*ddnv(ik,jspin2)
116
                      !+APW_LO
117
                      IF (input%l_useapw) THEN
118 119 120 121
                         apw_lo = CONJG(a(i,jspin1)*  uz(ik,jspin1) + b(i,jspin1)* udz(ik,jspin1) ) &
                              * (a(j,jspin1)* duz(ik,jspin1) + b(j,jspin1)*dudz(ik,jspin1) )&
                              +      (a(j,jspin1)*  uz(ik,jspin1) + b(j,jspin1)* udz(ik,jspin1) ) &
                              * CONJG(a(i,jspin1)* duz(ik,jspin1) + b(i,jspin1)*dudz(ik,jspin1) )
122
                         !            IF (i.lt.10) write (3,'(2i4,2f20.10)') i,j,apw_lo
Daniel Wortmann's avatar
Daniel Wortmann committed
123 124
                         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) 
125
                         ELSE 
Daniel Wortmann's avatar
Daniel Wortmann committed
126
                            hmat(s1,s2)%data_c(j,i0) = hmat(s1,s2)%data_c(j,i0) + 0.25 * apw_lo
127
                         ENDIF
128
                      ENDIF
Daniel Wortmann's avatar
Daniel Wortmann committed
129 130 131
                      !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)
132
                      ELSE 
Daniel Wortmann's avatar
Daniel Wortmann committed
133
                         smat(s1,s2)%data_c(j,i0) = smat(s1,s2)%data_c(j,i0) + sij
134
                      ENDIF
135 136
                   END IF
                ENDDO
Daniel Wortmann's avatar
Daniel Wortmann committed
137
                !Diagonal term of Overlapp matrix, Hamiltonian later
138
                sij = CONJG(a(i,jspin2))*a(i,jspin2) + CONJG(b(i,jspin2))*b(i,jspin2)*ddnv(ik,jspin2)
Daniel Wortmann's avatar
Daniel Wortmann committed
139
                IF (hmat(1,1)%l_real) THEN
140
                   smat(s2,s1)%data_r(j,i0) = smat(s1,s2)%data_r(j,i0) + REAL(sij)
141
                ELSE
142
                   smat(s2,s1)%data_c(j,i0) = smat(s1,s2)%data_c(j,i0) + sij
143
                ENDIF
144
             ENDDO
145 146 147
          ENDIF

          !--->    hamiltonian update
148
          DO  i = mpi%n_rank+1,lapw%nv(jspin1),mpi%n_size
Daniel Wortmann's avatar
Daniel Wortmann committed
149
             i0=(i-1)/mpi%n_size+1 !local column index
150 151 152 153 154
             ik = map2(i,jspin1)
             DO j = 1,MERGE(i,lapw%nv(jspin2),jspin1==jspin2)
                jk = map2(j,jspin2)
                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
155
                IF (hmat(1,1)%l_real) THEN
156
                   hmat(s2,s1)%data_r(j,i0) = hmat(s2,s1)%data_r(j,i0) + REAL(hij)
157
                ELSE
158
                   hmat(s2,s1)%data_c(j,i0) = hmat(s2,s1)%data_c(j,i0) + hij
159
                ENDIF
160 161 162 163 164 165 166 167
             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
168
 ENDDO
169 170 171

  END SUBROUTINE hsvac
END MODULE m_hsvac