vvacxc.f90 7.68 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37
MODULE m_vvacxc
  !     ********************************************************************
  !     calculates 2-dim star function coefficients of exchange-correlation*
  !     potential in the vacuum regions  and adds them to the corresponding*
  !     coeffs of the coulomb potential            c.l.fu, r.podloucky     *
  !     ********************************************************************
CONTAINS
  SUBROUTINE vvacxc(&
       &                  ifftd2,stars,vacuum,xcpot,input,noco,&
       &                  rhtxy,rht,cdomvxy,cdomvz,&
       &                  vxy,vz,&
       &                  excxy,excz)

    !     ********************************************************************
    !     instead of vvacxcor.f: the different exchange-correlation 
    !     potentials defined through the key icorr are called through 
    !     the driver subroutine vxcall.f, subroutines vectorized
    !     in case of TOTAL = .TRUE. calculates the ex.-corr. energy
    !     density through the driver subroutine excall.f
    !     ** r.pentcheva 08.05.96
    !     ********************************************************************

    USE m_xcall, ONLY : vxcall,excall
    USE m_fft2d
    USE m_types
    IMPLICIT NONE
    TYPE(t_xcpot),INTENT(IN)   :: xcpot
    TYPE(t_input),INTENT(IN)   :: input
    TYPE(t_vacuum),INTENT(IN)  :: vacuum
    TYPE(t_noco),INTENT(IN)    :: noco
    TYPE(t_stars),INTENT(IN)   :: stars
    !     ..
    !     .. Scalar Arguments ..
    INTEGER, INTENT (IN) :: ifftd2
    !     ..
    !     .. Array Arguments ..
    REAL,    INTENT (INOUT) :: rht(vacuum%nmzd,2,input%jspins)
38
    COMPLEX, INTENT (INOUT) :: rhtxy(vacuum%nmzxyd,stars%ng2-1,2,input%jspins)
39
    COMPLEX, INTENT (INOUT) :: cdomvz(vacuum%nmzd,2) 
40
    COMPLEX, INTENT (INOUT) :: cdomvxy(vacuum%nmzxyd,stars%ng2-1,2) 
41
    REAL,    INTENT (OUT) :: excz(vacuum%nmzd,2)
42
    COMPLEX, INTENT (OUT) :: excxy(vacuum%nmzxyd,stars%ng2-1,2)
43
    REAL,    INTENT (INOUT) :: vz(vacuum%nmzd,2,input%jspins)
44
    COMPLEX, INTENT (INOUT) :: vxy(vacuum%nmzxyd,stars%ng2-1,2,input%jspins) 
45 46 47 48 49 50 51
    !     ..
    !     .. Local Scalars ..
    INTEGER :: k,js,nt,irec2,nmzdiff,ivac,ip,i 
    REAL    :: rhti
    REAL    :: chdens,magmom 
    !     ..
    !     .. Local Arrays ..
52
    COMPLEX :: fgxy(stars%ng2-1)
53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94
    REAL    :: af2(0:ifftd2-1,input%jspins),bf2(0:ifftd2-1),fgz
    REAL,ALLOCATABLE :: mx(:),my(:) 
    !     warping region
    REAL    :: vxc(0:ifftd2-1,input%jspins),exc(0:ifftd2-1)
    REAL    :: vx (0:ifftd2-1,input%jspins)
    !     beyond warping region
    REAL    :: vxcz(vacuum%nmzd,input%jspins)
    REAL    :: vxz (vacuum%nmzd,input%jspins)

    IF (noco%l_noco) THEN
       ALLOCATE (mx(0:ifftd2-1),my(0:ifftd2-1)) 
    END IF

    nt = ifftd2
    rhti = 0.
    !
    !     the charge density in vacuum is expanded in 2-dim stars on a mesh 
    !     in z-direction . The G||.ne.zero-components expand from 1 to nmzxy
    !     the G||.eq.zero-components expand from 1 to nmz
    !     first we calculate vxc in the warping region
    !
    DO  ivac = 1,vacuum%nvac
       DO ip = 1,vacuum%nmzxy
          !
          !         transform charge density to real space: 2-dim FFT
          !
          DO js = 1,input%jspins
             CALL fft2d(&
                  &               stars,&
                  &               af2(0,js),bf2,&
                  &               rht(ip,ivac,js),rhti,rhtxy(ip,1,ivac,js),&
                  &               vacuum%nmzxyd,+1)
          END DO

          IF (noco%l_noco) THEN 

             CALL fft2d(&
                  &               stars,&
                  &               mx,my, &
                  &               REAL(cdomvz(ip,ivac)),AIMAG(cdomvz(ip,ivac)),&
                  &               cdomvxy(ip,1,ivac),&
                  &               vacuum%nmzxyd,1)
