visxcg.f90 13.9 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 9 10 11 12 13
MODULE m_visxcg
  USE m_juDFT
  !     ******************************************************
  !     subroutine generates the exchange-correlation potential
  !     in the interstitial region    c.l.fu
  !     including gradient corrections. t.a. 1996.
  !     ******************************************************
CONTAINS
14 15
  SUBROUTINE visxcg(ifftd,stars,sym, ifftxc3d, cell,den,&
       xcpot,input, obsolete,noco, vxc,vx,exc)
16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36

    !     ******************************************************
    !     instead of visxcor.f: the different exchange-correlation
    !     potentials defined through the key icorr are called through
    !     the driver subroutine vxcallg.f,for the energy density - excallg
    !     subroutines vectorized
    !     ** r.pentcheva 22.01.96
    !     *********************************************************
    !     in case of total = .true. calculates the ex-corr. energy
    !     density
    !     ** r.pentcheva 08.05.96
    !     ******************************************************************

    USE m_grdrsis
    USE m_prpxcfftmap
    USE m_mkgxyz3
    USE m_fft3d
    USE m_fft3dxc
    USE m_types
    IMPLICIT NONE

37
    CLASS(t_xcpot),INTENT(IN)     :: xcpot
38 39 40 41 42 43
    TYPE(t_obsolete),INTENT(IN)   :: obsolete
    TYPE(t_input),INTENT(IN)      :: input
    TYPE(t_noco),INTENT(IN)       :: noco
    TYPE(t_sym),INTENT(IN)        :: sym
    TYPE(t_stars),INTENT(IN)      :: stars
    TYPE(t_cell),INTENT(IN)       :: cell
Gregor Michalicek's avatar
Gregor Michalicek committed
44
    TYPE(t_potden),INTENT(IN)     :: den
45 46
    TYPE(t_potden),INTENT(INOUT)  :: vxc,vx,exc
    
47 48 49 50 51
    !     ..
    !     .. Scalar Arguments ..
    !     ..
    !     .. Array Arguments ..

52
    REAL rhmni ,d_15,sprsv
53 54 55 56 57 58 59 60 61 62 63
    !
    !----->  fft  information  for xc potential + energy
    !
    INTEGER, ALLOCATABLE :: igxc_fft(:)
    REAL,    ALLOCATABLE :: gxc_fft(:,:)
    !     ..
    !     .. Local Scalars ..
    INTEGER :: i ,k,js,nt,ifftxc3,idm,jdm,ndm,ig
    COMPLEX :: ci
    REAL    :: rhotot
    INTEGER :: ifftd,ifftxc3d
64
    TYPE(t_gradients)::grad
65 66 67 68 69 70 71 72 73

    !     ..
    !     .. Local Arrays ..
    COMPLEX, ALLOCATABLE :: fg3(:),cqpw(:,:), ph_wrk(:)
    REAL,    ALLOCATABLE :: bf3(:)
    REAL,    ALLOCATABLE :: rho(:,:),rhd1(:,:,:),rhd2(:,:,:)
    REAL,    ALLOCATABLE :: mx(:),my(:)
    REAL,    ALLOCATABLE :: magmom(:),dmagmom(:,:),ddmagmom(:,:,:) 
    ! 
74
    REAL, ALLOCATABLE :: v_x(:,:),v_xc(:,:),e_xc(:),vcon(:) 
75
  
76 77 78 79 80 81 82
    !     .. unused input (needed for other noco GGA-implementations) ..

    !ta+
    !.....------------------------------------------------------------------
    !------->          abbreviations
    !
    !     ph_wrk: work array containing phase * g_x,gy...... 
Gregor Michalicek's avatar
Gregor Michalicek committed
83
    !     den%pw: charge density stored as stars
84
    !     rho   : charge density stored in real space
85
    !     v_xc   : exchange-correlation potential in real space
