LCOV - code coverage report
Current view: top level - vgen - pw_tofrom_grid.F90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 104 106 98.1 %
Date: 2019-09-08 04:53:50 Functions: 4 4 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------
       2             : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
       3             : ! This file is part of FLEUR and available as free software under the conditions
       4             : ! of the MIT license as expressed in the LICENSE file in more detail.
       5             : !--------------------------------------------------------------------------------
       6             : MODULE m_pw_tofrom_grid
       7             :   USE m_types
       8             :   PRIVATE
       9             :   REAL,PARAMETER:: d_15=1.e-15
      10             : 
      11             :   INTEGER :: ifftd,ifftxc3
      12             :   !----->  fft  information  for xc potential + energy
      13             :   INTEGER, ALLOCATABLE :: igxc_fft(:)
      14             :   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         170 :   SUBROUTINE init_pw_grid(xcpot,stars,sym,cell)
      19             :     USE m_prpxcfftmap
      20             :     USE m_types
      21             :     IMPLICIT NONE
      22             :     CLASS(t_xcpot),INTENT(IN)     :: xcpot
      23             :     TYPE(t_stars),INTENT(IN)      :: stars
      24             :     TYPE(t_sym),INTENT(IN)        :: sym
      25             :     TYPE(t_cell),INTENT(IN)       :: cell
      26             :     
      27             :     !---> set up pointer for backtransformation of from g-vector in
      28             :     !     positive domain of xc density fftbox into stars.
      29             :     !     also the x,y,z components of the g-vectors are set up to calculate
      30             :     !     derivatives.
      31             :     !     in principle this can also be done in main program once.
      32             :     !     it is done here to save memory.
      33             :     !
      34             :       
      35         170 :     ifftd=27*stars%mx1*stars%mx2*stars%mx3
      36         170 :     ifftxc3  = stars%kxc1_fft*stars%kxc2_fft*stars%kxc3_fft
      37         170 :     IF (xcpot%needs_grad()) THEN
      38         161 :        CALL prp_xcfft_map(stars,sym, cell, igxc_fft,gxc_fft)
      39             :     ENDIF
      40             :        
      41         170 :   END SUBROUTINE init_pw_grid
      42             :   
      43         380 :   SUBROUTINE pw_to_grid(xcpot,jspins,l_noco,stars,cell,den_pw,grad,rho)
      44             :     !.....------------------------------------------------------------------
      45             :     !------->          abbreviations
      46             :     !
      47             :     !     ph_wrk: work array containing phase * g_x,gy...... 
      48             :     !     den%pw: charge density stored as stars
      49             :     !     rho   : charge density stored in real space
      50             :     !     v_xc   : exchange-correlation potential in real space
      51             :     !     exc   : exchange-correlation energy density in real space
      52             :     !     kxc1d  : dimension of the charge density fft box in the pos. domain
      53             :     !     kxc2d  : defined in dimens.f program (subroutine apws).1,2,3 indic
      54             :     !     kxc3d  ; a_1, a_2, a_3 directions.
      55             :     !     kq(i) : i=1,2,3 actual length of the fft-box for which fft is done
      56             :     !     nstr  : number of members (arms) of reciprocal lattice (g) vector
      57             :     !             of each star.
      58             :     !     nxc3_fft: number of stars in the  charge density  fft-box
      59             :     !     ng3   : number of 3 dim. stars in the charge density sphere define
      60             :     !             by gmax
      61             :     !     kmxxc_fft: number of g-vectors forming the nxc3_fft stars in the
      62             :     !               charge density or xc-density sphere
      63             :     !     kimax : number of g-vectors forming the ng3 stars in the gmax-sphe
      64             :     !     igfft : pointer from the g-sphere (stored as stars) to fft-grid
      65             :     !             and     from fft-grid to g-sphere (stored as stars)
      66             :     !     pgfft : contains the phases of the g-vectors of sph.
      67             :     !     isn   : isn = +1, fft transform for g-space to r-space
      68             :     !             isn = -1, vice versa
      69             :     !
      70             :     !-------------------------------------------------------------------
      71             :     USE m_grdrsis
      72             :     USE m_mkgxyz3
      73             :     USE m_fft3dxc
      74             :     USE m_fft3d
      75             :     USE m_types
      76             :     use m_constants
      77             :     IMPLICIT NONE
      78             :     CLASS(t_xcpot),INTENT(IN)     :: xcpot
      79             :     INTEGER,INTENT(IN)            :: jspins
      80             :     LOGICAL,INTENT(IN)            :: l_noco
      81             :     TYPE(t_stars),INTENT(IN)      :: stars
      82             :     TYPE(t_cell),INTENT(IN)       :: cell
      83             :     COMPLEX,INTENT(IN)            :: den_pw(:,:)
      84             :     TYPE(t_gradients),INTENT(OUT) :: grad
      85             :     REAL,ALLOCATABLE,INTENT(out),OPTIONAL   :: rho(:,:)
      86             :   
      87             : 
      88             :     INTEGER      :: js,i,idm,ig,ndm,jdm
      89             :     REAL         :: rhotot
      90             :     !     .. Local Arrays ..
      91         380 :     COMPLEX, ALLOCATABLE :: cqpw(:,:),ph_wrk(:)
      92         190 :     REAL,    ALLOCATABLE :: bf3(:)
      93         190 :     REAL,    ALLOCATABLE :: rhd1(:,:,:),rhd2(:,:,:)
      94         190 :     REAL,    ALLOCATABLE :: mx(:),my(:)
      95         190 :     REAL,    ALLOCATABLE :: magmom(:),dmagmom(:,:),ddmagmom(:,:,:) 
      96             :     
      97             :     ! Allocate arrays
      98         190 :     ALLOCATE( bf3(0:ifftd-1))
      99         190 :     IF (xcpot%needs_grad()) THEN
     100         181 :        IF (PRESENT(rho)) ALLOCATE(rho(0:ifftxc3-1,jspins))
     101         181 :        ALLOCATE( ph_wrk(0:ifftxc3-1),rhd1(0:ifftxc3-1,jspins,3))
     102         181 :        ALLOCATE( rhd2(0:ifftxc3-1,jspins,6) )
     103             :      ELSE
     104           9 :         IF (PRESENT(rho)) ALLOCATE(rho(0:ifftd-1,jspins))
     105             :      ENDIF
     106         190 :     IF (l_noco)  THEN
     107         124 :        IF (xcpot%needs_grad()) THEN
     108         124 :           ALLOCATE( mx(0:ifftxc3-1),my(0:ifftxc3-1),magmom(0:ifftxc3-1))
     109         124 :           ALLOCATE(dmagmom(0:ifftxc3-1,3),ddmagmom(0:ifftxc3-1,3,3) )
     110             :        ELSE
     111           0 :           ALLOCATE( mx(0:ifftd-1),my(0:ifftd-1),magmom(0:ifftd-1))
     112             :        ENDIF
     113             :     END IF
     114             : 
     115         190 :     IF (PRESENT(rho)) THEN
     116             :     !Put den_pw on grid and store into rho(:,1:2)
     117         474 :        DO js=1,jspins
     118         474 :           IF (xcpot%needs_grad()) THEN
     119             :              CALL fft3dxc(rho(0:,js),bf3, den_pw(:,js), stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft,&
     120         295 :                   stars%nxc3_fft,stars%kmxxc_fft,+1, stars%igfft(0:,1),igxc_fft,stars%pgfft,stars%nstr)
     121             :           ELSE
     122           9 :              CALL fft3d(rho(0,js),bf3, den_pw(:,js), stars,+1)
     123             :           ENDIF
     124             :        END DO
     125             : 
     126         170 :        IF (l_noco) THEN  
     127             :           !  Get mx,my on real space grid and recalculate rho and magmom
     128         124 :           IF (xcpot%needs_grad()) THEN
     129             :              CALL fft3dxc(mx,my, den_pw(:,3), stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft,&
     130         124 :                   stars%nxc3_fft,stars%kmxxc_fft,+1, stars%igfft(0:,1),igxc_fft,stars%pgfft,stars%nstr)
     131             :           ELSE
     132           0 :              CALL fft3d(mx,my, den_pw(:,3), stars,+1)
     133             :           ENDIF
     134      614524 :           DO i=0,MIN(SIZE(rho,1),size(mx))-1 
     135      614400 :              rhotot= 0.5*( rho(i,1) + rho(i,2) )
     136      614400 :              magmom(i)= SQRT(  (0.5*(rho(i,1)-rho(i,2)))**2 + mx(i)**2 + my(i)**2 )
     137      614400 :              rho(i,1)= rhotot+magmom(i)
     138      614524 :              rho(i,2)= rhotot-magmom(i)
     139             :           END DO
     140             :        ENDIF
     141             :     ENDIF
     142         190 :     IF (xcpot%needs_grad()) THEN  
     143             : 
     144             :     ! In collinear calculations all derivatives are calculated in g-spce,
     145             :     ! in non-collinear calculations the derivatives of |m| are calculated in real space. 
     146             : 
     147             :     !-->   for d(rho)/d(x,y,z) = rhd1(:,:,idm) (idm=1,2,3).
     148             :     !
     149             :     !         ph_wrk: exp(i*(g_x,g_y,g_z)*tau) * g_(x,y,z).
     150             : 
     151         181 :        ALLOCATE(cqpw(stars%ng3,jspins))
     152             : 
     153         536 :        cqpw(:,:)= ImagUnit*den_pw(:,:jspins)
     154             :    
     155        1267 :        DO idm=1,3
     156     1313202 :           DO ig = 0 , stars%kmxxc_fft - 1
     157     1313202 :              ph_wrk(ig) = stars%pgfft(ig) * gxc_fft(ig,idm)
     158             :           END DO
     159             : 
     160        2854 :           DO js=1,jspins
     161             :              CALL fft3dxc(rhd1(0:,js,idm),bf3, cqpw(:,js), stars%kxc1_fft,stars%kxc2_fft,&
     162        1608 :                   stars%kxc3_fft,stars%nxc3_fft,stars%kmxxc_fft,+1, stars%igfft(0:,1),igxc_fft,ph_wrk,stars%nstr)
     163             :           END DO
     164             :        END DO
     165             : 
     166         181 :        IF (l_noco) THEN
     167             : 
     168         124 :           CALL grdrsis(magmom,cell,stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft,dmagmom )
     169             : 
     170      614524 :           DO i=0,ifftxc3-1
     171     1843324 :              DO idm=1,3
     172     1843200 :                 rhotot= rhd1(i,1,idm)/2.+rhd1(i,2,idm)/2.
     173     1843200 :                 rhd1(i,1,idm)= rhotot+dmagmom(i,idm) 
     174     2457600 :                 rhd1(i,2,idm)= rhotot-dmagmom(i,idm) 
     175             :              END DO
     176             :           END DO
     177             :        END IF
     178             : 
     179             :        !-->   for dd(rho)/d(xx,xy,yy,zx,yz,zz) = rhd2(:,:,idm) (idm=1,2,3,4,5,6)
     180             :        !
     181             :        !         ph_wrk: exp(i*(g_x,g_y,g_z)*tau) * g_(x,y,z) * g_(x,y,z)
     182             : 
     183         536 :        cqpw(:,:)= -den_pw(:,:jspins)
     184             :    
     185             :        ndm = 0
     186        1267 :        DO idm = 1,3
     187        2896 :           DO jdm = 1,idm
     188        1086 :              ndm = ndm + 1
     189     2626404 :              DO ig = 0 , stars%kmxxc_fft-1
     190     2626404 :                 ph_wrk(ig) = stars%pgfft(ig)*gxc_fft(ig,idm)*gxc_fft(ig,jdm)
     191             :              ENDDO
     192             :              
     193        5889 :              DO js=1,jspins
     194             :                 CALL fft3dxc(rhd2(0:,js,ndm),bf3, cqpw(:,js), stars%kxc1_fft,stars%kxc2_fft,&
     195        3216 :                      stars%kxc3_fft,stars%nxc3_fft,stars%kmxxc_fft,+1, stars%igfft(0:,1),igxc_fft,ph_wrk,stars%nstr)
     196             :              END DO
     197             :           END DO ! jdm 
     198             :        END DO   ! idm 
     199             : 
     200         181 :        DEALLOCATE(cqpw)
     201             : 
     202         181 :        IF (l_noco) THEN
     203         868 :           DO idm = 1,3
     204         496 :              CALL grdrsis(dmagmom(0,idm),cell,stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft,ddmagmom(0,1,idm) )
     205             :           END DO
     206             :           ndm= 0
     207         868 :           DO idm = 1,3
     208        1984 :              DO jdm = 1,idm
     209         744 :                 ndm = ndm + 1  
     210     3687516 :                 DO i=0,ifftxc3-1
     211     3686400 :                    rhotot= rhd2(i,1,ndm)/2.+rhd2(i,2,ndm)/2.
     212     3686400 :                    rhd2(i,1,ndm)= rhotot + ( ddmagmom(i,jdm,idm) + ddmagmom(i,idm,jdm) )/2. 
     213     3687144 :                    rhd2(i,2,ndm)= rhotot - ( ddmagmom(i,jdm,idm) + ddmagmom(i,idm,jdm) )/2. 
     214             :                 END DO
     215             :              ENDDO !jdm
     216             :           ENDDO   !idm 
     217             :        END IF
     218         181 :        CALL xcpot%alloc_gradients(ifftxc3,jspins,grad)
     219             :  
     220             :        !
     221             :        !     calculate the quantities such as abs(grad(rho)),.. used in
     222             :        !     evaluating the gradient contributions to potential and energy.
     223             :        !
     224         181 :        IF (PRESENT(rho)) THEN
     225             :           CALL mkgxyz3 (rho,rhd1(0:,:,1),rhd1(0:,:,2),rhd1(0:,:,3),&
     226         161 :                rhd2(0:,:,1),rhd2(0:,:,3),rhd2(0:,:,6), rhd2(0:,:,5),rhd2(0:,:,4),rhd2(0:,:,2),grad)
     227             :        ELSE
     228             :           !Dummy rho (only possible if grad is used for libxc mode)
     229             :           CALL mkgxyz3 (RESHAPE((/0.0/),(/1,1/)),rhd1(0:,:,1),rhd1(0:,:,2),rhd1(0:,:,3),&
     230          20 :                rhd2(0:,:,1),rhd2(0:,:,3),rhd2(0:,:,6), rhd2(0:,:,5),rhd2(0:,:,4),rhd2(0:,:,2),grad)
     231             :        END IF
     232             :        
     233             :     ENDIF
     234         190 :     IF (PRESENT(rho)) THEN
     235         170 :        WHERE(ABS(rho) < d_15) rho = d_15
     236             :     ENDIF
     237             :    
     238         190 :   END SUBROUTINE pw_to_grid
     239             : 
     240             : 
     241         530 :   SUBROUTINE pw_from_grid(xcpot,stars,l_pw_w,v_in,v_out_pw,v_out_pw_w)
     242             :     USE m_fft3d
     243             :     USE m_fft3dxc
     244             :     USE m_types
     245             :     IMPLICIT NONE
     246             :     CLASS(t_xcpot),INTENT(in)     :: xcpot
     247             :     TYPE(t_stars),INTENT(IN)      :: stars
     248             :     REAL,INTENT(INOUT)            :: v_in(0:,:)
     249             :     LOGICAL,INTENT(in)            :: l_pw_w
     250             :     COMPLEX,INTENT(INOUT)         :: v_out_pw(:,:)
     251             :     COMPLEX,INTENT(INOUT),OPTIONAL:: v_out_pw_w(:,:)
     252             :     
     253             :     
     254             :     INTEGER              :: js,k,i
     255        1060 :     REAL,ALLOCATABLE     :: bf3(:),vcon(:) 
     256         530 :     COMPLEX, ALLOCATABLE :: fg3(:)
     257         530 :     ALLOCATE( bf3(0:ifftd-1),fg3(stars%ng3))
     258         530 :     ALLOCATE ( vcon(0:ifftd-1) )
     259        1368 :     DO js = 1,SIZE(v_in,2)
     260    19898569 :        bf3=0.0
     261         838 :        IF (xcpot%needs_grad()) THEN
     262             :           CALL fft3dxc(v_in(0:,js),bf3, fg3, stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft,&
     263         811 :                stars%nxc3_fft,stars%kmxxc_fft,-1, stars%igfft(0:,1),igxc_fft,stars%pgfft,stars%nstr)
     264             :        ELSE
     265          27 :           vcon(0:)=v_in(0:,js)
     266          27 :           CALL fft3d(v_in(0:,js),bf3, fg3, stars,-1)
     267             :        ENDIF
     268      482439 :        DO k = 1,MERGE(stars%nxc3_fft,stars%ng3,xcpot%needs_grad())
     269      482439 :           v_out_pw(k,js) = v_out_pw(k,js) + fg3(k)
     270             :        ENDDO
     271             : 
     272        1368 :        IF (l_pw_w) THEN
     273         778 :           IF (xcpot%needs_grad()) THEN
     274             :              !----> Perform fft transform: v_xc(star) --> vxc(r) 
     275             :              !     !Use large fft mesh for convolution
     276         751 :              fg3(stars%nxc3_fft+1:)=0.0
     277         751 :              CALL fft3d(vcon(0),bf3, fg3, stars,+1)
     278             :           ENDIF
     279             :           !
     280             :           !----> Convolute with step function
     281             :           !
     282    19025329 :           DO i=0,ifftd-1
     283         778 :              vcon(i)=stars%ufft(i)*vcon(i)
     284             :           ENDDO
     285    19025329 :           bf3=0.0
     286         778 :           CALL fft3d(vcon(0),bf3, fg3, stars,-1)
     287         778 :           fg3=fg3*stars%nstr
     288             :           !
     289             :           !----> add to warped coulomb potential
     290             :           !
     291      776523 :           DO k = 1,stars%ng3
     292      776523 :              v_out_pw_w(k,js) = v_out_pw_w(k,js) + fg3(k)
     293             :           ENDDO
     294             :        ENDIF
     295             :     END DO
     296         530 :   END SUBROUTINE pw_from_grid
     297             :     
     298         170 :   SUBROUTINE finish_pw_grid()
     299             :     IMPLICIT NONE
     300         170 :     IF (ALLOCATED(igxc_fft)) DEALLOCATE(igxc_fft,gxc_fft)
     301         170 :   END SUBROUTINE finish_pw_grid
     302             : 
     303             : END MODULE m_pw_tofrom_grid

Generated by: LCOV version 1.13