vgen_xcpot.F90 7.16 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 61
    TYPE(t_results),   INTENT(INOUT), OPTIONAL :: results

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


68
    CALL exc%init_potden_types(stars,atoms,sphhar,vacuum,1,.false.,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 71
    IF (PRESENT(results)) THEN
       CALL veff%init(stars,atoms,sphhar,vacuum,input%jspins,.FALSE.,1)
72
       ALLOCATE(veff%pw_w,mold=veff%pw)
Daniel Wortmann's avatar
Daniel Wortmann committed
73
    ENDIF
74

75
    ! exchange correlation potential
76

77
    ! vacuum region
78 79 80 81 82 83 84 85 86 87
    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
88
                CALL vvacxc(ifftd2,stars,vacuum,xcpot,input,noco,Den,vTot,exc)
89 90
             ELSE
                CALL judft_error("OneD broken")
91 92 93 94
                ! 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
95 96 97
          ELSE      ! GGA
             IF (oneD%odi%d1) THEN
                CALL judft_error("OneD broken")
98 99 100
                ! CALL vvacxcg(ifftd2,stars,vacuum,noco,oneD,&
                !              cell,xcpot,input,obsolete,workDen, ichsmrg,&
                !              vTot%vacxy,vTot%vacz,rhmn, exc%vacxy,exc%vacz)
101 102

             ELSE
103
                CALL vvacxcg(ifftd2,stars,vacuum,noco,oneD,cell,xcpot,input,obsolete,Den,vTot,exc)
104 105 106 107
             END IF
          END IF
          CALL timestop("Vxc in vacuum")
       END IF
108 109

       ! interstitial region
110 111
       CALL timestart("Vxc in interstitial")

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

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

124
       CALL timestop("Vxc in interstitial")
125
    END IF !irank==0
126 127


128 129 130
    !
    !     ------------------------------------------
    !     ----> muffin tin spheres region
131 132 133 134

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

136 137 138
    CALL vmt_xc(DIMENSION,mpi,sphhar,atoms, den,xcpot,input,sym,&
         obsolete, vTot,vx,exc)
    
139

140
    !
141

142
    ! add MT EXX potential to vr
143 144 145
    IF (mpi%irank == 0) THEN
       CALL timestop ("Vxc in MT")

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

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

159
          veff = vTot
160
          IF(xcpot%is_hybrid().AND.hybrid%l_subvxc) THEN
Matthias Redies's avatar
Matthias Redies committed
161
            DO ispin = 1, input%jspins
162
                CALL convol(stars,vx%pw_w(:,ispin),vx%pw(:,ispin),stars%ufft)
Matthias Redies's avatar
Matthias Redies committed
163 164 165 166 167 168 169
            END DO
            veff%pw   = vTot%pw   - xcpot%get_exchange_weight() * vx%pw
            veff%pw_w = vTot%pw_w - xcpot%get_exchange_weight() * vx%pw_w
            veff%mt   = vTot%mt   - xcpot%get_exchange_weight() * vx%mt
            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
170 171 172
          END IF

          results%te_veff = 0.0
173 174
          DO ispin = 1, input%jspins
             WRITE (6,FMT=8050) ispin
175
8050         FORMAT (/,10x,'density-effective potential integrals for spin ',i2,/)
176 177 178
             CALL int_nv(ispin,stars,vacuum,atoms,sphhar,cell,sym,input,oneD,veff,workden,results%te_veff)
          END DO

179 180
          WRITE (6,FMT=8060) results%te_veff
8060      FORMAT (/,10x,'total density-effective potential integral :', t40,f20.10)
181 182 183 184

          ! CALCULATE THE INTEGRAL OF n*exc

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

187 188
          WRITE (6,FMT=8070)
8070      FORMAT (/,10x,'charge density-energy density integrals',/)
189

190
          results%te_exc = 0.0
191
          CALL int_nv(1,stars,vacuum,atoms,sphhar,cell,sym,input,oneD,exc,workDen,results%te_exc)
192
          WRITE (6,FMT=8080) results%te_exc
193

194 195
8080      FORMAT (/,10x,'total charge density-energy density integral :', t40,f20.10)
       END IF
196
    END IF ! mpi%irank == 0
197

198
  END SUBROUTINE vgen_xcpot
199

200
END MODULE m_vgen_xcpot