86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114
    !     exc   : exchange-correlation energy density in real space
    !     kxc1d  : dimension of the charge density fft box in the pos. domain
    !     kxc2d  : defined in dimens.f program (subroutine apws).1,2,3 indic
    !     kxc3d  ; a_1, a_2, a_3 directions.
    !     kq(i) : i=1,2,3 actual length of the fft-box for which fft is done
    !     nstr  : number of members (arms) of reciprocal lattice (g) vector
    !             of each star.
    !     nxc3_fft: number of stars in the  charge density  fft-box
    !     ng3   : number of 3 dim. stars in the charge density sphere define
    !             by gmax
    !     kmxxc_fft: number of g-vectors forming the nxc3_fft stars in the
    !               charge density or xc-density sphere
    !     kimax : number of g-vectors forming the ng3 stars in the gmax-sphe
    !     ifftxc3d: elements (g-vectors) in the charge density  fft-box
    !     igfft : pointer from the g-sphere (stored as stars) to fft-grid
    !             and     from fft-grid to g-sphere (stored as stars)
    !     pgfft : contains the phases of the g-vectors of sph.
    !     isn   : isn = +1, fft transform for g-space to r-space
    !             isn = -1, vice versa
    !
    !-------------------------------------------------------------------
    !
    !---> set up pointer for backtransformation of from g-vector in
    !     positive domain of xc density fftbox into stars.
    !     also the x,y,z components of the g-vectors are set up to calculate
    !     derivatives.
    !     in principle this can also be done in main program once.
    !     it is done here to save memory.
    !
115 116
    ifftd=27*stars%mx1*stars%mx2*stars%mx3
    ifftxc3d = stars%kxc1_fft*stars%kxc2_fft*stars%kxc3_fft
117
    ALLOCATE ( igxc_fft(0:ifftxc3d-1),gxc_fft(0:ifftxc3d-1,3) )
118
    CALL prp_xcfft_map(stars,sym, cell, igxc_fft,gxc_fft)
119 120
    !
    ifftxc3=stars%kxc1_fft*stars%kxc2_fft*stars%kxc3_fft
121
   
122 123
    IF (stars%ng3.GT.stars%ng3) THEN
       WRITE(6,'(/'' stars%ng3.gt.stars%ng3. stars%ng3,stars%ng3='',2i6)') stars%ng3,stars%ng3
124 125 126 127 128 129 130 131 132
       CALL juDFT_error("ng3.gt.n3d",calledby="visxcg")
    ENDIF

    d_15=1.e-15
    !
    ci=CMPLX(0.,1.)
    !
    ! Allocate arrays
    ! ff
133 134
    ALLOCATE( bf3(0:ifftd-1),ph_wrk(0:ifftxc3d-1), rho(0:ifftxc3d-1,input%jspins),&
         rhd1(0:ifftxc3d-1,input%jspins,3), rhd2(0:ifftxc3d-1,input%jspins,6) )
135 136
    IF (noco%l_noco)  THEN
       ALLOCATE( mx(0:ifftxc3-1),my(0:ifftxc3-1),&
137
            magmom(0:ifftxc3-1), dmagmom(0:ifftxc3-1,3),ddmagmom(0:ifftxc3-1,3,3) )
138 139 140 141 142 143
    END IF


    !-->     transform charge density to real space

    DO js=1,input%jspins
144 145
       CALL fft3dxc(rho(0:,js),bf3, den%pw(:,js), stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft,&
            stars%nxc3_fft,stars%kmxxc_fft,+1, stars%igfft(0:,1),igxc_fft,stars%pgfft,stars%nstr)
146 147 148 149 150
    END DO

    IF (noco%l_noco) THEN  

       !       for off-diagonal parts the same
151 152
       CALL fft3dxc(mx,my, den%pw(:,3), stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft,&
            stars%nxc3_fft,stars%kmxxc_fft,+1, stars%igfft(0:,1),igxc_fft,stars%pgfft,stars%nstr)
153 154 155

       DO i=0,ifftxc3-1 
          rhotot= 0.5*( rho(i,1) + rho(i,2) )
156
          magmom(i)= SQRT(  (0.5*(rho(i,1)-rho(i,2)))**2 + mx(i)**2 + my(i)**2 )
