vgen.F90 32.5 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 7 8
MODULE m_vgen
  USE m_juDFT
CONTAINS
9
  SUBROUTINE vgen(hybrid,reap,input,xcpot,DIMENSION, atoms,sphhar,stars,&
10
       vacuum,sym, obsolete,cell,oneD,sliceplot,mpi, results,noco,den,vTot,vx,vCoul)
11
    !     ***********************************************************
12
    !     FLAPW potential generator                           *
13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
    !     ***********************************************************
    !     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
    !     ***********************************************************
#include"cpp_double.h"
    USE m_constants
    USE m_vmts
    USE m_intnv
    USE m_vmtxcg
    USE m_vmtxc
    USE m_vvacxc
    USE m_vvacxcg
    USE m_visxc
    USE m_visxcg
    USE m_vvac
    USE m_vvacis
    USE m_vvacxy
    USE m_vintcz
    USE m_checkdop
    USE m_wrtdop
36
    USE m_cdn_io
37 38 39 40 41 42 43 44 45 46 47 48 49 50
    USE m_qfix
    USE m_types
    USE m_od_vvac
    USE m_od_vvacis
    USE m_cylpts
    USE m_convol
    USE m_xyavden
    USE m_psqpw
    USE m_potmod
    USE m_intgr,         ONLY : intgr3
    USE m_cfft
    USE m_sphpts
    USE m_points
    USE m_fleur_vdw
51 52 53
#ifdef CPP_MPI
    USE m_mpi_bc_potden
#endif
54 55 56
    IMPLICIT NONE
    TYPE(t_results),INTENT(INOUT)   :: results
    TYPE(t_xcpot),INTENT(IN)        :: xcpot
57
    TYPE(t_hybrid),INTENT(IN)       :: hybrid
58 59 60 61 62 63 64 65 66 67 68 69 70
    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(INOUT)     :: input  !efield can be modified
    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(INOUT)     :: atoms !vr0 is updated
71
    TYPE(t_potden), INTENT(IN)      :: den
72
    TYPE(t_potden),INTENT(INOUT)    :: vTot,vx,vCoul
73 74
    !     ..
    !     .. Scalar Arguments ..
75 76 77 78
    LOGICAL, INTENT (IN) :: reap


    !     .. Local type instances ..
79
    TYPE(t_potden)    :: workDen
80 81 82
    !     .. Local Scalars ..
    COMPLEX vintcza,xint,rhobar
    INTEGER i,i3,irec2,irec3,ivac,j,js,k,k3,lh,n,nzst1
83
    INTEGER imz,imzxy,ichsmrg,ivfft,npd 
Daniel Wortmann's avatar
Daniel Wortmann committed
84
    INTEGER ifftd,ifftd2, ifftxc3d,datend
85
    INTEGER itypsym,itype,jsp,l,nat
86
    !      INTEGER i_sm,n_sm,i_sta,i_end
87
    REAL ani,g3,signum,z,rhmn,mfie,fermiEnergyTemp
Daniel Wortmann's avatar
Daniel Wortmann committed
88
    REAL sig1dh,vz1dh,zat_l(atoms%ntype),rdum,dpdot ! ,delta,deltb,corr
89
    LOGICAL l_pottot,l_vdw,l_qfix
90 91 92 93 94 95 96
    LOGICAL exi
    LOGICAL, PARAMETER :: l_xyav=.FALSE.
    !     ..
    !     .. Local Arrays ..
    COMPLEX, ALLOCATABLE :: alphm(:,:)
    COMPLEX, ALLOCATABLE :: excpw(:),excxy(:,:,:),vpw_w(:,:),psq(:)
    REAL,    ALLOCATABLE :: vbar(:),af1(:),bf1(:),xp(:,:)
97 98
    REAL,    ALLOCATABLE :: rhoc(:,:,:),rhoc_vx(:)
    REAL,    ALLOCATABLE :: tec(:,:), qintc(:,:)
99
    !.....potential
Daniel Wortmann's avatar
Daniel Wortmann committed
100 101
    REAL,    ALLOCATABLE :: vr_exx(:,:,:,:)
    REAL,    ALLOCATABLE :: veffr(:,:,:,:)
102
    COMPLEX, ALLOCATABLE :: vpw_exx(:,:),vpw_wexx(:,:)
Daniel Wortmann's avatar
Daniel Wortmann committed
103
    COMPLEX, ALLOCATABLE :: vxpw_w(:,:),veffpw_w(:,:)
104 105 106

    !.....energy density
    REAL,    ALLOCATABLE :: excz(:,:),excr(:,:,:)
107 108 109 110
#ifdef CPP_MPI
    include 'mpif.h'
    integer:: ierr
#endif
111 112 113 114 115 116 117 118 119 120 121 122 123 124
    !
    ! if you want to calculate potential gradients
    !
    !pg      COMPLEX vlm(-lmaxd:lmaxd,0:lmaxd,ntypd)
    !pg      REAL fint,f(jmtd)
    !pg      INTEGER ns, nl, l, jm, m, mb
    !     ..
    !     ..
    !     ..
    !     ----> note the following conventions:
    !     ivac=2: lower (negative z) vacuum
    !     ivac=1: upper (positive z) vacuum
    !     units: hartrees
    !
125

126 127 128 129
    CALL vTot%resetPotDen()
    CALL vCoul%resetPotDen()
    CALL vx%resetPotDen()

