pw_tofrom_grid.F90 15.7 KB
Newer Older
1 2 3 4 5 6
!--------------------------------------------------------------------------------
! 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.
!--------------------------------------------------------------------------------
MODULE m_pw_tofrom_grid
7 8 9
   USE m_types
   PRIVATE
   REAL,PARAMETER:: d_15=1.e-15
10

11 12 13 14
   INTEGER :: ifftd,ifftxc3
   !----->  fft  information  for xc potential + energy
   INTEGER, ALLOCATABLE :: igxc_fft(:)
   REAL,    ALLOCATABLE :: gxc_fft(:,:) !gxc_fft(ig,idm)
15
  
16
   PUBLIC :: init_pw_grid, pw_to_grid, pw_from_grid, finish_pw_grid
17
CONTAINS
18
  SUBROUTINE init_pw_grid(dograds,stars,sym,cell)
19 20
    USE m_prpxcfftmap
    USE m_types
21
    IMPLICIT NONE
22
    LOGICAL,INTENT(IN)            :: dograds
23 24 25 26
    TYPE(t_stars),INTENT(IN)      :: stars
    TYPE(t_sym),INTENT(IN)        :: sym
    TYPE(t_cell),INTENT(IN)       :: cell
    
27 28 29 30 31 32
      !---> 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.
33
      
34
    ifftd=27*stars%mx1*stars%mx2*stars%mx3
Matthias Redies's avatar
Matthias Redies committed
35
    ifftxc3  = stars%kxc1_fft*stars%kxc2_fft*stars%kxc3_fft
36
    IF (dograds) THEN
37 38
       CALL prp_xcfft_map(stars,sym, cell, igxc_fft,gxc_fft)
    ENDIF
39 40 41
       
  END SUBROUTINE init_pw_grid
  
42
  SUBROUTINE pw_to_grid(dograds,jspins,l_noco,stars,cell,den_pw,grad,xcpot,rho)
43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67
    !.....------------------------------------------------------------------
    !------->          abbreviations
    !
    !     ph_wrk: work array containing phase * g_x,gy...... 
    !     den%pw: charge density stored as stars
    !     rho   : charge density stored in real space
    !     v_xc   : exchange-correlation potential in real space
    !     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
    !     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
68 69
    !     l_rdm : if true, calculate noco gradients in local frame of density matrix
    !             improves convergence in case of GGA + noco (best with large G_max_xc)
70 71 72
    !
    !-------------------------------------------------------------------
    USE m_grdrsis
73
    USE m_polangle
74 75
    USE m_mkgxyz3
    USE m_fft3dxc
76
    USE m_fft3d
77
    USE m_types
78
    USE m_constants
79
    IMPLICIT NONE
80
    
81 82 83 84 85 86 87 88 89
    LOGICAL,INTENT(IN)                    :: dograds
    INTEGER,INTENT(IN)                    :: jspins
    LOGICAL,INTENT(IN)                    :: l_noco
    TYPE(t_stars),INTENT(IN)              :: stars
    TYPE(t_cell),INTENT(IN)               :: cell
    COMPLEX,INTENT(IN)                    :: den_pw(:,:)
    TYPE(t_gradients),INTENT(OUT)         :: grad
    CLASS(t_xcpot), INTENT(IN),OPTIONAL   :: xcpot
    REAL,ALLOCATABLE,INTENT(OUT),OPTIONAL :: rho(:,:)
90
  
91

92 93 94
    INTEGER      :: js,i,idm,ig,ndm,jdm,j
    REAL         :: rhotot,mmx,mmy,mmz,theta,phi
    COMPLEX      :: ci,rho21
95 96 97 98 99 100
    !     .. Local Arrays ..
    COMPLEX, ALLOCATABLE :: cqpw(:,:),ph_wrk(:)
    REAL,    ALLOCATABLE :: bf3(:)
    REAL,    ALLOCATABLE :: rhd1(:,:,:),rhd2(:,:,:)
    REAL,    ALLOCATABLE :: mx(:),my(:)
    REAL,    ALLOCATABLE :: magmom(:),dmagmom(:,:),ddmagmom(:,:,:) 
101 102 103 104 105 106 107
    REAL,    ALLOCATABLE :: rhodiag(:,:),der(:,:,:),dder(:,:,:,:)
    REAL,    ALLOCATABLE :: sinsqu(:),cossqu(:),sincos(:),rhdd(:,:,:,:)
    COMPLEX, ALLOCATABLE :: exi(:)
 
    LOGICAL, PARAMETER :: l_rdm=.true. 
    ci=cmplx(0.,1.)

