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
    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
128
         vpw_exx(stars%ng3,dimension%jspd),vpw_wexx(stars%ng3,dimension%jspd),&
Daniel Wortmann's avatar
Daniel Wortmann committed
129
         excz(vacuum%nmzd,2),excr(atoms%jmtd,0:sphhar%nlhd,atoms%ntype),&
Daniel Wortmann's avatar
Daniel Wortmann committed
130
         vpw_w(stars%ng3,dimension%jspd),vxpw_w(stars%ng3,dimension%jspd),psq(stars%ng3) )
131
    vTot%iter = den%iter
132

133
    CALL workDen%init(stars,atoms,sphhar,vacuum,oneD,DIMENSION%jspd,.FALSE.,POTDEN_TYPE_DEN)
134
    IF (noco%l_noco) THEN
135 136
       ALLOCATE (workDen%cdom(stars%ng3),workDen%cdomvz(vacuum%nmzd,2))
       ALLOCATE (workDen%cdomvxy(vacuum%nmzxyd,oneD%odi%n2d-1,2))
137
    ELSE
138
       ALLOCATE (workDen%cdom(1),workDen%cdomvz(1,1),workDen%cdomvxy(1,1,1))
139 140 141
    END IF
    !

142
    workDen = den
143 144 145 146 147 148 149 150 151
    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
152 153
       END IF

154 155 156 157 158 159
       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
160 161
          workDen%mt(:,0:,:,1)=workDen%mt(:,0:,:,1)+workDen%mt(:,0:,:,2)
          workDen%pw(:,1)=workDen%pw(:,1)+workDen%pw(:,2)
162
          IF (input%film) THEN
163 164
             workDen%vacxy(:,:,:,1)=workDen%vacxy(:,:,:,1)+workDen%vacxy(:,:,:,2)
             workDen%vacz(:,:,1) = workDen%vacz(:,:,1) + workDen%vacz(:,:,2)
165 166 167 168 169 170 171 172 173 174 175 176
          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
177 178

#ifdef CPP_MPI
179
       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
180 181
#endif

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

       IF (l_xyav) THEN        ! write out xy-averaged density & stop
189
          CALL xy_av_den(stars,vacuum, cell,psq,workDen%vacz)
190 191 192 193 194 195 196 197 198
          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
199
          CALL od_vvac(stars,vacuum,cell, psq,workDen%vacz, vTot%vacz)
200 201 202 203 204 205 206

          !---> 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,&
207
               oneD, workDen%vacz,workDen%vacxy,psq,vTot%vacz,sym, vTot%vacxy,vTot%pw)
208 209 210 211 212 213 214
          CALL timestop("coulomb potential")

          !+odim
       ELSEIF (input%film .AND. .NOT.oneD%odi%d1) THEN
          !     ----> potential in the  vacuum  region
          !       
          CALL timestart("p vac") 
215 216
          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)
217

218
          CALL vvacxy(stars,vacuum,cell,sym,input, workDen%vacxy, vTot%vacxy, alphm)
219 220 221 222 223 224 225 226 227
          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
228
          ivfft =  3*stars%mx3 
229 230 231 232 233 234 235 236 237
          !         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,&
238
                     z,irec2, psq,vTot%vacxy,vTot%vacz,rhobar,sig1dh,vz1dh,alphm)
239 240 241 242 243 244 245 246 247
                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
248
             !        --> 1-d fourier transform and store the coefficients in vTot%pw( ,1)
249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268
             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)
269
                      vTot%pw(irec3,1) = vTot%pw(irec3,1) + xint/nzst1
270 271 272 273 274
                   END IF
                ENDIF
             ENDDO
          ENDDO
       ELSEIF (.NOT.input%film) THEN
275 276
          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))       
277 278 279 280 281 282 283 284 285
       END IF

       CALL timestop("p int")

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

    CALL timestart("p vmts")
286
    CALL vmts(mpi, stars,sphhar,atoms, sym,cell,oneD, vTot%pw,workDen%mt, vTot%mt)