130 131
    ALLOCATE ( alphm(stars%ng2,2),excpw(stars%ng3),excxy(vacuum%nmzxyd,oneD%odi%n2d-1,2),&
         vbar(dimension%jspd),af1(3*stars%mx3),bf1(3*stars%mx3),xp(3,dimension%nspd),&
Daniel Wortmann's avatar
Daniel Wortmann committed
132
         vpw_exx(stars%ng3,dimension%jspd),vpw_wexx(stars%ng3,dimension%jspd),&
Daniel Wortmann's avatar
Daniel Wortmann committed
133
         excz(vacuum%nmzd,2),excr(atoms%jmtd,0:sphhar%nlhd,atoms%ntype),&
Daniel Wortmann's avatar
Daniel Wortmann committed
134
         vpw_w(stars%ng3,dimension%jspd),vxpw_w(stars%ng3,dimension%jspd),psq(stars%ng3) )
135
    vTot%iter = den%iter
136

137
    CALL workDen%init(stars,atoms,sphhar,vacuum,noco,oneD,input%jspins,.FALSE.,POTDEN_TYPE_DEN)
138

139
    workDen = den
140 141 142 143 144 145 146 147 148
    IF (mpi%irank == 0) THEN
       !
       ! --  total = .false. and reap = .false. means, that we now calculate
       !     the potential from the output density
       !
       IF ((.NOT.input%total).AND.(.NOT.reap)) THEN
          IF (noco%l_noco) THEN
             CALL juDFT_error("vgen:1",calledby ="vgen")
          ENDIF
149 150
       END IF

151 152 153 154 155 156
       WRITE (6,FMT=8000)
8000   FORMAT (/,/,t10,' p o t e n t i a l   g e n e r a t o r',/)
       vxpw_w = CMPLX(0.,0.)
       !     ---> perform spin summation of charge densities
       !     ---> for the calculation of the coulomb potentials
       IF ( (input%jspins.EQ.2).AND.(.NOT.l_xyav) ) THEN
157 158
          workDen%mt(:,0:,:,1)=workDen%mt(:,0:,:,1)+workDen%mt(:,0:,:,2)
          workDen%pw(:,1)=workDen%pw(:,1)+workDen%pw(:,2)
159
          IF (input%film) THEN
160 161
             workDen%vacxy(:,:,:,1)=workDen%vacxy(:,:,:,1)+workDen%vacxy(:,:,:,2)
             workDen%vacz(:,:,1) = workDen%vacz(:,:,1) + workDen%vacz(:,:,2)
162 163 164 165 166 167 168 169 170 171 172 173
          END IF
       END IF
       !
       !     ************** coulomb potential ***********************
       IF (l_xyav) THEN
          zat_l(:) = 0.          ! for xy-averaged densities do not
       ELSE                     ! include nuclear charge
          zat_l(:) = atoms%zatom(:)
       ENDIF

       !     ----> create pesudo-charge density coefficients
    ENDIF ! (mpi%irank == 0)
Daniel Wortmann's avatar
Daniel Wortmann committed
174 175

#ifdef CPP_MPI
176
    CALL MPI_BCAST(workDen%mt,atoms%jmtd*(1+sphhar%nlhd)*atoms%ntype*dimension%jspd,MPI_DOUBLE_PRECISION,0,mpi%mpi_comm,ierr)
Daniel Wortmann's avatar
Daniel Wortmann committed
177 178
#endif

179 180
    CALL timestart("psqpw")      
    CALL psqpw(mpi, atoms,sphhar,stars,vacuum, dimension,cell,input,sym,oneD,&
181
         workDen%pw,workDen%mt,workDen%vacz,l_xyav, psq)
182 183 184 185
    CALL timestop("psqpw")
    IF (mpi%irank == 0) THEN

       IF (l_xyav) THEN        ! write out xy-averaged density & stop
186
          CALL xy_av_den(stars,vacuum, cell,psq,workDen%vacz)
187 188 189 190 191 192 193 194 195
          CALL juDFT_error("xy-averaged density calculated", calledby ="vgen")
       ENDIF

       !     ------------------------------------------
       IF (oneD%odi%d1) THEN
          !-odim
          CALL timestart("coulomb potential")

          !---> generates the m=0,gz=0 component of the vacuum potential
196
          CALL od_vvac(stars,vacuum,cell, psq,workDen%vacz, vTot%vacz)
197 198 199 200 201 202 203

          !---> generation of the vacuum warped potential components and
          !---> interstitial pw potential
          !---> vvacxy_5.F is a symmetrized potential generator

          CALL od_vvacis(oneD%odi%n2d,dimension,vacuum,oneD%odi%nq2,&
               oneD%odi%kv,cell,oneD%odi%M,stars,oneD%odi%nst2,&
204
               oneD, workDen%vacz,workDen%vacxy,psq,vTot%vacz,sym, vTot%vacxy,vTot%pw)
205 206 207 208 209 210 211
          CALL timestop("coulomb potential")

          !+odim
       ELSEIF (input%film .AND. .NOT.oneD%odi%d1) THEN
          !     ----> potential in the  vacuum  region
          !       
          CALL timestart("p vac") 
212 213
          CALL vvac(vacuum,stars, cell,sym,input, psq,workDen%vacz, vTot%vacz,rhobar,sig1dh,vz1dh)
          CALL vvacis(stars,vacuum, sym,cell, psq, input, vTot%vacxy)
214

215
          CALL vvacxy(stars,vacuum,cell,sym,input, workDen%vacxy, vTot%vacxy, alphm)
216 217 218 219 220 221 222 223 224
          CALL timestop("p vac")
       END IF
       !     ------------------------------------------
       !     ----> potential in the  interstitial  region
       CALL timestart("p int")
       WRITE (6,FMT=8010)