108
    ! Allocate arrays
109
    ALLOCATE( bf3(0:ifftd-1))
110
    IF (dograds) THEN
Matthias Redies's avatar
Matthias Redies committed
111 112 113
       IF (PRESENT(rho)) ALLOCATE(rho(0:ifftxc3-1,jspins))
       ALLOCATE( ph_wrk(0:ifftxc3-1),rhd1(0:ifftxc3-1,jspins,3))
       ALLOCATE( rhd2(0:ifftxc3-1,jspins,6) )
114
     ELSE
115
        IF (PRESENT(rho)) ALLOCATE(rho(0:ifftd-1,jspins))
116
     ENDIF
117
    IF (l_noco)  THEN
118
       IF (dograds) THEN
Daniel Wortmann's avatar
Daniel Wortmann committed
119
          ALLOCATE( mx(0:ifftxc3-1),my(0:ifftxc3-1),magmom(0:ifftxc3-1))
120
          IF (l_rdm) THEN
121 122
           ALLOCATE( rhodiag(0:ifftxc3-1,jspins),der(0:ifftxc3-1,3,4),dder(0:ifftxc3-1,3,3,4),rhdd(0:ifftxc3-1,2,3,3) )
           ALLOCATE( sinsqu(0:ifftxc3-1),cossqu(0:ifftxc3-1),sincos(0:ifftxc3-1),exi(0:ifftxc3-1) )
123 124 125
          ELSE
            ALLOCATE( dmagmom(0:ifftxc3-1,3),ddmagmom(0:ifftxc3-1,3,3) )
          ENDIF
Daniel Wortmann's avatar
Daniel Wortmann committed
126 127 128
       ELSE
          ALLOCATE( mx(0:ifftd-1),my(0:ifftd-1),magmom(0:ifftd-1))
       ENDIF
129
    END IF
Matthias Redies's avatar
Matthias Redies committed
130

131 132 133
    IF (PRESENT(rho)) THEN
    !Put den_pw on grid and store into rho(:,1:2)
       DO js=1,jspins
134
          IF (dograds) THEN
135 136 137 138 139 140
             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)
          ELSE
             CALL fft3d(rho(0,js),bf3, den_pw(:,js), stars,+1)
          ENDIF
       END DO
141

142 143
       IF (l_noco) THEN  
          !  Get mx,my on real space grid and recalculate rho and magmom
144
          IF (dograds) THEN
145 146 147 148 149
             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)
          ELSE
             CALL fft3d(mx,my, den_pw(:,3), stars,+1)
          ENDIF
Daniel Wortmann's avatar
Daniel Wortmann committed
150
          DO i=0,MIN(SIZE(rho,1),size(mx))-1 
151 152
             rhotot= 0.5*( rho(i,1) + rho(i,2) )
             magmom(i)= SQRT(  (0.5*(rho(i,1)-rho(i,2)))**2 + mx(i)**2 + my(i)**2 )
153
             IF (l_rdm.AND.dograds) THEN
154 155 156
                rhodiag(i,1) = rho(i,1)
                rhodiag(i,2) = rho(i,2)
             ENDIF
157 158
             rho(i,1)= rhotot+magmom(i)
             rho(i,2)= rhotot-magmom(i)
159
             IF (l_rdm.AND.dograds) THEN              ! prepare rotation matrix
160 161 162 163 164 165 166 167 168
                mmx = 2.0 * mx(i)
                mmy = 2.0 * my(i)
                mmz = rhodiag(i,1) - rhodiag(i,2)
                CALL pol_angle(mmx,mmy,mmz,theta,phi)
                cossqu(i) = cos(0.5*theta)**2
                sinsqu(i) = sin(0.5*theta)**2
                sincos(i) = 2.0*sin(0.5*theta)*cos(0.5*theta)
                exi(i)    = exp(-ci*phi)
             ENDIF
169
          END DO
170
       ENDIF
171
    ENDIF
172
    IF (dograds) THEN
173 174 175 176 177 178 179 180

    ! 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).

181
       ALLOCATE(cqpw(stars%ng3,jspins))
182

183
       cqpw(:,:)= ImagUnit*den_pw(:,:jspins)
184 185 186 187 188 189
   
       DO idm=1,3
          DO ig = 0 , stars%kmxxc_fft - 1
             ph_wrk(ig) = stars%pgfft(ig) * gxc_fft(ig,idm)
          END DO

190
          DO js=1,jspins
191 192 193 194 195
             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)
          END DO
       END DO

196
       IF (l_noco) THEN
197