287 288 289 290 291 292 293 294 295 296 297 298 299 300 301
    !     --------------------------------------------
    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,&
302
                     vacuum,cell,oneD, vTot%pw,vTot%mt,vTot%vacxy,vTot%vacz)
303 304 305 306 307 308 309 310 311 312
             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,&
313
                  sphhar,stars,sym, vacuum,cell,oneD, vTot%pw,vTot%mt,vTot%vacxy,vTot%vacz)
314 315 316 317 318 319 320
             !+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,&
321
                  sphhar,stars,sym, vacuum,cell,oneD, vTot%pw,vTot%mt,vTot%vacxy,vTot%vacz)
322 323 324 325 326 327 328 329 330 331
             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,
332
       !     >                 area,ig2,kv3,amat,vTot%pw,vTot%vacz(1,1,1))
333 334 335 336 337 338 339 340 341 342 343
       !         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
344
             atoms%vr0(n)=vTot%mt(atoms%jri(n),0,n,1)
345 346 347 348 349 350 351 352 353 354 355 356
          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')
          !
357
          CALL convol(stars, vpw_w, vTot%pw, stars%ufft)
358 359 360 361 362
          !
          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,&
363
               workDen%pw,vpw_w, workDen%vacxy,vTot%vacxy, workDen%vacz,vTot%vacz, workDen%mt,vTot%mt, results%te_vcoul)
364 365 366 367 368 369 370 371 372 373 374 375 376

          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,&
377
                  cell,sym,oneD,vacuum, workDen%pw(:,1),workDen%mt(:,:,:,1), vpw_w(:,1),vTot%mt(:,:,:,:))
378 379 380 381 382 383 384 385 386 387 388
             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
389
          workDen = den
390 391
          vTot%mt(:,0:,:,2) = vTot%mt(:,0:,:,1)
          vTot%pw(:,2) = vTot%pw(:,1)
392
          IF (input%film) THEN
393 394
             vTot%vacxy(:,:,:,2) = vTot%vacxy(:,:,:,1)
             vTot%vacz(:,:,2)=vTot%vacz(:,:,1)
395 396 397 398 399
          END IF
       END IF
       IF (input%total) THEN
          DO js = 1,input%jspins
             ! to enable a GW calculation,
400
             vpw_w(1:stars%ng3,js)=vpw_w(1:stars%ng3,js)/stars%nstr(1:stars%ng3)     ! the PW-coulomb part is not
401 402
             ! used otherwise anyway.
          ENDDO
403 404 405 406 407 408
          vCoul%iter = vTot%iter
          vCoul%mt = vTot%mt
          vCoul%pw = vpw_w
          vCoul%vacz = vTot%vacz
          vCoul%vacxy = vTot%vacxy

409 410 411 412 413 414 415 416 417
          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,&
418
               11, vTot%iter,vTot%mt,vTot%pw,vTot%vacz,vTot%vacxy)
419 420
          CLOSE(11)
       END IF
421
    ENDIF !irank==0
422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438

       !     ******** 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
439
    IF (mpi%irank == 0) THEN
440 441 442 443
       IF (input%film) THEN

          CALL timestart("Vxc in vacuum")

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

447
          IF (.NOT.xcpot%is_gga()) THEN  ! LDA
448 449 450 451

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

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

             ELSE
                CALL judft_error("OneD broken")
                !           CALL vvacxc(&
                !     &                 stars,oneD%M,vacuum,odi%n2d,dimension,ifftd2,&
                !     &                 xcpot,input,odi%nq2,&
459
                !     &                 odi%nst2,workDen%vacxy,workDen%vacz,workDen%cdomvxy,workDen%cdomvz,noco,&
460
                !     &                 odi%kimax2%igf,odl%pgf,&
461
                !     &                 vTot%vacxy,vTot%vacz,&
462 463 464 465 466 467 468 469 470 471 472
                !     &                 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,&
473
                     workDen%vacxy,workDen%vacz,workDen%cdomvxy,workDen%cdomvz, vTot%vacxy,vTot%vacz,rhmn, excxy,excz)