8010   FORMAT (/,5x,'coulomb potential in the interstitial region:')
       IF (input%film .AND. .NOT.oneD%odi%d1) THEN
          !           -----> create v(z) for each 2-d reciprocal vector
225
          ivfft =  3*stars%mx3 
226 227 228 229 230 231 232 233 234
          !         ivfft = 2*mx3 - 1
          ani = 1.0/REAL(ivfft)
          DO  irec2 = 1,stars%ng2
             i = 0
             DO i3 = 0,ivfft - 1
                i = i + 1
                z = cell%amat(3,3)*i3*ani
                IF (z.GT.cell%amat(3,3)/2.) z = z - cell%amat(3,3)
                vintcza = vintcz(stars,vacuum,cell,sym,input,&
235
                     z,irec2, psq,vTot%vacxy,vTot%vacz,rhobar,sig1dh,vz1dh,alphm)
236 237 238 239 240 241 242 243 244
                af1(i) = REAL(vintcza)
                bf1(i) = AIMAG(vintcza)
             ENDDO
             !                z = (i_sm-1)*ani
             !                IF (z > 0.5) z = z - 1.0
             !                af1(i_sm) = af1(i_sm) + z * delta
             !                bf1(i_sm) = bf1(i_sm) + z * deltb
             !              ENDDO
             !            ENDIF
245
             !        --> 1-d fourier transform and store the coefficients in vTot%pw( ,1)
246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265
             CALL cfft(af1,bf1,ivfft,ivfft,ivfft,-1)
             !            delta = ivfft * delta * 2 / fpi ! * amat(3,3)**2 * ani
             i = 0
             DO  i3 = 0,ivfft - 1
                k3 = i3
                IF (k3 > FLOOR(ivfft/2.) ) k3 = k3 - ivfft
                i = i + 1
                IF ((k3.GE.-stars%mx3).AND.(k3.LE.stars%mx3)) THEN
                   irec3 = stars%ig(stars%kv2(1,irec2),stars%kv2(2,irec2),k3)

                   !                 IF ( (irec2 == 1).AND.(i3 > 0) ) THEN                 ! smooth potential
                   !                   corr = 2.0*mod(abs(k3),2) - 1.0
                   !                   bf1(i) = bf1(i) + delta * corr / k3
                   !                 ENDIF

                   !       ----> only stars within g_max sphere (shz oct.97)
                   IF (irec3.NE.0) THEN
                      !
                      xint = CMPLX(af1(i),bf1(i))*ani
                      nzst1 = stars%nstr(irec3)/stars%nstr2(irec2)
266
                      vTot%pw(irec3,1) = vTot%pw(irec3,1) + xint/nzst1
267 268 269 270 271
                   END IF
                ENDIF
             ENDDO
          ENDDO
       ELSEIF (.NOT.input%film) THEN
272 273
          vTot%pw(1,1) = CMPLX(0.0,0.0)
          vTot%pw(2:stars%ng3,1)=fpi_const*psq(2:stars%ng3)/(stars%sk3(2:stars%ng3)*stars%sk3(2:stars%ng3))       
274 275 276 277 278 279 280 281 282
       END IF

       CALL timestop("p int")

    ENDIF ! mpi%irank == 0
    !     --------------------------------------------
    !     ---> potential in the muffin-tin spheres

    CALL timestart("p vmts")
283
    CALL vmts(mpi, stars,sphhar,atoms, sym,cell,oneD, vTot%pw,workDen%mt, vTot%mt)
284 285 286 287 288 289 290 291 292 293 294 295 296 297 298
    !     --------------------------------------------
    CALL timestop("p vmts")
    IF (mpi%irank == 0) THEN
       !     ---> check continuity of coulomb potential

       IF (input%vchk) THEN
          CALL timestart("checking")
          !           ----> vacuum boundaries
          IF (input%film .AND. .NOT.oneD%odi%d1) THEN
             npd = MIN(dimension%nspd,25)
             CALL points(xp,npd)
             DO ivac = 1,vacuum%nvac
                signum = 3. - 2.*ivac
                xp(3,:npd) = signum*cell%z1/cell%amat(3,3)
                CALL checkdop(xp,npd,0,0,ivac,1,1,.FALSE.,dimension,atoms, sphhar,stars,sym,&
299
                     vacuum,cell,oneD, vTot%pw,vTot%mt,vTot%vacxy,vTot%vacz)
300 301 302 303 304 305 306 307 308 309
             ENDDO
          ELSEIF (oneD%odi%d1) THEN
             !-odim
             npd = MIN(dimension%nspd,25)
             CALL cylpts(xp,npd,cell%z1)
             !           DO j = 1,npd
             !              xp(1,j) = xp(1,j)/amat(1,1)
             !              xp(2,j) = xp(2,j)/amat(2,2)
             !           ENDDO
             CALL checkdop(xp,npd,0,0,vacuum%nvac,1,1,.FALSE.,dimension,atoms,&
310
                  sphhar,stars,sym, vacuum,cell,oneD, vTot%pw,vTot%mt,vTot%vacxy,vTot%vacz)
311 312 313 314 315 316 317
             !+odim
          END IF
          !           ----> m.t. boundaries
          nat = 1
          DO  n = 1,atoms%ntype
             CALL sphpts(xp,dimension%nspd,atoms%rmt(n),atoms%pos(1,nat))
             CALL checkdop(xp,dimension%nspd,n,nat,0,-1,1,.FALSE.,dimension,atoms,&
318
                  sphhar,stars,sym, vacuum,cell,oneD, vTot%pw,vTot%mt,vTot%vacxy,vTot%vacz)
