vgen_xcpot.F90 7.47 KB
Newer Older
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_vgen_xcpot
7

8
  USE m_juDFT
9

10
CONTAINS
11 12

  SUBROUTINE vgen_xcpot(hybrid,input,xcpot,dimension,atoms,sphhar,stars,vacuum,sym,&
13
                        obsolete,cell,oneD,sliceplot,mpi,noco,den,denRot,vTot,vx,results)
14

15 16 17 18 19 20 21 22 23
    !     ***********************************************************
    !     FLAPW potential generator                           *
    !     ***********************************************************
    !     calculates the density-potential integrals needed for the
    !     total energy
    !     TE_VCOUL  :   charge density-coulomb potential integral
    !     TE_VEFF:   charge density-effective potential integral
    !     TE_EXC :   charge density-ex-corr.energy density integral
    !     ***********************************************************
24 25

    USE m_types
26 27
    USE m_constants
    USE m_intnv
28
    USE m_vmt_xc
29 30
    USE m_vvacxc
    USE m_vvacxcg
31
    USE m_vis_xc
32
    USE m_checkdopall
33 34 35
    USE m_cdn_io
    USE m_convol

36
    IMPLICIT NONE
37

38

39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54
    CLASS(t_xcpot),    INTENT(IN)              :: xcpot
    TYPE(t_hybrid),    INTENT(IN)              :: hybrid
    TYPE(t_mpi),       INTENT(IN)              :: mpi
    TYPE(t_dimension), INTENT(IN)              :: dimension
    TYPE(t_oneD),      INTENT(IN)              :: oneD
    TYPE(t_obsolete),  INTENT(IN)              :: obsolete
    TYPE(t_sliceplot), INTENT(IN)              :: sliceplot
    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_sphhar),    INTENT(IN)              :: sphhar
    TYPE(t_atoms),     INTENT(IN)              :: atoms 
    TYPE(t_potden),    INTENT(IN)              :: den,denRot
55
    TYPE(t_potden),    INTENT(INOUT)           :: vTot,vx
56 57 58 59 60
    TYPE(t_results),   INTENT(INOUT), OPTIONAL :: results

    ! Local type instances
    TYPE(t_potden) :: workDen,exc,veff
    ! Local Scalars
61
    INTEGER ifftd,ifftd2,ifftxc3d,ispin,i
62 63 64 65 66 67
#ifdef CPP_MPI
    include 'mpif.h'
    integer:: ierr
#endif


Daniel Wortmann's avatar
Daniel Wortmann committed
68
    CALL exc%init(stars,atoms,sphhar,vacuum,noco,1,1) !one spin only
69
    ALLOCATE(exc%pw_w(stars%ng3,1));exc%pw_w=0.0
Daniel Wortmann's avatar
Daniel Wortmann committed
70
    IF (PRESENT(results)) THEN
Daniel Wortmann's avatar
Daniel Wortmann committed
71
       CALL veff%init(stars,atoms,sphhar,vacuum,noco,input%jspins,1)
72
#ifndef CPP_OLDINTEL
73
       ALLOCATE(veff%pw_w,mold=veff%pw)
74 75 76
#else
       ALLOCATE( veff%pw_w(size(veff%pw,1),size(veff%pw,2)) )
#endif
Daniel Wortmann's avatar
Daniel Wortmann committed
77
    ENDIF
78

79
    ! exchange correlation potential
80

81
    ! vacuum region
82 83 84 85 86 87 88 89 90 91
    IF (mpi%irank == 0) THEN
       IF (input%film) THEN
          CALL timestart("Vxc in vacuum")

          ifftd2 = 9*stars%mx1*stars%mx2
          IF (oneD%odi%d1) ifftd2 = 9*stars%mx3*oneD%odi%M

          IF (.NOT.xcpot%is_gga()) THEN  ! LDA

             IF (.NOT.oneD%odi%d1) THEN
92
                CALL vvacxc(ifftd2,stars,vacuum,xcpot,input,noco,Den,vTot,exc)
93 94
             ELSE
                CALL judft_error("OneD broken")
95 96 97 98
                ! CALL vvacxc(stars,oneD%M,vacuum,odi%n2d,dimension,ifftd2,&
                !             xcpot,input,odi%nq2,odi%nst2,den,noco,odi%kimax2%igf,&
                !             odl%pgf,vTot%vacxy,vTot%vacz,excxy,excz)
             END IF
99 100 101
          ELSE      ! GGA
             IF (oneD%odi%d1) THEN
                CALL judft_error("OneD broken")
102 103 104
                ! CALL vvacxcg(ifftd2,stars,vacuum,noco,oneD,&
                !              cell,xcpot,input,obsolete,workDen, ichsmrg,&
                !              vTot%vacxy,vTot%vacz,rhmn, exc%vacxy,exc%vacz)
105 106

             ELSE