198 199 200 201 202 203
          IF (l_rdm) THEN

             CALL grdrsis( rhodiag(0,1),cell,stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft,der(0,1,1) )
             CALL grdrsis( rhodiag(0,2),cell,stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft,der(0,1,2) )
             CALL grdrsis( mx(0),cell,stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft,der(0,1,3) )
             CALL grdrsis( my(0),cell,stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft,der(0,1,4) )
204

205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226
             DO i=0,ifftxc3-1 ! project on magnetization axis
                DO idm=1,3
                   rho21 = der(i,idm,3) + ci * der(i,idm,4)
                   rhd1(i,1,idm) = cossqu(i) * der(i,idm,1) + &
                                   sincos(i) * real( exi(i)*rho21 ) + &
                                   sinsqu(i) * der(i,idm,2)
                   rhd1(i,2,idm) = sinsqu(i) * der(i,idm,1) - &
                                   sincos(i) * real( exi(i)*rho21 ) + &
                                   cossqu(i) * der(i,idm,2)
                ENDDO
             ENDDO

          ELSE

             CALL grdrsis(magmom,cell,stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft,dmagmom )

             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
227
             END DO
228 229 230

          END IF

231 232 233 234 235 236
       END IF

       !-->   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)

237
       cqpw(:,:)= -den_pw(:,:jspins)
238 239 240 241 242 243 244 245 246
   
       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
             
247
             DO js=1,jspins
248 249 250 251 252 253 254 255
                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)
             END DO
          END DO ! jdm 
       END DO   ! idm 

       DEALLOCATE(cqpw)

256
       IF (l_noco) THEN
257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304

          IF (l_rdm) THEN
             DO idm = 1,3
               CALL grdrsis( der(0,idm,1),cell,stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft, dder(0,1,idm,1) )
               CALL grdrsis( der(0,idm,2),cell,stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft, dder(0,1,idm,2) )
               CALL grdrsis( der(0,idm,3),cell,stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft, dder(0,1,idm,3) )
               CALL grdrsis( der(0,idm,4),cell,stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft, dder(0,1,idm,4) )
               DO i=0,ifftxc3-1 ! project on magnetization axis
                 DO jdm=1,3
                   rho21 = dder(i,jdm,idm,3) + ci * dder(i,jdm,idm,4)
                   rhdd(i,1,jdm,idm) = cossqu(i) * dder(i,jdm,idm,1) + &
                                   sincos(i) * real( exi(i)*rho21 ) + &
                                   sinsqu(i) * dder(i,jdm,idm,2)
                   rhdd(i,2,jdm,idm) = sinsqu(i) * dder(i,jdm,idm,1) - &
                                   sincos(i) * real( exi(i)*rho21 ) + &
                                   cossqu(i) * dder(i,jdm,idm,2)
                 ENDDO
               ENDDO
             ENDDO
             DO j=1,2
               DO i=0,ifftxc3-1
                 rhd2(i,j,1) = rhdd(i,j,1,1)
                 rhd2(i,j,2) = 0.5*(rhdd(i,j,1,2)+rhdd(i,j,2,1)) ! xy - averaging should be unneccessary
                 rhd2(i,j,3) = rhdd(i,j,2,2)
                 rhd2(i,j,4) = 0.5*(rhdd(i,j,1,3)+rhdd(i,j,3,1)) ! zx
                 rhd2(i,j,5) = 0.5*(rhdd(i,j,2,3)+rhdd(i,j,3,2)) ! yz
                 rhd2(i,j,6) = rhdd(i,j,3,3)
               ENDDO
             ENDDO
             DEALLOCATE (rhodiag,der,rhdd,dder,sinsqu,cossqu,sincos,exi)
          ELSE

             DO idm = 1,3
                CALL grdrsis(dmagmom(0,idm),cell,stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft,ddmagmom(0,1,idm) )
             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.
                      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. 
                   END DO
                ENDDO !jdm
             ENDDO   !idm 
             DEALLOCATE(dmagmom,ddmagmom)
          END IF
305
       END IF
306

307 308 309
       IF (PRESENT(xcpot)) THEN
          CALL xcpot%alloc_gradients(ifftxc3,jspins,grad)
       END IF
310 311

       ! 
312 313
       !     calculate the quantities such as abs(grad(rho)),.. used in
       !     evaluating the gradient contributions to potential and energy.
314
       ! 
315 316 317 318 319
       IF (PRESENT(rho)) THEN
          CALL mkgxyz3 (rho,rhd1(0:,:,1),rhd1(0:,:,2),rhd1(0:,:,3),&
               rhd2(0:,:,1),rhd2(0:,:,3),rhd2(0:,:,6), rhd2(0:,:,5),rhd2(0:,:,4),rhd2(0:,:,2),grad)
       ELSE
          !Dummy rho (only possible if grad is used for libxc mode)