319 320 321 322 323 324 325 326 327 328
             nat = nat + atoms%neq(n)
          ENDDO
          CALL timestop("checking")
       END IF
       !
       !========TOTAL==============================================
       !
       !      IF (l_xyav) THEN        ! write out xy-averaged potential & stop
       !        CALL xy_av_den(
       !     >                 n3d,k3d,nq3,nmzd,nmz,dvac,delz,
329
       !     >                 area,ig2,kv3,amat,vTot%pw,vTot%vacz(1,1,1))
330 331 332 333 334 335 336 337 338 339 340
       !         CALL juDFT_error("xy-averaged potential calculated",calledby="vgen")
       !      ENDIF

       IF (input%total) THEN
          CALL timestart("int_nv")

          !
          !      ---> AVERAGE COULOMB POTENTIAL ON THE SPHERE 
          !          FOR CALCULATING THE MADELUNG TERM in totale.f
          !           r=Rmt
          DO n=1,atoms%ntype
341
             atoms%vr0(n)=vTot%mt(atoms%jri(n),0,n,1)
342 343 344 345 346 347 348 349 350 351 352 353
          ENDDO
          !
          !     CALCULATE THE INTEGRAL OF n*Vcoulomb
          !
          WRITE (6,FMT=8020)
          WRITE (16,FMT=8020)