157 158 159 160 161 162
          rho(i,1)= rhotot+magmom(i)
          rho(i,2)= rhotot-magmom(i)
       END DO

    ENDIF

163
    IF (.not.xcpot%is_gga()) GOTO 100  
164 165 166 167 168 169 170 171

    ! In collinear calculations all derivatives are calculated in g-spce,
    ! in non-collinear calculations the derivatives of |m| are calculated in real space. 

    !-->   for d(rho)/d(x,y,z) = rhd1(:,:,idm) (idm=1,2,3).
    !
    !         ph_wrk: exp(i*(g_x,g_y,g_z)*tau) * g_(x,y,z).

172
    ALLOCATE(cqpw(stars%ng3,input%jspins))
173 174 175

    DO js= 1,input%jspins
       DO i = 1,stars%ng3
Gregor Michalicek's avatar
Gregor Michalicek committed
176
          cqpw(i,js)= ci*den%pw(i,js)
177 178 179 180 181 182 183 184 185 186
       END DO
    END DO

    DO idm=1,3

       DO ig = 0 , stars%kmxxc_fft - 1
          ph_wrk(ig) = stars%pgfft(ig) * gxc_fft(ig,idm)
       END DO

       DO js=1,input%jspins
187 188
          CALL fft3dxc(rhd1(0:,js,idm),bf3, cqpw(:,js), stars%kxc1_fft,stars%kxc2_fft,&
               stars%kxc3_fft,stars%nxc3_fft,stars%kmxxc_fft,+1, stars%igfft(0:,1),igxc_fft,ph_wrk,stars%nstr)
189 190 191 192 193 194
       END DO

    END DO

    IF (noco%l_noco) THEN

195
       CALL grdrsis(magmom,cell,stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft,obsolete, dmagmom )
196 197 198 199 200 201 202 203 204 205 206

       DO i=0,ifftxc3-1
          DO idm=1,3
             rhotot= rhd1(i,1,idm)/2.+rhd1(i,2,idm)/2.
             rhd1(i,1,idm)= rhotot+dmagmom(i,idm) 
             rhd1(i,2,idm)= rhotot-dmagmom(i,idm) 
          END DO
       END DO

    END IF

207
   
208 209 210 211 212 213
    !-->   for dd(rho)/d(xx,xy,yy,zx,yz,zz) = rhd2(:,:,idm) (idm=1,2,3,4,5,6)
    !
    !         ph_wrk: exp(i*(g_x,g_y,g_z)*tau) * g_(x,y,z) * g_(x,y,z)

    DO i = 1,stars%ng3
       DO js=1,input%jspins 
Gregor Michalicek's avatar
Gregor Michalicek committed
214
          cqpw(i,js)= -den%pw(i,js)
215 216 217 218 219 220 221 222 223 224 225 226 227
       END DO
    END DO

    ndm = 0
    DO idm = 1,3
       DO jdm = 1,idm
          ndm = ndm + 1

          DO ig = 0 , stars%kmxxc_fft-1
             ph_wrk(ig) = stars%pgfft(ig)*gxc_fft(ig,idm)*gxc_fft(ig,jdm)
          ENDDO

          DO js=1,input%jspins
228 229
             CALL fft3dxc(rhd2(0:,js,ndm),bf3, cqpw(:,js), stars%kxc1_fft,stars%kxc2_fft,&
                  stars%kxc3_fft,stars%nxc3_fft,stars%kmxxc_fft,+1, stars%igfft(0:,1),igxc_fft,ph_wrk,stars%nstr)
230 231 232 233 234 235 236 237 238
          END DO
       END DO ! jdm 
    END DO   ! idm 

    DEALLOCATE(cqpw)

    IF (noco%l_noco) THEN

       DO idm = 1,3
239
          CALL grdrsis(dmagmom(0,idm),cell,stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft,obsolete, ddmagmom(0,1,idm) )
