LCOV - code coverage report
Current view: top level - init - prp_qfft.f90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 43 64 67.2 %
Date: 2019-09-08 04:53:50 Functions: 1 1 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             : 
       7             :       MODULE m_prpqfft
       8             :       use m_juDFT
       9             :       CONTAINS
      10          38 :       SUBROUTINE prp_qfft(&
      11             :      &                    stars,&
      12             :      &                    cell,noco,&
      13             :      &                    input)
      14             : !*********************************************************************
      15             : !     This subroutine prepares the necessary variables and checks
      16             : !     to calculate the plane wave chargedensity in the interstitial
      17             : !     (subroutine pwden.f) by fast fourier transform (FFT).
      18             : !     In certain cases rkmax is slightly reajusted to perform
      19             : !     a quick FFT based on powers   (2**P) * (3**Q) * (5**R) only.
      20             : !
      21             : !     dimensions kd(i)d for charge denisity box are checked.
      22             : !        s. bluegel, JRCAT, Feb. 97
      23             : !
      24             : !     subroutine boxdim added to subroutine to treat non-orthogonal
      25             : !     lattice systems.
      26             : !        s.bluegel, IFF, 18.Nov.97
      27             : !*********************************************************************
      28             : !
      29             :       USE m_ifft, ONLY : ifft235
      30             :       USE m_boxdim
      31             :       USE m_types
      32             :       IMPLICIT NONE
      33             : !
      34             : !---> Arguments
      35             : !
      36             :       TYPE(t_stars),INTENT(INOUT)  :: stars
      37             :       TYPE(t_cell),INTENT(IN)      :: cell
      38             :       TYPE(t_noco),INTENT(IN)      :: noco
      39             :       TYPE(t_input),INTENT(INOUT)  :: input  !<rkmax might be modified
      40             : !
      41             : !
      42             : !---> local variables
      43             : !
      44             :       INTEGER ksfft,mq1,mq2,mq3,istr,iofile,ng2_fft,kmxq2_fft
      45             :       REAL    arltv1,arltv2,arltv3,rknew
      46             : !
      47             : !---> intrinsic functions
      48             : !    
      49             :       INTRINSIC int,min
      50             : !
      51             : !---> data statement
      52             : !
      53             :       REAL gmaxp
      54             :       DATA gmaxp /2.0/
      55             : !
      56             : !------->          BACKGROUND       
      57             : !
      58             : !        Determine the limits  for the charge density fft-box
      59             : !        The charge density is the  square of the wavefunction.
      60             : !        Since the largest G-vector of the wavefunction is
      61             : !        given by |G + k| < rkmax, the largest G -vector gmaxq
      62             : !        contributing the charge-density is exactly :
      63             : !        gmaxq = gmaxp*rkmax,  with gmaxp =2. 
      64             : !        Therefore the box dimension must be  determined by :
      65             : !        m(i) >= gmaxp * rkmax*a(i)/(2*pi), where a(i) is the 
      66             : !        magnitude of the ith basis vector. (this is only true
      67             : !        for <a(i)|a(j)> = 0 for i .ne. j., otherwise see boxdim)
      68             : !        REMEMBER: gmax used through out the FLAPW code can be
      69             : !        larger than gmaxp*rkmax. 
      70             : !
      71             : !
      72             : !------->          ABBREVIATIONS
      73             : !
      74             : !   ksfft=(0,1) : KEY OF SELECTING FFT-PRDOGRAM AND RADIX-TYPE
      75             : !                      0  PROGRAM, RADIX-2 ONLY
      76             : !                      1  PROGRAM, RADIX-2, RADIX-3,RADIX-5
      77             : !   iofile      : device number for in and output
      78             : !   gmax        : actually used gmax troughtout the FLAPW code
      79             : !   gmaxq       : cut-off wavevector for charge density
      80             : !   rkmax       : cut-off for |g+k|
      81             : !   gmaxp       : gmaxp = gmaxq/rkmax, ideal: gmaxp=2
      82             : !   kq1,2,3d    : dimensions of the charge density fft-box
      83             : !                 ( in positive fft domain )
      84             : !   mq1,2,3     : actual size of the charge density fft-box
      85             : !                 ( in positive fft domain )
      86             : !   nq3_fft     : number of stars in the charge density sphere
      87             : !   kmxq_fft    : number of g-vectors in the charge density sphere
      88             : !   arltv(i)    : length of reciprical lattice vector along
      89             : !                 direction (i)
      90             : !
      91             : !---> Determine dimensions of fft-box of size mq1, mq2, mq3,
      92             : !     for which |G(mq1,mq2,mq3)| < Gmaxp*Kmax
      93             : !
      94             :       CALL boxdim(&
      95             :      &            cell%bmat,&
      96          38 :      &            arltv1,arltv2,arltv3)
      97             : !
      98          38 :       mq1 = int( gmaxp*input%rkmax/arltv1 ) + 1
      99          38 :       mq2 = int( gmaxp*input%rkmax/arltv2 ) + 1
     100          38 :       mq3 = int( gmaxp*input%rkmax/arltv3 ) + 1
     101             : 
     102             : !---> add + 1 in spin spiral calculation, to make sure that all G's are 
     103             : !---> still within the FFT-box after being shifted by the spin spiral
     104             : !---> q-vector.
     105          38 :       IF (noco%l_ss) THEN
     106           1 :          mq1 = mq1 + 1
     107           1 :          mq2 = mq2 + 1
     108           1 :          mq3 = mq3 + 1
     109             :       ENDIF         
     110             : !
     111             : !----> fft done in positive domain 
     112             : !      (Remember number of grid points not 2*m+1 because of cyclic 
     113             : !       boundary condition)
     114             : !
     115          38 :       mq1 = mq1 + mq1
     116          38 :       mq2 = mq2 + mq2
     117          38 :       mq3 = mq3 + mq3
     118             : !
     119             : !---> fft's are usually fastest for low primes
     120             : !     (restrict kqid to: kwid=  (2**P) * (3**Q) * (5**R)
     121             : !
     122          38 :       iofile = 6
     123          38 :       ksfft  = 1
     124          38 :       stars%kq1_fft = ifft235 (iofile,ksfft,mq1,gmaxp)
     125          38 :       stars%kq2_fft = ifft235 (iofile,ksfft,mq2,gmaxp)
     126          38 :       stars%kq3_fft = ifft235 (iofile,ksfft,mq3,gmaxp)
     127             : !+gb
     128             : !      kq1_fft = mq1
     129             : !      kq2_fft = mq2
     130             : !      kq3_fft = mq3
     131             : !
     132             : !---> For mq = 2**P, FFT is very fast. If mq very close to 2**P
     133             : !     choose this, even 2**P < mq . Therefore:
     134             : !     factor 0.5 added by S.B. 21.Aug.97 
     135             : !
     136             :       rknew = min( real( stars%kq1_fft )*arltv1, real( stars%kq2_fft )*arltv2, &
     137          38 :      &             real( stars%kq3_fft )*arltv3)
     138          38 :       rknew = 0.5 * rknew / gmaxp
     139          38 :       IF (rknew.LT.input%rkmax) THEN
     140           0 :          WRITE (6,'('' rkmax and true gmax recalculated '')')
     141           0 :          WRITE (6,2100) input%rkmax, rknew, rknew*rknew
     142           0 :          WRITE (6,2200) gmaxp*rknew, gmaxp*rknew*gmaxp*rknew
     143           0 :          input%rkmax = rknew
     144             :       ENDIF
     145             : !
     146             : !-----> gmax => gmaxp*rkmax 
     147             : !       (otherwise too few elements in arrays defined in strng1)
     148             : !
     149          38 :       IF (2*input%rkmax.GT.stars%gmax) THEN
     150           0 :          WRITE (6,'('' gmax must be at least 2*rkmax'')')
     151           0 :          WRITE (6,'('' increase gmax , or reduce rkmax'')')
     152           0 :          WRITE (6,'('' rkmax ='',f10.3,''  gmax ='',f10.3)') input%rkmax,stars%gmax
     153             :          CALL juDFT_error("rkmax,gmax",calledby ="prp_qfft",hint&
     154           0 :      &        ="gmax must be at least 2*rkmax")
     155             :       ENDIF
     156             : !
     157             : !     mq1 = kq1_fft
     158             : !     mq2 = kq2_fft
     159             : !     mq3 = kq3_fft
     160             : !
     161             : 
     162             : !
     163             : !-----> how many stars are in charge density sphere?
     164             : !       assume stars are ordered according to length
     165             : !
     166             : !---> 3d stars
     167          38 :       stars%ng3_fft = 0
     168      106903 :       DO istr = 1 , stars%ng3
     169      106903 :          IF ( stars%sk3(istr).LE.gmaxp*input%rkmax ) THEN
     170       21630 :             stars%ng3_fft = istr
     171             :          ENDIF
     172             :       ENDDO
     173             : !---> 2d stars
     174          38 :       ng2_fft = 0
     175        3822 :       DO istr = 1 , stars%ng2
     176        3822 :          IF ( stars%sk2(istr).LE.gmaxp*input%rkmax ) THEN
     177        1212 :             ng2_fft = istr
     178             :          ENDIF
     179             :       ENDDO
     180             : !
     181          38 :       IF ( stars%ng3_fft.EQ.0 ) THEN
     182           0 :          WRITE(6,'('' presumably ng3 too small '')')
     183             :          WRITE(6,'('' sk3max, gmaxp*rkmax '', 2f10.6)') &
     184           0 :      &                stars%sk3(stars%ng3),gmaxp*input%rkmax
     185           0 :          CALL juDFT_error("presumably ng3 too small","prp_qfft")
     186             :       ENDIF
     187             : !
     188          38 :       IF ( stars%ng3_fft.GT.stars%ng3 ) THEN
     189           0 :          WRITE(6,'('' nq3_fft > n3d '')')
     190           0 :          WRITE(6,'('' nq3_fft, n3d '',2i10)') stars%ng3_fft, stars%ng3
     191           0 :           CALL juDFT_error("nq3_fft > n3d",calledby="prp_qfft")
     192             :       ENDIF
     193             : !
     194             : !-----> check that all nq3_fft stars fit into the charge density FFT-box
     195             : !
     196       21668 :       DO istr = 1 , stars%ng3_fft
     197             :         IF ( ( 2*stars%kv3(1,istr).GT.stars%kq1_fft ) .OR.&
     198       21630 :      &       ( 2*stars%kv3(2,istr).GT.stars%kq2_fft ) .OR.&
     199          38 :      &       ( 2*stars%kv3(3,istr).GT.stars%kq3_fft ) ) THEN
     200           0 :           WRITE (6,'('' not all nq3_fft stars in chg. den. FFT box'')')
     201           0 :           WRITE (6,'('' inconsistency in def.s see also strgn1'')')
     202             :           WRITE (6,'('' mq1d,mq2d,mq3d,kv1,kv2,kv3 '',6i5)')&
     203           0 :      &                  stars%kq1_fft,stars%kq2_fft,stars%kq3_fft,2*stars%kv3(1,istr),2*stars%kv3(2,istr),&
     204           0 :      &                                 2*stars%kv3(3,istr)
     205             :           CALL juDFT_error("not all nq3_fft stars in chg. den. FFT box",&
     206           0 :      &         calledby ="prp_qfft")
     207             :         ENDIF
     208             :       ENDDO
     209             : !
     210             : !-----> how many G-vectors belong to these nq3_fft stars                    
     211             : !
     212             : !---> 3d vectors
     213          38 :       stars%kmxq_fft = 0
     214       21668 :       DO istr = 1 , stars%ng3_fft
     215       21668 :          stars%kmxq_fft = stars%kmxq_fft + stars%nstr(istr)
     216             :       ENDDO
     217          38 :       IF ( stars%kmxq_fft .GT. stars%kq1_fft*stars%kq2_fft*stars%kq3_fft ) THEN
     218             :          WRITE (6,'('' array dimensions in later subroutines too'',&
     219           0 :      &             '' small'',2i10)') stars%kmxq_fft,stars%kq1_fft*stars%kq2_fft*stars%kq3_fft
     220             :       ENDIF
     221             : !---> 2d vectors
     222          38 :       kmxq2_fft = 0
     223        1250 :       DO istr = 1 , ng2_fft
     224        1250 :          kmxq2_fft = kmxq2_fft + stars%nstr2(istr)
     225             :       ENDDO
     226          38 :       IF ( kmxq2_fft .GT. stars%kq1_fft*stars%kq2_fft ) THEN
     227             :          WRITE (6,'('' array dimensions in later subroutines too'',&
     228           0 :      &             '' small'',2i10)') kmxq2_fft,stars%kq1_fft*stars%kq2_fft
     229             :       ENDIF
     230             : 
     231             : !
     232             :  2100 FORMAT (/,1x,'old rkmax   =',f10.5, '(a.u.)**(-1) ==>  new rkmax '&
     233             :      &      ,'  =',f10.5,'(a.u.)**(-1) ',/,t38,'==>  new E_cut(wf)   =',&
     234             :      &            f10.5,' Ry')
     235             :  2200 FORMAT (/,1x,'true gmax   =',f10.5, '(a.u.)**(-1)',&
     236             :      &       /,t38,'==>  new E_cut(chg)  =', f10.5,' Ry')
     237             : 
     238          38 :       END SUBROUTINE prp_qfft
     239             :       END MODULE m_prpqfft

Generated by: LCOV version 1.13