8020      FORMAT (/,10x,'density-coulomb potential integrals',/)
          !
          !       interstitial first
          !
          !       convolute ufft and pot: F(G) = \sum_(G') U(G - G') V(G')
          !
354
          CALL convol(stars, vpw_w, vTot%pw, stars%ufft)
355 356 357 358 359
          !
          IF (input%jspins.EQ.2) CALL CPP_BLAS_ccopy(stars%ng3,vpw_w(1,1),1,vpw_w(1,input%jspins),1)
          !
          results%te_vcoul = 0.0
          CALL int_nv(stars,vacuum,atoms,sphhar, cell,sym,input,oneD,&
360
               workDen%pw,vpw_w, workDen%vacxy,vTot%vacxy, workDen%vacz,vTot%vacz, workDen%mt,vTot%mt, results%te_vcoul)
361 362 363 364 365 366 367 368 369 370 371 372 373

          WRITE (6,FMT=8030) results%te_vcoul
          WRITE (16,FMT=8030) results%te_vcoul
8030      FORMAT (/,10x,'total density-coulomb potential integral :', t40,f20.10)

          CALL timestop("int_nv")

          INQUIRE(file='vdW_kernel_table',exist=l_vdw)
          IF (l_vdw) THEN

             CALL timestart("fleur_vdW")
             ! calculate vdW contribution to potential
             CALL fleur_vdW(mpi,atoms,sphhar,stars, input,dimension,&
374
                  cell,sym,oneD,vacuum, workDen%pw(:,1),workDen%mt(:,:,:,1), vpw_w(:,1),vTot%mt(:,:,:,:))
375 376 377 378 379 380 381 382 383 384 385
             CALL timestop("fleur_vdW")

          ENDIF

       END IF
       !
       !==========END TOTAL===================================================
       !
       !     ----> reload the density for calculating vxc (for spin-pol. case)
       !
       IF (input%jspins.EQ.2) THEN
386
          workDen = den
387 388
          vTot%mt(:,0:,:,2) = vTot%mt(:,0:,:,1)
          vTot%pw(:,2) = vTot%pw(:,1)
389
          IF (input%film) THEN
390 391
             vTot%vacxy(:,:,:,2) = vTot%vacxy(:,:,:,1)
             vTot%vacz(:,:,2)=vTot%vacz(:,:,1)
392 393 394 395 396
          END IF
       END IF
       IF (input%total) THEN
          DO js = 1,input%jspins
             ! to enable a GW calculation,
397
             vpw_w(1:stars%ng3,js)=vpw_w(1:stars%ng3,js)/stars%nstr(1:stars%ng3)     ! the PW-coulomb part is not
398 399
             ! used otherwise anyway.
          ENDDO
400 401 402 403 404 405
          vCoul%iter = vTot%iter
          vCoul%mt = vTot%mt
          vCoul%pw = vpw_w
          vCoul%vacz = vTot%vacz
          vCoul%vacxy = vTot%vacxy

406 407 408 409 410 411 412 413 414
          DO js = 1,input%jspins
             DO i = 1,stars%ng3
                vpw_w(i,js)=vpw_w(i,js)*stars%nstr(i)
             ENDDO
          ENDDO
       END IF
       IF (sliceplot%plpot) THEN
          OPEN (11,file='potcoul_pl',form='unformatted',status='unknown')
          CALL wrtdop(stars,vacuum,atoms,sphhar, input,sym,&
415
               11, vTot%iter,vTot%mt,vTot%pw,vTot%vacz,vTot%vacxy)
416 417
          CLOSE(11)
       END IF
418
    ENDIF !irank==0
419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435

       !     ******** exchange correlation potential******************
       !+ta
       !     rhmn: rho-min.
       !     ichsmrg: i-charge-small-region.
       !          0(not watched)
       !          1(in muffin-tin), 2(interstitial), 3(warped-vac),4(vacuum)

       ichsmrg=0
       rhmn=1.e+10
       !-ta
       excpw(:) = CMPLX(0.,0.)
       excz(:,:) = 0.0
       excxy(:,:,:) = CMPLX(0.,0.)
       excr(:,:,:) = 0.0

       !     ---> vacuum region
436
    IF (mpi%irank == 0) THEN
437 438 439 440
       IF (input%film) THEN

          CALL timestart("Vxc in vacuum")

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

444
          IF (.NOT.xcpot%is_gga()) THEN  ! LDA
445 446 447 448

             IF (.NOT.oneD%odi%d1) THEN

                CALL vvacxc(ifftd2,stars,vacuum,xcpot,input,noco,&
449
                     workDen%vacxy,workDen%vacz,workDen%cdomvxy,workDen%cdomvz, vTot%vacxy,vTot%vacz, excxy,excz)
450 451 452 453 454 455

             ELSE
                CALL judft_error("OneD broken")
                !           CALL vvacxc(&
                !     &                 stars,oneD%M,vacuum,odi%n2d,dimension,ifftd2,&
                !     &                 xcpot,input,odi%nq2,&
456
                !     &                 odi%nst2,workDen%vacxy,workDen%vacz,workDen%cdomvxy,workDen%cdomvz,noco,&
457
                !     &                 odi%kimax2%igf,odl%pgf,&
458
                !     &                 vTot%vacxy,vTot%vacz,&
459 460 461 462 463 464 465 466 467 468 469
                !     &                 excxy,excz)

             ENDIF

          ELSE      ! GGA

             IF (oneD%odi%d1) THEN
                CALL judft_error("OneD broken")

                CALL vvacxcg(ifftd2,stars,vacuum,noco,oneD,&
                     cell,xcpot,input,obsolete, ichsmrg,&
470
                     workDen%vacxy,workDen%vacz,workDen%cdomvxy,workDen%cdomvz, vTot%vacxy,vTot%vacz,rhmn, excxy,excz)
471 472 473 474

             ELSE
                CALL vvacxcg(ifftd2,stars,vacuum,noco,oneD,&
                     cell,xcpot,input,obsolete, ichsmrg,&
475
                     workDen%vacxy,workDen%vacz,workDen%cdomvxy,workDen%cdomvz, vTot%vacxy,vTot%vacz,rhmn, excxy,excz)
476 477 478 479 480 481 482 483 484 485 486

             END IF

          END IF
          CALL timestop("Vxc in vacuum")
          !+odim
       END IF
       !     ----------------------------------------
       !     ---> interstitial region
       CALL timestart("Vxc in interstitial")

487
       ifftd=27*stars%mx1*stars%mx2*stars%mx3
488

489
       IF ( (.NOT. obsolete%lwb) .OR. ( .not.xcpot%is_gga() ) ) THEN
490 491
          ! no White-Bird-trick

492
          ifftxc3d = stars%kxc1_fft*stars%kxc2_fft*stars%kxc3_fft
493

494
          IF ( .NOT.xcpot%is_gga() ) THEN
495 496
             ! LDA

497
             CALL visxc(ifftd,stars,noco,xcpot,input, workDen%pw,workDen%cdom,&
498
                  vTot%pw,vpw_w,vx%pw,vxpw_w, excpw)
499 500 501

          ELSE ! GGA

502
             CALL visxcg(ifftd,stars,sym, ifftxc3d, cell, workDen%pw,workDen%cdom, xcpot,input,&
503
                  obsolete,noco, rhmn,ichsmrg, vTot%pw,vpw_w,vx%pw,vxpw_w, excpw)
504 505 506 507 508 509

          END IF

       ELSE
          ! White-Bird-trick

Daniel Wortmann's avatar
Daniel Wortmann committed
510 511
          WRITE(6,'(a)') "W+B trick cancelled out. visxcwb uses at present common block cpgft3.",&
             "visxcwb needs to be reprogrammed according to visxcg.f"
512 513
          CALL juDFT_error("visxcwb",calledby ="vgen")
          !sb       CALL visxcwb(
514
          !sb  >                 workDen%pw,kimax,igfft,pgfft,ufft,
515 516 517
          !sb  >                 icorr,total,krla,
          !sb  >                 igrd,ndvgrd,idsprs,isprsv,
          !sb  >                 idsprsi,chng,sprsv,lwb,rhmn,ichsmrg,
518
          !sb  =                 vTot%pw,vpw_w,
519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549
          !sb  <                 excpw)

       END IF
       !
       ! --> on output vpw_w contains the warped effective potential and
       !               excpw the warped XC-energy density
       !

       !
       !     add IR EXX potential to vpw_w
       !
       INQUIRE(file='vpw_wexx',exist=exi)
       IF( exi ) THEN 
          WRITE(*,*) 'Read in vpw_wexx...'
          OPEN (351,file='vpw_wexx',form='formatted')
          DO js = 1,input%jspins
             DO i = 1,stars%ng3
                READ(351,'(2f30.15)') vpw_wexx(i,js) 
                vpw_w(i,js) = vpw_w(i,js) + vpw_wexx(i,js)*stars%nstr(i)
             END DO
          END DO
          CLOSE(351)
       END IF

       INQUIRE(file='vpw_exx',exist=exi)
       IF( exi ) THEN 
          WRITE(*,*) 'Read in vpw_exx...'
          OPEN (351,file='vpw_exx',form='formatted')
          DO js = 1,input%jspins
             DO i = 1,stars%ng3
                READ(351,'(2f30.15)') vpw_exx(i,js) 
550
                vTot%pw(i,js) = vTot%pw(i,js) + vpw_exx(i,js)
551 552 553 554 555 556 557 558 559 560 561 562
             END DO
          END DO
          CLOSE(351)
       END IF
       CALL timestop("Vxc in interstitial")

       !      ---> evaluate the interstitial average potential
       vbar(:) = vpw_w(1,:)*cell%omtil/cell%volint
       !-gu
       WRITE (6,FMT=8040) (vbar(js),js=1,input%jspins)
       WRITE (16,FMT=8040) (vbar(js),js=1,input%jspins)
8040   FORMAT (/,5x,'interstitial potential average (vbar) =',2f10.6)
563
    ENDIF !irank==0
564 565 566 567
       !
       !     ------------------------------------------
       !     ----> muffin tin spheres region

568
    IF (mpi%irank == 0) THEN
569
       CALL timestart ("Vxc in MT")
570
    END IF
571 572 573
#ifdef CPP_MPI
       CALL MPI_BCAST(atoms%vr0,atoms%ntype,MPI_DOUBLE_PRECISION,0,mpi%mpi_comm,ierr)
       CALL MPI_BCAST(input%efield%vslope,1,MPI_DOUBLE_COMPLEX,0,mpi%mpi_comm,ierr)
574
       CALL MPI_BCAST(workDen%mt,atoms%jmtd*(1+sphhar%nlhd)*atoms%ntype*dimension%jspd,MPI_DOUBLE_PRECISION,0,mpi%mpi_comm,ierr)
575
       CALL MPI_BCAST(vTot%mt,atoms%jmtd*(1+sphhar%nlhd)*atoms%ntype*dimension%jspd,MPI_DOUBLE_PRECISION,0,mpi%mpi_comm,ierr)
576 577 578
       CALL MPI_BCAST(rhmn,1,MPI_DOUBLE_PRECISION,0,mpi%mpi_comm,ierr)
       CALL MPI_BCAST(ichsmrg,1,MPI_INTEGER,0,mpi%mpi_comm,ierr)
#endif
579
       IF (xcpot%is_gga()) THEN
580
          CALL vmtxcg(dimension,mpi,sphhar,atoms, workDen%mt,xcpot,input,sym,&
581
               obsolete, vx%mt,vTot%mt,rhmn,ichsmrg, excr)
582
       ELSE
583
          CALL vmtxc(DIMENSION,sphhar,atoms, workDen%mt,xcpot,input,sym, vTot%mt, excr,vx%mt)
584 585 586 587 588 589 590
       ENDIF


       !
       ! add MT EXX potential to vr
       !

591
    IF (mpi%irank == 0) THEN
592 593
       INQUIRE(file='vr_exx',exist=exi)
       IF( exi ) THEN
Daniel Wortmann's avatar
Daniel Wortmann committed
594
          ALLOCATE( vr_exx(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,dimension%jspd) )
595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613
          OPEN (350,file='vr_exx',form='formatted')
          DO jsp = 1,dimension%jspd
             DO  itype = 1,atoms%ntype
                itypsym = atoms%ntypsy( SUM(atoms%neq(:itype-1)) + 1 )
                DO lh = 0,sphhar%nlh(itypsym)
                   l = sphhar%llh(lh,itypsym)
                   DO i = 1,atoms%jmtd
                      READ(350,'(f30.15)') vr_exx(i,lh,itype,jsp)
                   END DO

                   IF( l .EQ. 0 ) THEN
                      vr_exx(:,lh,itype,jsp) = vr_exx(:,lh,itype,jsp)*sfp_const/atoms%rmsh(:,itype)
                   END IF

                END DO
             END DO
          END DO
          CLOSE(350)

614
          vTot%mt = vTot%mt + vr_exx
615 616 617 618 619 620 621 622 623 624 625 626 627 628 629

       END IF
       CALL timestop ("Vxc in MT")
       !     ------------------------------------------
       !     ---> check continuity of total potential

       IF (input%vchk) THEN
          !           ----> vacuum boundaries
          IF (input%film .AND. .NOT.oneD%odi%d1) THEN
             npd = MIN(dimension%nspd,25)
             CALL points(xp,npd)
             DO ivac = 1,vacuum%nvac
                signum = 3. - 2.*ivac
                xp(3,:npd) = signum*cell%z1/cell%amat(3,3)
                CALL checkdop(xp,npd,0,0,ivac,1,1,.FALSE.,dimension,atoms, sphhar,stars,sym,&
630
                     vacuum,cell,oneD, vTot%pw,vTot%mt,vTot%vacxy,vTot%vacz)
631 632 633 634 635 636 637 638 639 640
             ENDDO ! ivac = 1,vacuum%nvac
          ELSEIF (oneD%odi%d1) THEN
             !-odim
             npd = MIN(dimension%nspd,25)
             CALL cylpts(xp,npd,cell%z1)
             !           DO j = 1,npd
             !              xp(1,j) = xp(1,j)/amat(1,1)
             !              xp(2,j) = xp(2,j)/amat(2,2)
             !           ENDDO
             CALL checkdop(xp,npd,0,0,vacuum%nvac,1,1,.FALSE.,dimension,atoms,&
641
                  sphhar,stars,sym, vacuum,cell,oneD, vTot%pw,vTot%mt,vTot%vacxy,vTot%vacz)
642 643 644 645 646 647 648
             !+odim
          END IF
          !           ----> m.t. boundaries
          nat = 1
          DO n = 1, atoms%ntype
             CALL sphpts(xp,dimension%nspd,atoms%rmt(n),atoms%pos(1,nat))
             CALL checkdop(xp,dimension%nspd,n,nat,0,-1,1,.FALSE.,dimension,&
649
                  atoms,sphhar,stars,sym, vacuum,cell,oneD, vTot%pw,vTot%mt,vTot%vacxy,vTot%vacz)
650 651 652 653 654
             nat = nat + atoms%neq(n)
          ENDDO ! n = 1, atoms%ntype
       END IF


655
       CALL pot_mod(atoms,sphhar,vacuum,stars, input, vTot%mt,vTot%vacxy,vTot%vacz,vTot%pw,vpw_w)
656 657 658 659 660
       !
       !============TOTAL======================================
       !
       IF (input%total) THEN

661
          IF (noco%l_noco) THEN ! load workDen%pw,workDen%vacz,workDen%vacxyxy from 'cdn'-file
662
             CALL readDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,CDN_ARCHIVE_TYPE_CDN_const,CDN_INPUT_DEN_const,&
663
                              0,fermiEnergyTemp,l_qfix,workDen)
