vvacxc.f90 5.92 KB
Newer Older
1 2 3 4 5 6 7
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
8
  SUBROUTINE vvacxc(ifftd2,stars,vacuum,xcpot,input,noco,den,vxc,exc)
9 10 11 12 13 14 15 16 17 18 19 20 21

    !     ********************************************************************
    !     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_fft2d
    USE m_types
    IMPLICIT NONE
22
    CLASS(t_xcpot),INTENT(IN)     :: xcpot
Gregor Michalicek's avatar
Gregor Michalicek committed
23 24 25 26
    TYPE(t_input),INTENT(IN)     :: input
    TYPE(t_vacuum),INTENT(IN)    :: vacuum
    TYPE(t_noco),INTENT(IN)      :: noco
    TYPE(t_stars),INTENT(IN)     :: stars
27 28
    TYPE(t_potden),INTENT(IN)    :: den
    TYPE(t_potden),INTENT(INOUT) :: vxc,exc
29 30 31 32
    !     ..
    !     .. Scalar Arguments ..
    INTEGER, INTENT (IN) :: ifftd2
    !     ..
33
 
34 35 36 37 38 39 40
    !     ..
    !     .. Local Scalars ..
    INTEGER :: k,js,nt,irec2,nmzdiff,ivac,ip,i 
    REAL    :: rhti
    REAL    :: chdens,magmom 
    !     ..
    !     .. Local Arrays ..
41
    COMPLEX :: fgxy(stars%ng2-1)
42 43 44
    REAL    :: af2(0:ifftd2-1,input%jspins),bf2(0:ifftd2-1),fgz
    REAL,ALLOCATABLE :: mx(:),my(:) 
    !     warping region
45 46
    REAL    :: v_xc(0:ifftd2-1,input%jspins),e_xc(0:ifftd2-1)
    REAL    :: v_x (0:ifftd2-1,input%jspins)
47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68
    !     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
69 70
             CALL fft2d(stars, af2(0,js),bf2, den%vacz(ip,ivac,js),rhti,&
                  den%vacxy(ip,1,ivac,js), vacuum%nmzxyd,+1)
71 72 73 74
          END DO

          IF (noco%l_noco) THEN 

75 76
             CALL fft2d(stars, mx,my, den%vacz(ip,ivac,3),den%vacz(ip,ivac,4),&
                                 den%vacxy(ip,1,ivac,3), vacuum%nmzxyd,1)
77
             DO i=0,9*stars%mx1*stars%mx2-1 
78
                chdens= (af2(i,1)+af2(i,2))/2.  
79
                magmom= mx(i)**2 + my(i)**2 + ((af2(i,1)-af2(i,2))/2.)**2 
80 81 82 83 84 85 86 87 88
                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
          !
89
          CALL xcpot%get_vxc(input%jspins,af2,v_xc,v_x)
90 91 92 93 94 95

          DO  js = 1,input%jspins
             !
             !            ----> 2-d back fft to g space
             !
             bf2=0.0
96
             CALL fft2d(stars, v_xc(0,js),bf2, fgz,rhti,fgxy, 1,-1)
97 98 99 100
             !
             !            ----> and add vxc to coulomb potential
             !            the G||.eq.zero component is added to vz
             !
101
             vxc%vacz(ip,ivac,js) = fgz + vxc%vacz(ip,ivac,js)
102
             !
103
             !            the G||.ne.zero components are added to vxc%vacxy
104 105
             !
             DO irec2 = 1,stars%ng2-1
106
                vxc%vacxy(ip,irec2,ivac,js) = vxc%vacxy(ip,irec2,ivac,js) + fgxy(irec2)
107 108 109 110 111
             ENDDO
          ENDDO
          !
          !i        calculate the exchange-correlation energy density in  real space
          !
112
          IF (ALLOCATED(exc%vacz)) THEN
113
             call xcpot%get_exc(input%jspins,af2,e_xc, mt_call=.False.)
114 115 116 117
             !
             !     ----> 2-d back fft to g space
             !
             bf2=0.0
118
             CALL fft2d(stars, e_xc,bf2, exc%vacz(ip,ivac,1),rhti,exc%vacxy(ip,1,ivac,1), vacuum%nmzxyd,-1)
119 120 121 122 123 124 125 126 127 128 129
          ENDIF

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

       DO k=1,nmzdiff

          DO js=1,input%jspins
Gregor Michalicek's avatar
Gregor Michalicek committed
130
             af2(k-1,js) = den%vacz(vacuum%nmzxy+k,ivac,js)
131 132 133 134
          ENDDO

          IF (noco%l_noco) THEN

Gregor Michalicek's avatar
Gregor Michalicek committed
135 136
             mx(0)= den%vacz(vacuum%nmzxy+k,ivac,3)
             my(0)= den%vacz(vacuum%nmzxy+k,ivac,4)
137
             chdens= (af2(k-1,1)+af2(k-1,2))/2.
138
             magmom= mx(0)**2 + my(0)**2 + ((af2(k-1,1)-af2(k-1,2))/2.)**2
139 140 141 142 143 144 145 146
             magmom= SQRT(magmom)
             af2(k-1,1)= chdens + magmom
             af2(k-1,2)= chdens - magmom

          END IF

       ENDDO

147
       CALL xcpot%get_vxc(input%jspins,af2(0:nmzdiff-1,:),vxz,vxcz)
148
      
149 150 151
       !+gu
       DO  js=1,input%jspins
          DO k=vacuum%nmzxy+1,vacuum%nmz
152
             vxc%vacz(k,ivac,js) = vxc%vacz(k,ivac,js) + vxcz(k-vacuum%nmzxy,js)
153 154 155
          ENDDO
       ENDDO
       !
156
       WRITE (6,FMT=8020) ivac, (vxc%vacz(vacuum%nmz,ivac,js),js=1,input%jspins)
157 158 159 160
8020   FORMAT (/,5x,'vacuum zero for vacuum',i3,' = ',2f10.5)
       !
       !        calculate the ex.-corr. energy density now beyond warping region
       !
161
       IF (ALLOCATED(exc%vacz)) THEN
162
          CALL xcpot%get_exc(input%jspins,af2(0:nmzdiff-1,:),exc%vacz(vacuum%nmzxy+1:,ivac,1), mt_call=.False.)
163 164 165 166 167 168 169 170
       END IF
    ENDDO
    IF (noco%l_noco) THEN 
       DEALLOCATE (mx,my)
    END IF

  END SUBROUTINE vvacxc
END MODULE m_vvacxc