95
             DO i=0,9*stars%mx1*stars%mx2-1 
96 97 98 99 100 101 102 103 104 105 106 107 108
                chdens= (af2(i,1)+af2(i,2))/2.  
                magmom= mx(i)**2 + my(i)**2 +&
                     &                ((af2(i,1)-af2(i,2))/2.)**2 
                magmom= SQRT(magmom) 
                af2(i,1)= chdens + magmom 
                af2(i,2)= chdens - magmom
             END DO

          END IF
          !
          !         calculate the exchange-correlation potential in  real space
          !
          CALL vxcall&
Daniel Wortmann's avatar
Daniel Wortmann committed
109
               &               (6,xcpot%icorr,xcpot%krla,input%jspins,&
110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140
               &                ifftd2,nt,af2,&
               &                vx,vxc)   

          DO  js = 1,input%jspins
             !
             !            ----> 2-d back fft to g space
             !
             bf2=0.0
             CALL fft2d(&
                  &                 stars,&
                  &                 vxc(0,js),bf2,&
                  &                 fgz,rhti,fgxy,&
                  &                 1,-1)
             !
             !            ----> and add vxc to coulomb potential
             !            the G||.eq.zero component is added to vz
             !
             vz(ip,ivac,js) = fgz + vz(ip,ivac,js)
             !
             !            the G||.ne.zero components are added to vxy
             !
             DO irec2 = 1,stars%ng2-1
                vxy(ip,irec2,ivac,js) = vxy(ip,irec2,ivac,js) +&
                     &                                              fgxy(irec2)
             ENDDO
          ENDDO
          !
          !i        calculate the exchange-correlation energy density in  real space
          !
          IF (input%total) THEN   
             CALL excall&
Daniel Wortmann's avatar
Daniel Wortmann committed
141
                  &                 (6,xcpot%icorr,xcpot%krla,input%jspins,&
142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182
                  &                  ifftd2,nt,af2,&
                  &                  exc)   
             !
             !     ----> 2-d back fft to g space
             !
             bf2=0.0
             CALL fft2d(&
                  &                 stars,&
                  &                 exc,bf2,&
                  &                 excz(ip,ivac),rhti,excxy(ip,1,ivac),&
                  &                 vacuum%nmzxyd,-1)
          ENDIF

       ENDDO
       !
       !        calculate vxc for z now beyond warping region 
       !
       nmzdiff = vacuum%nmz - vacuum%nmzxy

       DO k=1,nmzdiff

          DO js=1,input%jspins
             af2(k-1,js) = rht(vacuum%nmzxy+k,ivac,js)
          ENDDO

          IF (noco%l_noco) THEN

             mx(0)= REAL(cdomvz(vacuum%nmzxy+k,ivac))
             my(0)= AIMAG(cdomvz(vacuum%nmzxy+k,ivac))
             chdens= (af2(k-1,1)+af2(k-1,2))/2.
             magmom= mx(0)**2 + my(0)**2 +&
                  &               ((af2(k-1,1)-af2(k-1,2))/2.)**2
             magmom= SQRT(magmom)
             af2(k-1,1)= chdens + magmom
             af2(k-1,2)= chdens - magmom

          END IF

       ENDDO

       CALL vxcall&
Daniel Wortmann's avatar
Daniel Wortmann committed
183
            &              (6,xcpot%icorr,xcpot%krla,input%jspins,&
184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200
            &               vacuum%nmzd,nmzdiff,af2,&
            &               vxz,vxcz)
       !+gu
       DO  js=1,input%jspins
          DO k=vacuum%nmzxy+1,vacuum%nmz
             vz(k,ivac,js) = vz(k,ivac,js) + vxcz(k-vacuum%nmzxy,js)
          ENDDO
       ENDDO
       !
       WRITE (6,FMT=8020) ivac, (vz(vacuum%nmz,ivac,js),js=1,input%jspins)
       WRITE (16,FMT=8020) ivac, (vz(vacuum%nmz,ivac,js),js=1,input%jspins)
8020   FORMAT (/,5x,'vacuum zero for vacuum',i3,' = ',2f10.5)
       !
       !        calculate the ex.-corr. energy density now beyond warping region
       !
       IF (input%total) THEN
          CALL excall&
Daniel Wortmann's avatar
Daniel Wortmann committed
201
               &                   (6,xcpot%icorr,xcpot%krla,input%jspins,&
202
               &                    vacuum%nmzd,nmzdiff,af2,&
Daniel Wortmann's avatar
Daniel Wortmann committed
203
               &                    excz(vacuum%nmzxy+1,ivac))
204 205 206 207 208 209 210 211
       END IF
    ENDDO
    IF (noco%l_noco) THEN 
       DEALLOCATE (mx,my)
    END IF

  END SUBROUTINE vvacxc
END MODULE m_vvacxc