664 665 666 667 668
          ENDIF
          !
          !     CALCULATE THE INTEGRAL OF n1*Veff1 + n2*Veff2
          !     Veff = Vcoulomb + Vxc
          !
669
          ALLOCATE( veffpw_w(stars%ng3,dimension%jspd) )
Daniel Wortmann's avatar
Daniel Wortmann committed
670
          ALLOCATE( veffr(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,dimension%jspd) )
671 672
          IF( xcpot%is_hybrid() ) THEN
             veffpw_w = vpw_w - xcpot%get_exchange_weight() * vxpw_w
673
             veffr    = vTot%mt    - xcpot%get_exchange_weight() * vx%mt
674 675
          ELSE
             veffpw_w = vpw_w
676
             veffr    = vTot%mt
677 678 679
          END IF

          !HF     kinetic energy correction for core states
680
          IF ( hybrid%l_addhf.AND.(xcpot%is_hybrid()  )  ) THEN
681 682 683 684
             ALLOCATE(rhoc(atoms%jmtd,atoms%ntype,DIMENSION%jspd), rhoc_vx(atoms%jmtd))
             ALLOCATE(tec(atoms%ntype,DIMENSION%jspd),qintc(atoms%ntype,DIMENSION%jspd))
             CALL readCoreDensity(input,atoms,dimension,rhoc,tec,qintc)
             DEALLOCATE(tec,qintc)