474 475 476 477

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

             END IF

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

490
       ifftd=27*stars%mx1*stars%mx2*stars%mx3
491

492
       IF ( (.NOT. obsolete%lwb) .OR. ( .not.xcpot%is_gga() ) ) THEN
493 494
          ! no White-Bird-trick

495
          ifftxc3d = stars%kxc1_fft*stars%kxc2_fft*stars%kxc3_fft
496

497
          IF ( .NOT.xcpot%is_gga() ) THEN
498 499
             ! LDA

500
             CALL visxc(ifftd,stars,noco,xcpot,input, workDen%pw,workDen%cdom,&
501
                  vTot%pw,vpw_w,vx%pw,vxpw_w, excpw)
502 503 504

          ELSE ! GGA

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

          END IF

       ELSE
          ! White-Bird-trick

Daniel Wortmann's avatar
Daniel Wortmann committed
513 514
          WRITE(6,'(a)') "W+B trick cancelled out. visxcwb uses at present common block cpgft3.",&
             "visxcwb needs to be reprogrammed according to visxcg.f"
515 516
          CALL juDFT_error("visxcwb",calledby ="vgen")
          !sb       CALL visxcwb(
517
          !sb  >                 workDen%pw,kimax,igfft,pgfft,ufft,
518 519 520
          !sb  >                 icorr,total,krla,
          !sb  >                 igrd,ndvgrd,idsprs,isprsv,
          !sb  >                 idsprsi,chng,sprsv,lwb,rhmn,ichsmrg,
521
          !sb  =                 vTot%pw,vpw_w,
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 550 551 552
          !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) 
553
                vTot%pw(i,js) = vTot%pw(i,js) + vpw_exx(i,js)
554 555 556 557 558 559 560 561 562 563 564 565
             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)
566
    ENDIF !irank==0
567 568 569 570
       !
       !     ------------------------------------------
       !     ----> muffin tin spheres region

571
    IF (mpi%irank == 0) THEN
572
       CALL timestart ("Vxc in MT")
573
    END IF
574 575 576
#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)
577
       CALL MPI_BCAST(workDen%mt,atoms%jmtd*(1+sphhar%nlhd)*atoms%ntype*dimension%jspd,MPI_DOUBLE_PRECISION,0,mpi%mpi_comm,ierr)
578
       CALL MPI_BCAST(vTot%mt,atoms%jmtd*(1+sphhar%nlhd)*atoms%ntype*dimension%jspd,MPI_DOUBLE_PRECISION,0,mpi%mpi_comm,ierr)
579 580 581
       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
582
       IF (xcpot%is_gga()) THEN
583
          CALL vmtxcg(dimension,mpi,sphhar,atoms, workDen%mt,xcpot,input,sym,&
584
               obsolete, vx%mt,vTot%mt,rhmn,ichsmrg, excr)
585
       ELSE
586
          CALL vmtxc(DIMENSION,sphhar,atoms, workDen%mt,xcpot,input,sym, vTot%mt, excr,vx%mt)
587 588 589 590 591 592 593
       ENDIF


       !
       ! add MT EXX potential to vr
       !

594
    IF (mpi%irank == 0) THEN
595 596
       INQUIRE(file='vr_exx',exist=exi)
       IF( exi ) THEN
Daniel Wortmann's avatar
Daniel Wortmann committed
597
          ALLOCATE( vr_exx(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,dimension%jspd) )
598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616
          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)

617
          vTot%mt = vTot%mt + vr_exx
618 619 620 621 622 623 624 625 626 627 628 629 630 631 632

       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,&
633
                     vacuum,cell,oneD, vTot%pw,vTot%mt,vTot%vacxy,vTot%vacz)
634 635 636 637 638 639 640 641 642 643
             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,&