240 241 242 243 244 245 246 247 248
       END DO

       ndm= 0
       DO idm = 1,3
          DO jdm = 1,idm
             ndm = ndm + 1  

             DO i=0,ifftxc3-1
                rhotot= rhd2(i,1,ndm)/2.+rhd2(i,2,ndm)/2.
249 250
                rhd2(i,1,ndm)= rhotot + ( ddmagmom(i,jdm,idm) + ddmagmom(i,idm,jdm) )/2. 
                rhd2(i,2,ndm)= rhotot - ( ddmagmom(i,jdm,idm) + ddmagmom(i,idm,jdm) )/2. 
251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273
             END DO

          ENDDO !jdm
       ENDDO   !idm 

    END IF

100 CONTINUE


    DEALLOCATE ( ph_wrk )
    IF (noco%l_noco) THEN 
       DEALLOCATE(mx,my,magmom,dmagmom,ddmagmom) 
    END IF
    !
    DO js=1,input%jspins 
       DO i=0,ifftxc3-1
          rho(i,js)=MAX(rho(i,js),d_15)
       ENDDO
    END DO
    bf3=0.0
    ! allocate the other arrays 
    !
274 275
    CALL xcpot%alloc_gradients(ifftxc3d,input%jspins,grad)
 
276 277 278 279
    !
    !     calculate the quantities such as abs(grad(rho)),.. used in
    !     evaluating the gradient contributions to potential and energy.
    !
280 281
    CALL mkgxyz3 (ifftxc3d,input%jspins,ifftxc3,input%jspins,rho, rhd1(0,1,1),rhd1(0,1,2),rhd1(0,1,3),&
         rhd2(0,1,1),rhd2(0,1,3),rhd2(0,1,6), rhd2(0,1,5),rhd2(0,1,4),rhd2(0,1,2), grad)
282 283

    DEALLOCATE ( rhd1,rhd2 )
284 285
    ALLOCATE ( v_xc(0:ifftxc3d-1,input%jspins) )
    ALLOCATE ( v_x (0:ifftxc3d-1,input%jspins) )
286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301
    !
    !     calculate the exchange-correlation potential in  real space
    !
    nt=ifftxc3
    !
    !      rhmni: rho_minimum_interstitial.

    rhmni=10.e+10

    DO js=1,input%jspins
       DO i=0,ifftxc3-1
          rho(i,js)=MAX(rho(i,js),d_15)
          rhmni=MIN(rhmni,rho(i,js))
       ENDDO
    ENDDO

302 303
  
    IF (rhmni.LT.obsolete%chng) THEN
304
       WRITE(6,'(/'' rhmn.lt.obsolete%chng in visxc. rhmn,obsolete%chng='',2d9.2)') rhmni,obsolete%chng
305 306
       !          CALL juDFT_error("visxcg: rhmn.lt.chng",calledby="visxcg")
    ENDIF
307
    CALL xcpot%get_vxc(input%jspins,rho,v_xc,v_x,grad)
308 309
    !
    !----> back fft to g space
310
    !----> perform back  fft transform: v_xc(r) --> vxc(star)
311
    !
312
    ALLOCATE(fg3(stars%ng3))
313 314 315

    DO js = 1,input%jspins
       bf3=0.0
316 317
       CALL fft3dxc(v_xc(0:,js),bf3, fg3, stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft,&
            stars%nxc3_fft,stars%kmxxc_fft,-1, stars%igfft(0:,1),igxc_fft,stars%pgfft,stars%nstr)
318 319
       !
       DO k = 1,stars%nxc3_fft
320
          vxc%pw(k,js) = vxc%pw(k,js) + fg3(k)
321 322 323 324 325 326 327
       ENDDO

       !
       !====>  INCLUDING TOTAL ENERGY
       !
       IF (input%total) THEN
          !
328
          !----> Perform fft transform: v_xc(star) --> vxc(r) 
329 330 331 332
          !     !Use large fft mesh for convolution
          !
          fg3(stars%nxc3_fft+1:)=0.0
          ALLOCATE ( vcon(0:ifftd-1) )
333
          CALL fft3d(vcon(0),bf3, fg3, stars,+1)