685 686 687 688 689 690 691 692 693 694
          END IF
          !HF

          results%te_veff = 0.0
          DO 370 js = 1,input%jspins
             WRITE (6,FMT=8050) js
             WRITE (16,FMT=8050) js
8050         FORMAT (/,10x,'density-effective potential integrals for spin ',i2,/)

             CALL int_nv(stars,vacuum,atoms,sphhar, cell,sym,input,oneD,&
695 696
                  workDen%pw(:,js),veffpw_w(:,js), workDen%vacxy(:,:,:,js),vTot%vacxy(:,:,:,js),&
                  workDen%vacz(:,:,js),vTot%vacz(:,:,js), workDen%mt(1,0,1,js),veffr(1,0,1,js), results%te_veff)
697 698

             !HF
699
             IF (hybrid%l_addhf.and.( xcpot%is_hybrid()  )  ) THEN
700 701 702
                nat = 1
                DO n = 1,atoms%ntype
                   DO j = 1, atoms%jri(n)
Daniel Wortmann's avatar
Daniel Wortmann committed
703
                      rhoc_vx(j) = rhoc(j,n,js) * vx%mt(j,0,n,js) / sfp_const
704 705
                   END DO
                   CALL intgr3(rhoc_vx,atoms%rmsh(1,n),atoms%dx(n),atoms%jri(n),dpdot)
706 707
                   IF( xcpot%is_hybrid() ) THEN
                      results%te_veff = results%te_veff + xcpot%get_exchange_weight() * dpdot*atoms%neq(n)
708 709 710 711 712 713 714 715 716
                   END IF
                   nat     = nat + atoms%neq(n)
                END DO
                !           Skip over parts in cdnc not used
             END IF
             !HF

370          CONTINUE

717
          
718 719 720 721 722 723 724 725 726 727 728 729
             DEALLOCATE( veffpw_w,veffr )
             WRITE (6,FMT=8060) results%te_veff
             WRITE (16,FMT=8060) results%te_veff
8060         FORMAT (/,10x,'total density-effective potential integral :', t40,f20.10)
             !
             !     CALCULATE THE INTEGRAL OF n*exc
             !
             !     ---> perform spin summation of charge densities
             !     ---> for the calculation of Exc
             IF (input%jspins.EQ.2) THEN
                nat = 1
                DO n = 1,atoms%ntype
730 731 732
                   workDen%mt(:atoms%jri(n),0:sphhar%nlh(atoms%ntypsy(nat)),n,1) = &
                      workDen%mt(:atoms%jri(n),0:sphhar%nlh(atoms%ntypsy(nat)),n,1) + &
                      workDen%mt(:atoms%jri(n),0:sphhar%nlh(atoms%ntypsy(nat)),n,input%jspins)
733 734 735

                   nat = nat + atoms%neq(n)
                ENDDO
736
                workDen%pw(:stars%ng3,1) = workDen%pw(:stars%ng3,1) + workDen%pw(:stars%ng3,input%jspins)
737
                IF (input%film) THEN