107
                CALL vvacxcg(ifftd2,stars,vacuum,noco,oneD,cell,xcpot,input,obsolete,Den,vTot,exc)
108 109 110 111
             END IF
          END IF
          CALL timestop("Vxc in vacuum")
       END IF
112 113

       ! interstitial region
114 115
       CALL timestart("Vxc in interstitial")

116
       
117 118
       IF ( (.NOT. obsolete%lwb) .OR. ( .not.xcpot%is_gga() ) ) THEN
          ! no White-Bird-trick
119
          CALL vis_xc(stars,sym,cell,den,xcpot,input,noco,vTot,vx,exc)
120 121 122 123

       ELSE
          ! White-Bird-trick
          WRITE(6,'(a)') "W+B trick cancelled out. visxcwb uses at present common block cpgft3.",&
124
                         "visxcwb needs to be reprogrammed according to visxcg.f"
125 126
          CALL juDFT_error("visxcwb",calledby ="vgen")
       END IF
127

128
       CALL timestop("Vxc in interstitial")
129
    END IF !irank==0
130 131


132 133 134
    !
    !     ------------------------------------------
    !     ----> muffin tin spheres region
135 136 137 138

    IF (mpi%irank == 0) THEN
       CALL timestart ("Vxc in MT")
    END IF
139

140 141 142
    CALL vmt_xc(DIMENSION,mpi,sphhar,atoms, den,xcpot,input,sym,&
         obsolete, vTot,vx,exc)
    
143

144
    !
145

146
    ! add MT EXX potential to vr
147 148 149
    IF (mpi%irank == 0) THEN
       CALL timestop ("Vxc in MT")

150 151
       ! check continuity of total potential
       IF (input%vchk) CALL checkDOPAll(input,dimension,sphhar,stars,atoms,sym,vacuum,oneD,cell,vTot,1)
152

153
       ! TOTAL 
154
       IF (PRESENT(results)) THEN
155 156
          ! CALCULATE THE INTEGRAL OF n1*Veff1 + n2*Veff2
          ! Veff = Vcoulomb + Vxc
157 158
          IF (noco%l_noco) THEN 
             workDen = denRot
159
          ELSE
160 161
             workden = den
          END IF
162

163
          veff = vTot
164
          IF(xcpot%is_hybrid().AND.hybrid%l_subvxc) THEN
165
             DO ispin = 1, input%jspins
166
                CALL convol(stars,vx%pw_w(:,ispin),vx%pw(:,ispin),stars%ufft)
167
             END DO
168
             veff%pw   = vTot%pw   - xcpot%get_exchange_weight() * vx%pw
169 170
             veff%pw_w = vTot%pw_w - xcpot%get_exchange_weight() * vx%pw_w
             veff%mt   = vTot%mt   - xcpot%get_exchange_weight() * vx%mt
171 172 173
             exc%pw   = exc%pw   - xcpot%get_exchange_weight() * exc%pw
             exc%pw_w = exc%pw_w - xcpot%get_exchange_weight() * exc%pw_w
             exc%mt   = exc%mt   - xcpot%get_exchange_weight() * exc%mt
174 175
          END IF

176 177 178 179 180 181 182
          DO ispin = 1, input%jspins
             DO i = 1, stars%ng3
                vx%pw(i,ispin) = vx%pw(i,ispin) / stars%nstr(i)
                vx%pw_w(i,ispin) = vx%pw_w(i,ispin) / stars%nstr(i)
             END DO
          END DO

183
          results%te_veff = 0.0
184 185
          DO ispin = 1, input%jspins
             WRITE (6,FMT=8050) ispin
186
8050         FORMAT (/,10x,'density-effective potential integrals for spin ',i2,/)
187 188 189
             CALL int_nv(ispin,stars,vacuum,atoms,sphhar,cell,sym,input,oneD,veff,workden,results%te_veff)
          END DO

190 191
          WRITE (6,FMT=8060) results%te_veff
8060      FORMAT (/,10x,'total density-effective potential integral :', t40,f20.10)
192 193 194 195

          ! CALCULATE THE INTEGRAL OF n*exc

          ! perform spin summation of charge densities for the calculation of Exc
196
          CALL  workden%sum_both_spin()
197

198 199
          WRITE (6,FMT=8070)
8070      FORMAT (/,10x,'charge density-energy density integrals',/)
200

201
          results%te_exc = 0.0
202
          CALL int_nv(1,stars,vacuum,atoms,sphhar,cell,sym,input,oneD,exc,workDen,results%te_exc)
203
          WRITE (6,FMT=8080) results%te_exc
204

205 206
8080      FORMAT (/,10x,'total charge density-energy density integral :', t40,f20.10)
       END IF
207
    END IF ! mpi%irank == 0
208

209
  END SUBROUTINE vgen_xcpot
210

211
END MODULE m_vgen_xcpot