320 321
          !CALL mkgxyz3 (RESHAPE((/0.0/),(/1,1/)),rhd1(0:,:,1),rhd1(0:,:,2),rhd1(0:,:,3),&
          !     rhd2(0:,:,1),rhd2(0:,:,3),rhd2(0:,:,6), rhd2(0:,:,5),rhd2(0:,:,4),rhd2(0:,:,2),grad)
322 323 324 325 326

          IF (dograds.and.(.not.PRESENT(xcpot))) THEN
             ALLOCATE(grad%gr(3,ifftxc3,1))
          END IF

327
          CALL mkgxyz3 (0*rhd1(0:,:,1),rhd1(0:,:,1),rhd1(0:,:,2),rhd1(0:,:,3),&
328 329
               rhd2(0:,:,1),rhd2(0:,:,3),rhd2(0:,:,6), rhd2(0:,:,5),rhd2(0:,:,4),rhd2(0:,:,2),grad)
       END IF
330 331
       
    ENDIF
332 333 334
    IF (PRESENT(rho)) THEN
       WHERE(ABS(rho) < d_15) rho = d_15
    ENDIF
335 336 337 338
   
  END SUBROUTINE pw_to_grid


339
  SUBROUTINE pw_from_grid(dograds,stars,l_pw_w,v_in,v_out_pw,v_out_pw_w)
340 341 342 343
    USE m_fft3d
    USE m_fft3dxc
    USE m_types
    IMPLICIT NONE
344
    LOGICAL,INTENT(IN)            :: dograds
345
    TYPE(t_stars),INTENT(IN)      :: stars
346
    REAL,INTENT(INOUT)            :: v_in(0:,:)
347
    LOGICAL,INTENT(in)            :: l_pw_w
348
    COMPLEX,INTENT(INOUT)         :: v_out_pw(:,:)
Matthias Redies's avatar
Matthias Redies committed
349
    COMPLEX,INTENT(INOUT),OPTIONAL:: v_out_pw_w(:,:)
350 351 352 353 354 355
    
    
    INTEGER              :: js,k,i
    REAL,ALLOCATABLE     :: bf3(:),vcon(:) 
    COMPLEX, ALLOCATABLE :: fg3(:)
    ALLOCATE( bf3(0:ifftd-1),fg3(stars%ng3))
356
    ALLOCATE ( vcon(0:ifftd-1) )
357 358
    DO js = 1,SIZE(v_in,2)
       bf3=0.0
359
       IF (dograds) THEN
360 361 362 363 364 365
          CALL fft3dxc(v_in(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)
       ELSE
          vcon(0:)=v_in(0:,js)
          CALL fft3d(v_in(0:,js),bf3, fg3, stars,-1)
       ENDIF
366
       DO k = 1,MERGE(stars%nxc3_fft,stars%ng3,dograds)
367
          v_out_pw(k,js) = v_out_pw(k,js) + fg3(k)
368 369 370
       ENDDO

       IF (l_pw_w) THEN
371
          IF (dograds) THEN
372 373 374 375 376
             !----> Perform fft transform: v_xc(star) --> vxc(r) 
             !     !Use large fft mesh for convolution
             fg3(stars%nxc3_fft+1:)=0.0
             CALL fft3d(vcon(0),bf3, fg3, stars,+1)
          ENDIF
377 378 379 380 381 382 383
          !
          !----> Convolute with step function
          !
          DO i=0,ifftd-1
             vcon(i)=stars%ufft(i)*vcon(i)
          ENDDO
          bf3=0.0
384 385
          CALL fft3d(vcon(0),bf3, fg3, stars,-1)
          fg3=fg3*stars%nstr
386 387 388
          !
          !----> add to warped coulomb potential
          !
389 390 391 392 393
          IF (PRESENT(v_out_pw_w)) THEN
             DO k = 1,stars%ng3
                v_out_pw_w(k,js) = v_out_pw_w(k,js) + fg3(k)
             ENDDO
          END IF
394 395 396 397 398 399
       ENDIF
    END DO
  END SUBROUTINE pw_from_grid
    
  SUBROUTINE finish_pw_grid()
    IMPLICIT NONE
400
    IF (ALLOCATED(igxc_fft)) DEALLOCATE(igxc_fft,gxc_fft)
401 402 403
  END SUBROUTINE finish_pw_grid

END MODULE m_pw_tofrom_grid