738 739 740 741 742
                   workDen%vacxy(:vacuum%nmzxy,:oneD%odi%nq2 - 1,:vacuum%nvac,1) = &
                        workDen%vacxy(:vacuum%nmzxy,:oneD%odi%nq2 - 1,:vacuum%nvac,1) + &
                        workDen%vacxy(:vacuum%nmzxy,:oneD%odi%nq2 - 1,:vacuum%nvac,input%jspins)
                   workDen%vacz(:vacuum%nmz,:vacuum%nvac,1) = workDen%vacz(:vacuum%nmz,:vacuum%nvac,1) +&
                        workDen%vacz(:vacuum%nmz,:vacuum%nvac,input%jspins)
743 744 745 746 747 748 749 750
                END IF
             END IF
             WRITE (6,FMT=8070)
             WRITE (16,FMT=8070)
8070         FORMAT (/,10x,'charge density-energy density integrals',/)

             results%te_exc = 0.0
             CALL int_nv(stars,vacuum,atoms,sphhar, cell,sym,input,oneD,&
751
                  workDen%pw(:,1),excpw(1), workDen%vacxy,excxy, workDen%vacz,excz, workDen%mt,excr, results%te_exc)
752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775
             WRITE (6,FMT=8080) results%te_exc
             WRITE (16,FMT=8080) results%te_exc

8080         FORMAT (/,10x,'total charge density-energy density integral :', t40,f20.10)

          END IF
          !
          !==========END TOTAL============================================
          !
          !           ---> store v(l=0) component as r*v(l=0)/sqrt(4pi)

          DO js = 1,input%jspins

             l_pottot = .FALSE.                  ! adds a B-field, if file
             IF (dimension%jspd.EQ.2) THEN                       !   mfee is present
                INQUIRE (file='mfee',exist=l_pottot)
                IF (l_pottot) THEN
                   OPEN (88,file='mfee',form='formatted',status='unknown')
                   REWIND 88
                   DO n=1,atoms%ntype
                      READ (88,*) i,mfie
                      WRITE (*,*) 'type,field:',i,mfie
                      IF (i/=n)  CALL juDFT_error("wrong types in mfee", calledby="vgen")
                      IF (js.EQ.1) THEN
776
                         vTot%mt(:atoms%jri(n),0,n,js) = vTot%mt(:atoms%jri(n),0,n,js) - mfie/2.
777
                      ELSE
778
                         vTot%mt(:atoms%jri(n),0,n,js) = vTot%mt(:atoms%jri(n),0,n,js) + mfie/2.
779 780 781 782 783 784 785 786 787
                      ENDIF
                   ENDDO
                   CLOSE (88)
                ENDIF
             ENDIF



             DO n = 1,atoms%ntype
788
                vTot%mt(:atoms%jri(n),0,n,js)  = atoms%rmsh(:atoms%jri(n),n)*vTot%mt(:atoms%jri(n),0,n,js)/sfp_const
Daniel Wortmann's avatar
Daniel Wortmann committed
789
                vx%mt(:atoms%jri(n),0,n,js) = atoms%rmsh(:atoms%jri(n),n)*vx%mt(:atoms%jri(n),0,n,js)/sfp_const
790 791 792 793 794 795 796 797 798 799 800 801
             ENDDO

          ENDDO     ! js =1,input%jspins


          IF ((.NOT.reap).OR.(noco%l_noco)) THEN
             IF (input%total) THEN
                OPEN (9,file='nrp',form='unformatted',status='unknown')
             ELSE
                OPEN (9,file='nrp',form='unformatted',position='append')
             ENDIF
             CALL wrtdop(stars,vacuum,atoms,sphhar, input,sym,&
802
                  9, vTot%iter,vTot%mt,vTot%pw,vTot%vacz,vTot%vacxy)
803 804 805 806 807 808 809 810
             CLOSE(9)
          ENDIF

          !     **************** reanalyze vpw *************************
          !                        call cpu_time(cp0)
          IF (input%total) THEN
             DO js=1,input%jspins
                DO i=1,stars%ng3
811
                   vTot%pw(i,js)=vpw_w(i,js)/stars%nstr(i)
812 813 814 815 816
                ENDDO
             ENDDO

             DO js=1,input%jspins
                DO i=1,stars%ng3
Daniel Wortmann's avatar
Daniel Wortmann committed
817
                   vx%pw(i,js)=vxpw_w(i,js)/stars%nstr(i)
818 819
                ENDDO
             ENDDO
820
             vx%iter = vTot%iter
821 822
          END IF

823 824 825 826 827 828 829 830 831 832 833
          IF (vacuum%nvac==1) THEN
             vTot%vacz(:,2,:)  = vTot%vacz(:,1,:)
             vCoul%vacz(:,2,:) = vCoul%vacz(:,1,:)
             IF (sym%invs) THEN
                vTot%vacxy(:,:,2,:)  = cmplx(vTot%vacxy(:,:,1,:))
                vCoul%vacxy(:,:,2,:) = cmplx(vCoul%vacxy(:,:,1,:))
             ELSE
                vTot%vacxy(:,:,2,:)  = vTot%vacxy(:,:,1,:)
                vCoul%vacxy(:,:,2,:) = vCoul%vacxy(:,:,1,:)
             ENDIF
          ENDIF
834 835
       ENDIF ! mpi%irank == 0

836
       ! broadcast potentials
837
#ifdef CPP_MPI
838 839
       CALL mpi_bc_potden(mpi,stars,sphhar,atoms,input,vacuum,oneD,noco,vTot)
       CALL mpi_bc_potden(mpi,stars,sphhar,atoms,input,vacuum,oneD,noco,vCoul)
840
       CALL mpi_bc_potden(mpi,stars,sphhar,atoms,input,vacuum,oneD,noco,vx)
841
#endif
842 843 844

     END SUBROUTINE vgen
   END MODULE m_vgen