pw_tofrom_grid.F90 15.8 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 119 120 121
       IF (l_rdm) THEN
          ALLOCATE (rhodiag(0:SIZE(rho,1)-1,jspins))
          ALLOCATE (sinsqu(0:SIZE(rho,1)-1),cossqu(0:SIZE(rho,1)-1),sincos(0:SIZE(rho,1)-1),exi(0:SIZE(rho,1)-1))
       END IF
122
       IF (dograds) THEN
Daniel Wortmann's avatar
Daniel Wortmann committed
123
          ALLOCATE( mx(0:ifftxc3-1),my(0:ifftxc3-1),magmom(0:ifftxc3-1))
124
          IF (l_rdm) THEN
125 126
           ALLOCATE( der(0:ifftxc3-1,3,4),dder(0:ifftxc3-1,3,3,4),rhdd(0:ifftxc3-1,2,3,3) )
           
127 128 129
          ELSE
            ALLOCATE( dmagmom(0:ifftxc3-1,3),ddmagmom(0:ifftxc3-1,3,3) )
          ENDIF
Daniel Wortmann's avatar
Daniel Wortmann committed
130 131 132
       ELSE
          ALLOCATE( mx(0:ifftd-1),my(0:ifftd-1),magmom(0:ifftd-1))
       ENDIF
133
    END IF
Matthias Redies's avatar
Matthias Redies committed
134

135 136 137
    IF (PRESENT(rho)) THEN
    !Put den_pw on grid and store into rho(:,1:2)
       DO js=1,jspins
138
          IF (dograds) THEN
139 140 141 142 143 144
             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
145

146 147
       IF (l_noco) THEN  
          !  Get mx,my on real space grid and recalculate rho and magmom
148
          IF (dograds) THEN
149 150 151 152 153
             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
154
          DO i=0,MIN(SIZE(rho,1),size(mx))-1 
155 156
             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 )
157 158 159 160
             IF (l_rdm) THEN
                rhodiag(i,1) = rho(i,1)
                rhodiag(i,2) = rho(i,2)
             ENDIF
161 162
             rho(i,1)= rhotot+magmom(i)
             rho(i,2)= rhotot-magmom(i)
163 164 165 166 167 168 169 170 171 172
             IF (l_rdm) THEN              ! prepare rotation matrix
                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
173
          END DO
174
       ENDIF
175
    ENDIF
176
    IF (dograds) THEN
177 178 179 180 181 182 183 184

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

185
       ALLOCATE(cqpw(stars%ng3,jspins))
186

187
       cqpw(:,:)= ImagUnit*den_pw(:,:jspins)
188 189 190 191 192 193
   
       DO idm=1,3
          DO ig = 0 , stars%kmxxc_fft - 1
             ph_wrk(ig) = stars%pgfft(ig) * gxc_fft(ig,idm)
          END DO

194
          DO js=1,jspins
195 196 197 198 199
             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

200
       IF (l_noco) THEN
201

202 203 204 205 206 207
          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) )
208

209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230
             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
231
             END DO
232 233 234

          END IF

235 236 237 238 239 240
       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)

241
       cqpw(:,:)= -den_pw(:,:jspins)
242 243 244 245 246 247 248 249 250
   
       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
             
251
             DO js=1,jspins
252 253 254 255 256 257 258 259
                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)

260
       IF (l_noco) THEN
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 305 306 307 308

          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
309
       END IF
310

311 312 313
       IF (PRESENT(xcpot)) THEN
          CALL xcpot%alloc_gradients(ifftxc3,jspins,grad)
       END IF
314 315

       ! 
316 317
       !     calculate the quantities such as abs(grad(rho)),.. used in
       !     evaluating the gradient contributions to potential and energy.
318
       ! 
319 320 321 322 323
       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)
324 325
          !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)
326 327 328 329 330

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

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


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

       IF (l_pw_w) THEN
375
          IF (dograds) THEN
376 377 378 379 380
             !----> 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
381 382 383 384 385 386 387
          !
          !----> Convolute with step function
          !
          DO i=0,ifftd-1
             vcon(i)=stars%ufft(i)*vcon(i)
          ENDDO
          bf3=0.0
388 389
          CALL fft3d(vcon(0),bf3, fg3, stars,-1)
          fg3=fg3*stars%nstr
390 391 392
          !
          !----> add to warped coulomb potential
          !
393 394 395 396 397
          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
398 399 400 401 402 403
       ENDIF
    END DO
  END SUBROUTINE pw_from_grid
    
  SUBROUTINE finish_pw_grid()
    IMPLICIT NONE
404
    IF (ALLOCATED(igxc_fft)) DEALLOCATE(igxc_fft,gxc_fft)
405 406 407
  END SUBROUTINE finish_pw_grid

END MODULE m_pw_tofrom_grid