vgen_xcpot.F90 7.37 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 13 14

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

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 28 29 30 31 32 33 34
    USE m_constants
    USE m_intnv
    USE m_vmtxcg
    USE m_vmtxc
    USE m_vvacxc
    USE m_vvacxcg
    USE m_visxc
    USE m_visxcg
    USE m_checkdopall
35 36 37
    USE m_cdn_io
    USE m_convol

38
    IMPLICIT NONE
39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62

    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
    TYPE(t_potden),    INTENT(INOUT)           :: vTot,vXC
    TYPE(t_results),   INTENT(INOUT), OPTIONAL :: results

    ! Local type instances
    TYPE(t_potden) :: workDen,exc,veff
    ! Local Scalars
    INTEGER ifftd,ifftd2,ifftxc3d,ispin
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
Daniel Wortmann's avatar
Daniel Wortmann committed
69 70 71
    ALLOCATE(exc%pw_w(stars%ng3,1))
    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
       ifftd = 27*stars%mx1*stars%mx2*stars%mx3
113

114
       IF ((.NOT.obsolete%lwb).OR.(.NOT.xcpot%is_gga())) THEN
115 116 117
          ! no White-Bird-trick
          ifftxc3d = stars%kxc1_fft*stars%kxc2_fft*stars%kxc3_fft

118 119
          IF (.NOT.xcpot%is_gga()) THEN  ! LDA
             CALL visxc(ifftd,stars,noco,xcpot,input,den,vTot,vXC,exc)
120
          ELSE ! GGA
121
             CALL visxcg(ifftd,stars,sym,ifftxc3d,cell,den,xcpot,input,obsolete,noco,vTot,vXC,exc)
122 123 124 125
          END IF
       ELSE
          ! White-Bird-trick
          WRITE(6,'(a)') "W+B trick cancelled out. visxcwb uses at present common block cpgft3.",&
126
                         "visxcwb needs to be reprogrammed according to visxcg.f"
127 128
          CALL juDFT_error("visxcwb",calledby ="vgen")
       END IF
129

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

133 134
    ! muffin tin spheres region
    IF (mpi%irank == 0) CALL timestart ("Vxc in MT")
135
#ifdef CPP_MPI
136
    CALL MPI_BCAST(den%mt,atoms%jmtd*(1+sphhar%nlhd)*atoms%ntype*dimension%jspd,MPI_DOUBLE_PRECISION,0,mpi%mpi_comm,ierr)
137
#endif
138
    IF (xcpot%is_gga()) THEN
139
       CALL vmtxcg(dimension,mpi,sphhar,atoms,den,xcpot,input,sym,obsolete,vTot,vXC,exc)
140
    ELSE
141
       CALL vmtxc(DIMENSION,sphhar,atoms,den,xcpot,input,sym,vTot,exc,vXC)
142
    ENDIF
143

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

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

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

          veff = vTot
          IF(xcpot%is_hybrid()) THEN
             DO ispin = 1, input%jspins
                CALL convol(stars,vXC%pw_w(:,ispin),vXC%pw(:,ispin),stars%ufft)
             END DO
             veff%pw   = vTot%pw   - xcpot%get_exchange_weight() * vXC%pw
             veff%pw_w = vTot%pw_w - xcpot%get_exchange_weight() * vXC%pw_w
             veff%mt   = vTot%mt   - xcpot%get_exchange_weight() * vXC%mt
169 170 171
          END IF

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

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

          ! CALCULATE THE INTEGRAL OF n*exc

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

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

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

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

197
  END SUBROUTINE vgen_xcpot
198

199
END MODULE m_vgen_xcpot