644
                  sphhar,stars,sym, vacuum,cell,oneD, vTot%pw,vTot%mt,vTot%vacxy,vTot%vacz)
645 646 647 648 649 650 651
             !+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,&
652
                  atoms,sphhar,stars,sym, vacuum,cell,oneD, vTot%pw,vTot%mt,vTot%vacxy,vTot%vacz)
653 654 655 656 657
             nat = nat + atoms%neq(n)
          ENDDO ! n = 1, atoms%ntype
       END IF


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

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

          !HF     kinetic energy correction for core states
683
          IF ( hybrid%l_addhf.AND.(xcpot%is_hybrid()  )  ) THEN
684 685 686 687
             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)
688 689 690 691 692 693 694 695 696 697
          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,&
698 699
                  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)
700 701

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

370          CONTINUE

720
          
721 722 723 724 725 726 727 728 729 730 731 732
             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
733
                   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)
734 735 736

                   nat = nat + atoms%neq(n)
                ENDDO
737
                workDen%pw(:stars%ng3,1) = workDen%pw(:stars%ng3,1) + workDen%pw(:stars%ng3,input%jspins)
738
                IF (input%film) THEN
739 740 741 742 743
                   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)
744 745 746 747 748 749 750 751
                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,&
752
                  workDen%pw(:,1),excpw(1), workDen%vacxy,excxy, workDen%vacz,excz, workDen%mt,excr, results%te_exc)
753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776
             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
777
                         vTot%mt(:atoms%jri(n),0,n,js) = vTot%mt(:atoms%jri(n),0,n,js) - mfie/2.
778
                      ELSE
779
                         vTot%mt(:atoms%jri(n),0,n,js) = vTot%mt(:atoms%jri(n),0,n,js) + mfie/2.
780 781 782 783 784 785 786 787 788
                      ENDIF
                   ENDDO
                   CLOSE (88)
                ENDIF
             ENDIF



             DO n = 1,atoms%ntype
789
                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
790
                vx%mt(:atoms%jri(n),0,n,js) = atoms%rmsh(:atoms%jri(n),n)*vx%mt(:atoms%jri(n),0,n,js)/sfp_const
791 792 793 794 795 796 797 798 799 800 801 802
             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,&
803
                  9, vTot%iter,vTot%mt,vTot%pw,vTot%vacz,vTot%vacxy)
804 805 806 807 808 809 810 811
             CLOSE(9)
          ENDIF

          !     **************** reanalyze vpw *************************
          !                        call cpu_time(cp0)
          IF (input%total) THEN
             DO js=1,input%jspins
                DO i=1,stars%ng3
812
                   vTot%pw(i,js)=vpw_w(i,js)/stars%nstr(i)
813 814
                ENDDO
             ENDDO
815
             IF (vacuum%nvac==1) THEN
816
                vTot%vacz(:,2,:)=vTot%vacz(:,1,:)
817
                IF (sym%invs) THEN
818
                   vTot%vacxy(:,:,2,:)= cmplx(vTot%vacxy(:,:,1,:))
819
                ELSE
820
                   vTot%vacxy(:,:,2,:)=vTot%vacxy(:,:,1,:)
821 822
                ENDIF
             ENDIF
823 824 825

             DO js=1,input%jspins
                DO i=1,stars%ng3
Daniel Wortmann's avatar
Daniel Wortmann committed
826
                   vx%pw(i,js)=vxpw_w(i,js)/stars%nstr(i)
827 828
                ENDDO
             ENDDO
829
             vx%iter = vTot%iter
830 831 832 833
          END IF

       ENDIF ! mpi%irank == 0

834
       ! broadcast potentials
835
#ifdef CPP_MPI
836 837
       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)
838
       CALL mpi_bc_potden(mpi,stars,sphhar,atoms,input,vacuum,oneD,noco,vx)
839
#endif
840 841 842

     END SUBROUTINE vgen
   END MODULE m_vgen