334 335 336 337 338 339 340
          !
          !----> Convolute with step function
          !
          DO i=0,ifftd-1
             vcon(i)=stars%ufft(i)*vcon(i)
          ENDDO
          bf3=0.0
341
          CALL fft3d(vcon(0),bf3, fg3, stars,-1,.FALSE.)
342 343 344 345 346
          DEALLOCATE ( vcon )
          !
          !----> add to warped coulomb potential
          !
          DO k = 1,stars%ng3
347
             vxc%pw_w(k,js) = vxc%pw_w(k,js) + fg3(k)
348 349 350 351 352
          ENDDO

       ENDIF

       bf3=0.0
353 354
       CALL fft3dxc(v_x(0:,js),bf3, fg3, stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft, &
            stars%nxc3_fft,stars%kmxxc_fft,-1, stars%igfft(0:,1),igxc_fft,stars%pgfft,stars%nstr)
355 356
       !
       DO k = 1,stars%nxc3_fft
357
          vx%pw(k,js) = vx%pw(k,js) + fg3(k)
358 359 360 361 362 363 364
       ENDDO

       !
       !====>   INCLUDING TOTAL ENERGY
       !
       IF (input%total) THEN
          !
365
          !---->  Perform fft transform: v_xc(star) --> vxc(r) 
366 367 368 369 370
          !       !Use large fft mesh for convolution
          !

          fg3(stars%nxc3_fft+1:)=0.0
          ALLOCATE ( vcon(0:ifftd-1) )
371
          CALL fft3d(vcon(0),bf3, fg3, stars,+1)
372 373 374 375 376 377 378
          !
          !----> Convolute with step function
          !
          DO i=0,ifftd-1
             vcon(i)=stars%ufft(i)*vcon(i)
          ENDDO
          bf3=0.0
379
          CALL fft3d(vcon(0),bf3, fg3, stars,-1,.FALSE.)
380 381 382 383 384
          DEALLOCATE ( vcon )
          !
          !----> add to warped exchange-potential
          !
          DO k = 1,stars%ng3
385
             vx%pw_w(k,js) = vx%pw_w(k,js) + fg3(k)
386 387 388 389 390
          ENDDO

       ENDIF

    ENDDO
391
    DEALLOCATE ( v_x,v_xc )
392 393 394
    !
    !     calculate the ex.-cor energy density in real space
    !
395
    IF (ALLOCATED(exc%pw_w)) THEN
396
       ALLOCATE ( e_xc(0:ifftxc3d-1) )
397
       CALL xcpot%get_exc(input%jspins,rho,e_xc,grad)
398 399 400 401
       !
       !---->   perform back  fft transform: exc(r) --> exc(star)
       !
       bf3=0.0
402 403
       CALL fft3dxc(e_xc,bf3, fg3, stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft,&
            stars%nxc3_fft,stars%kmxxc_fft,-1, stars%igfft(0:,1),igxc_fft,stars%pgfft,stars%nstr)
404
       DEALLOCATE ( e_xc )
405 406 407 408 409 410 411
       !
       !---->   Perform fft transform: exc(star) --> exc(r) 
       !        !Use large fft mesh for convolution
       !
       fg3(stars%nxc3_fft+1:)=0.0
       bf3=0.0
       ALLOCATE ( vcon(0:ifftd-1) )
412
       CALL fft3d(vcon,bf3,fg3, stars,+1)
413 414 415 416 417 418 419 420

       DO i=0,ifftd-1
          vcon(i)=stars%ufft(i)*vcon(i)
       ENDDO
       !
       !         ---> back fft to g space
       !
       bf3=0.0
421
       CALL fft3d(vcon,bf3,exc%pw_w(:,1), stars,-1,.FALSE.)
422 423 424 425 426 427
       DEALLOCATE ( vcon )
       !
    ENDIF

    DEALLOCATE(fg3)
    DEALLOCATE ( bf3,rho,igxc_fft,gxc_fft )
428
  
429 430 431 432 433



  END SUBROUTINE visxcg
END MODULE m_visxcg