pw_tofrom_grid.F90 11.5 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 68 69 70 71 72
    !.....------------------------------------------------------------------
    !------->          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
    !
    !-------------------------------------------------------------------
    USE m_grdrsis
    USE m_mkgxyz3
    USE m_fft3dxc
73
    USE m_fft3d
74
    USE m_types
75
    USE m_constants
76
    IMPLICIT NONE
77
    
78 79 80 81 82 83 84 85 86
    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(:,:)
87
  
88 89 90 91 92 93 94 95 96 97 98

    INTEGER      :: js,i,idm,ig,ndm,jdm
    REAL         :: rhotot
    !     .. Local Arrays ..
    COMPLEX, ALLOCATABLE :: cqpw(:,:),ph_wrk(:)
    REAL,    ALLOCATABLE :: bf3(:)
    REAL,    ALLOCATABLE :: rhd1(:,:,:),rhd2(:,:,:)
    REAL,    ALLOCATABLE :: mx(:),my(:)
    REAL,    ALLOCATABLE :: magmom(:),dmagmom(:,:),ddmagmom(:,:,:) 
    
    ! Allocate arrays
99
    ALLOCATE( bf3(0:ifftd-1))
100
    IF (dograds) THEN
Matthias Redies's avatar
Matthias Redies committed
101 102 103
       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) )
104
     ELSE
105
        IF (PRESENT(rho)) ALLOCATE(rho(0:ifftd-1,jspins))
106
     ENDIF
107
    IF (l_noco)  THEN
108
       IF (dograds) THEN
Daniel Wortmann's avatar
Daniel Wortmann committed
109 110 111 112 113
          ALLOCATE( mx(0:ifftxc3-1),my(0:ifftxc3-1),magmom(0:ifftxc3-1))
          ALLOCATE(dmagmom(0:ifftxc3-1,3),ddmagmom(0:ifftxc3-1,3,3) )
       ELSE
          ALLOCATE( mx(0:ifftd-1),my(0:ifftd-1),magmom(0:ifftd-1))
       ENDIF
114
    END IF
Matthias Redies's avatar
Matthias Redies committed
115

116 117 118
    IF (PRESENT(rho)) THEN
    !Put den_pw on grid and store into rho(:,1:2)
       DO js=1,jspins
119
          IF (dograds) THEN
120 121 122 123 124 125
             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
126

127 128
       IF (l_noco) THEN  
          !  Get mx,my on real space grid and recalculate rho and magmom
129
          IF (dograds) THEN
130 131 132 133 134
             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
135
          DO i=0,MIN(SIZE(rho,1),size(mx))-1 
136 137 138 139 140
             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 )
             rho(i,1)= rhotot+magmom(i)
             rho(i,2)= rhotot-magmom(i)
          END DO
141
       ENDIF
142
    ENDIF
143
    IF (dograds) THEN
144 145 146 147 148 149 150 151

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

152
       ALLOCATE(cqpw(stars%ng3,jspins))
153

154
       cqpw(:,:)= ImagUnit*den_pw(:,:jspins)
155 156 157 158 159 160
   
       DO idm=1,3
          DO ig = 0 , stars%kmxxc_fft - 1
             ph_wrk(ig) = stars%pgfft(ig) * gxc_fft(ig,idm)
          END DO

161
          DO js=1,jspins
162 163 164 165 166
             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

167
       IF (l_noco) THEN
168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183

          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
          END DO
       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)

184
       cqpw(:,:)= -den_pw(:,:jspins)
185 186 187 188 189 190 191 192 193
   
       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
             
194
             DO js=1,jspins
195 196 197 198 199 200 201 202
                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)

203
       IF (l_noco) THEN
204 205 206 207 208 209 210 211 212 213 214 215 216 217 218
          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 
       END IF
219

220 221 222
       IF (PRESENT(xcpot)) THEN
          CALL xcpot%alloc_gradients(ifftxc3,jspins,grad)
       END IF
223 224

       ! 
225 226
       !     calculate the quantities such as abs(grad(rho)),.. used in
       !     evaluating the gradient contributions to potential and energy.
227
       ! 
228 229 230 231 232
       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)
233 234
          !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)
235 236 237 238 239

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

240
          CALL mkgxyz3 (0*rhd1(0:,:,1),rhd1(0:,:,1),rhd1(0:,:,2),rhd1(0:,:,3),&
241 242
               rhd2(0:,:,1),rhd2(0:,:,3),rhd2(0:,:,6), rhd2(0:,:,5),rhd2(0:,:,4),rhd2(0:,:,2),grad)
       END IF
243 244
       
    ENDIF
245 246 247
    IF (PRESENT(rho)) THEN
       WHERE(ABS(rho) < d_15) rho = d_15
    ENDIF
248 249 250 251
   
  END SUBROUTINE pw_to_grid


252
  SUBROUTINE pw_from_grid(dograds,stars,l_pw_w,v_in,v_out_pw,v_out_pw_w)
253 254 255 256
    USE m_fft3d
    USE m_fft3dxc
    USE m_types
    IMPLICIT NONE
257
    LOGICAL,INTENT(IN)            :: dograds
258
    TYPE(t_stars),INTENT(IN)      :: stars
259
    REAL,INTENT(INOUT)            :: v_in(0:,:)
260
    LOGICAL,INTENT(in)            :: l_pw_w
261
    COMPLEX,INTENT(INOUT)         :: v_out_pw(:,:)
262
    COMPLEX,INTENT(INOUT),OPTIONAL:: v_out_pw_w(:,:)
263 264 265 266 267 268
    
    
    INTEGER              :: js,k,i
    REAL,ALLOCATABLE     :: bf3(:),vcon(:) 
    COMPLEX, ALLOCATABLE :: fg3(:)
    ALLOCATE( bf3(0:ifftd-1),fg3(stars%ng3))
269
    ALLOCATE ( vcon(0:ifftd-1) )
270 271
    DO js = 1,SIZE(v_in,2)
       bf3=0.0
272
       IF (dograds) THEN
273 274 275 276 277 278
          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
279
       DO k = 1,MERGE(stars%nxc3_fft,stars%ng3,dograds)
280
          v_out_pw(k,js) = v_out_pw(k,js) + fg3(k)
281 282 283
       ENDDO

       IF (l_pw_w) THEN
284
          IF (dograds) THEN
285 286 287 288 289
             !----> 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
290 291 292 293 294 295 296
          !
          !----> Convolute with step function
          !
          DO i=0,ifftd-1
             vcon(i)=stars%ufft(i)*vcon(i)
          ENDDO
          bf3=0.0
297 298
          CALL fft3d(vcon(0),bf3, fg3, stars,-1)
          fg3=fg3*stars%nstr
299 300 301 302
          !
          !----> add to warped coulomb potential
          !
          DO k = 1,stars%ng3
303
             v_out_pw_w(k,js) = v_out_pw_w(k,js) + fg3(k)
304 305 306 307 308 309 310
          ENDDO
       ENDIF
    END DO
  END SUBROUTINE pw_from_grid
    
  SUBROUTINE finish_pw_grid()
    IMPLICIT NONE
311
    IF (ALLOCATED(igxc_fft)) DEALLOCATE(igxc_fft,gxc_fft)
312 313 314
  END SUBROUTINE finish_pw_grid

END MODULE m_pw_tofrom_grid