LCOV - code coverage report
Current view: top level - init - prp_xcfft.f90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 41 60 68.3 %
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_prpxcfft
       8             :       use m_juDFT
       9             :       CONTAINS
      10          38 :       SUBROUTINE prp_xcfft(&
      11             :      &                     stars,&
      12             :      &                     input,cell,&
      13             :      &                     xcpot)
      14             : !*********************************************************************
      15             : !     this subroutine prepares the necessary variables and checks
      16             : !     to calculate the xc-potential and energy in the interstitial
      17             : !     (subroutine visxc(g).f) by fast fourier transform (fft).
      18             : !     in certain cases gmaxxc 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, iff , oct. 97
      23             : !     subroutine boxdim added to subroutine to treat non-othogonal
      24             : !     lattice systems
      25             : !        s.bluegel, IFF, 18.Nov.97
      26             : !*********************************************************************
      27             : !
      28             :       USE m_ifft, ONLY : ifft235
      29             :       USE m_types
      30             :       USE m_boxdim
      31             :       IMPLICIT NONE
      32             :       TYPE(t_stars),INTENT(INOUT)   :: stars
      33             :       TYPE(t_input),INTENT(IN)      :: input
      34             :       TYPE(t_cell),INTENT(IN)       :: cell
      35             :       CLASS(t_xcpot),INTENT(INOUT)  :: xcpot
      36             : 
      37             : !
      38             : !---> local variables
      39             : !
      40             :       INTEGER ksfft,mxc1,mxc2,mxc3,istr,iofile
      41             :       REAL    arltv1,arltv2,arltv3,gmxxc_new
      42             : !
      43             : !---> intrinsic functions
      44             : !
      45             :       INTRINSIC int,min
      46             : !
      47             : !------->          background
      48             : !
      49             : !        determine the limits  for the xc energy/pot fft-box.
      50             : !        since the xc functional has a non-linear dependence 
      51             : !        of the charge density in many cases a gmax to determine
      52             : !        the xc-potential or energy of gmax=2*rkmax=gmaxq is 
      53             : !        in many cases not sufficent. Therefore, we use a larger
      54             : !        gmax called gmaxxc to evaluate the xc-potential and energy
      55             : !        in the interstitial region. This leads to a larger xc-fft-box.
      56             : !        the box dimension is  determined by :
      57             : !        m(i) >= gmaxxc*a(i)/(2*pi), where a(i) is the
      58             : !        magnitude of the ith basis vector. (this is only true
      59             : !        for <a(i)|a(j)> = 0 for i .ne. j., otherwise see boxdim)
      60             : !        remember: gmax used through out the flapw code must be
      61             : !        larger than or equal  gmax >= gmaxxc >= gmaxp*rkmax.
      62             : !
      63             : !
      64             : !------->          abbreviations
      65             : !
      66             : !   ksfft=(0,1) : key of selecting fft-prdogram and radix-type
      67             : !                      0  program, radix-2 only
      68             : !                      1  program, radix-2, radix-3,radix-5
      69             : !   iofile      : device number for in and output
      70             : !   rkmax       : cut-off for |g+k|
      71             : !   gmax        : actually used largest g-vector used througout 
      72             : !               : the calculation
      73             : !   gmaxxc      : cut-off wavevector for fast fourier transform of 
      74             : !                 xc-potential and energy
      75             : !                 gmaxxc usually gmaxxc >= 2*rkmax
      76             : !   kxc1,2,3d   : dimensions of the xc-pot/energy fft-box
      77             : !                 ( in positive fft domain )
      78             : !   kxc1,2,3_fft: actual size of the xc-pot/energy fft-box
      79             : !                 ( in positive fft domain )
      80             : !   nxc3_fft     : number of stars in the xc-pot/energy sphere
      81             : !   kmxxc_fft    : number of g-vectors in the xc-pot/energy sphere
      82             : !   arltv(i)    : length of reciprical lattice vector along
      83             : !                 direction (i)
      84             : !
      85          38 :       IF (.NOT.xcpot%needs_grad()) xcpot%gmaxxc=stars%gmax
      86          38 :       WRITE (6,'('' gmaxxc should be: 2*kmax <= gmaxxc <= gmax '')')
      87          38 :       IF ( abs( xcpot%gmaxxc - stars%gmax ) .le. 10.0**(-6) ) THEN
      88             :         WRITE (6,'('' concerning memory, you may want to choose'',&
      89          10 :      &              '' a smaller value for gmax'')')
      90             :       END IF
      91          38 :       IF ( xcpot%gmaxxc .LE. 10.0**(-6) ) THEN
      92             :          WRITE (6,'(" gmaxxc=0 : gmaxxc=gmax choosen as default",&
      93           0 :      &              " value")')
      94             :          WRITE (6,'(" concerning memory, you may want to choose",&
      95           0 :      &              " a smaller value for gmax")')
      96           0 :          xcpot%gmaxxc=stars%gmax
      97             :       END IF
      98          38 :       IF ( xcpot%gmaxxc .LE. 2*input%rkmax ) THEN
      99             :          WRITE (6,'('' concerning accuracy and total energy'',&
     100             :      &              '' convergence, you may want'',/,&
     101           0 :      &              '' to choose a larger gmaxxc '')')
     102             :       END IF
     103          38 :       write (6,'('' gmaxxc ='',f10.6)') xcpot%gmaxxc
     104             : !
     105             : !---> Determine dimensions of fft-box of size mxc1, mxc2, mxc3,
     106             : !     for which |G(mxc1,mxc2,mxc3)| < Gmaxxc
     107             : !
     108             :       CALL boxdim(&
     109             :      &            cell%bmat,&
     110          38 :      &            arltv1,arltv2,arltv3)!keep
     111             : !
     112          38 :       mxc1 = int( xcpot%gmaxxc/arltv1 ) + 1
     113          38 :       mxc2 = int( xcpot%gmaxxc/arltv2 ) + 1
     114          38 :       mxc3 = int( xcpot%gmaxxc/arltv3 ) + 1
     115             : !
     116             : !----> fft done in positive domain
     117             : !      (remember number of grid points not 2*m+1 because of cyclic
     118             : !       boundary condition)
     119             : !
     120          38 :       mxc1 = mxc1 + mxc1
     121          38 :       mxc2 = mxc2 + mxc2
     122          38 :       mxc3 = mxc3 + mxc3
     123             : !
     124             : !---> fft's are usually fastest for low primes
     125             : !     (restrict kqid to: kwid=  (2**p) * (3**q) * (5**r)
     126             : !
     127          38 :       iofile = 6
     128          38 :       ksfft  = 1
     129          38 :       stars%kxc1_fft = ifft235 (iofile,ksfft,mxc1,2.0)
     130          38 :       stars%kxc2_fft = ifft235 (iofile,ksfft,mxc2,2.0)
     131          38 :       stars%kxc3_fft = ifft235 (iofile,ksfft,mxc3,2.0)
     132             : !
     133             : !---> for mxc = 2**p, fft is very fast. if mq very close to 2**p
     134             : !     choose this, even 2**p < mxc . therefore:
     135             : !
     136             :       gmxxc_new = min( real( stars%kxc1_fft )*arltv1, &
     137          38 :      &                 real( stars%kxc2_fft )*arltv2, real( stars%kxc3_fft )*arltv3)
     138          38 :       gmxxc_new = 0.5 * gmxxc_new
     139             : 
     140          38 :       IF (gmxxc_new.LT.xcpot%gmaxxc) THEN
     141           0 :          WRITE (6,'('' gmaxxc recalculated '')')
     142           0 :          WRITE (6,2100) xcpot%gmaxxc, gmxxc_new, gmxxc_new*gmxxc_new
     143           0 :          xcpot%gmaxxc = gmxxc_new
     144             :       ENDIF
     145             : !
     146             : !-----> gmax => gmaxxc
     147             : !       (otherwise too few elements in arrays defined in strng1)
     148             : !
     149          38 :       IF (gmxxc_new.GT.stars%gmax) THEN
     150          11 :          WRITE (6,'('' gmax must be at least gmxxc_new'')')
     151          11 :          WRITE (6,'('' increase gmax , or reduce gmaxxc'')')
     152             :          WRITE (6,'('' gmxxc_new ='',f10.3,''  gmax ='',f10.3)') &
     153          11 :      &                gmxxc_new, stars%gmax
     154             : !cc          CALL juDFT_error("gmxxc_new.gt.gmax",calledby="prp_xcfft")
     155             :       ENDIF
     156             : !
     157             : !------> check dimensions of fft chargedensity box used in pwden.f
     158             : !
     159             :        IF (stars%kxc1_fft.GT.stars%kxc1_fft .OR. stars%kxc2_fft.gt.stars%kxc2_fft .OR. &
     160             :      &                            stars%kxc3_fft.gt.stars%kxc3_fft) THEN
     161             :           WRITE (6,'('' box dim. for fft too small'')')
     162             :           WRITE (6,2110) stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft,stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft
     163             :           CALL juDFT_error("mxc[1,2,3]d>kxc[1,2,3]d ",calledby&
     164             :      &         ="prp_xcfft")
     165             :        ENDIF
     166             :  2110  FORMAT (' stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft,stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft ',6i5)
     167             : !
     168             : !-----> how many stars are in charge density sphere?
     169             : !       assume stars are ordered according to length
     170             : !
     171          38 :       stars%nxc3_fft = 0
     172      106903 :       DO istr = 1 , stars%ng3
     173      106903 :          IF ( stars%sk3(istr).LE.xcpot%gmaxxc ) THEN
     174       65447 :             stars%nxc3_fft = istr
     175             :          ENDIf
     176             :       ENDDO
     177             : !
     178          38 :       IF ( stars%nxc3_fft.EQ.0 ) THEN
     179           0 :          WRITE (6,'('' presumably ng3 too small '')')
     180             :          WRITE (6,'('' sk3max, gmaxxc '', 2f10.6)')&
     181           0 :      &                stars%sk3(stars%ng3),xcpot%gmaxxc
     182           0 :          CALL juDFT_error("nxc3_fft==0",calledby ="prp_xcfft")
     183             :       ENDIF
     184             : !
     185          38 :       IF ( stars%nxc3_fft.GT.stars%ng3 ) THEN
     186           0 :          WRITE (6,'('' nxc3_fft > n3d '')')
     187           0 :          WRITE (6,'('' nxc3_fft, n3d '',2i10)') stars%nxc3_fft, stars%ng3
     188           0 :          CALL juDFT_error("nxc3_fft>n3d ",calledby ="prp_xcfft")
     189             :       ENDIF
     190             : !
     191             : !-----> check that all nxc3_fft stars fit into the xc-pot/energy fft-box
     192             : !
     193       65485 :       DO istr = 1 , stars%nxc3_fft
     194             :         IF ( ( 2*stars%kv3(1,istr).gt.stars%kxc1_fft ) .OR.&
     195       65447 :      &       ( 2*stars%kv3(2,istr).gt.stars%kxc2_fft ) .OR.&
     196          38 :      &       ( 2*stars%kv3(3,istr).gt.stars%kxc3_fft ) ) THEN
     197           0 :           WRITE(6,'('' not all nxc3_fft stars in xc-pot/eng fft box'')')
     198           0 :           WRITE(6,'('' inconsistency in def.s see also strgn1'')')
     199             :           WRITE(6,'('' kxc1_fft,kxc2_fft,kxc3_fft,kv1,kv2,kv3 '',6i5)')&
     200           0 :      &                 stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft,2*stars%kv3(1,istr),&
     201           0 :      &                 2*stars%kv3(2,istr),2*stars%kv3(3,istr)
     202             :           CALL juDFT_error("2*stars([1,2,3],istr)>mxc[1,2,3]d",calledby&
     203           0 :      &         ="prp_xcfft")
     204             :         ENDIF
     205             :       ENDDO
     206             : !
     207             : !-----> how many g-vectors belong to these nxc3_fft stars
     208             : !
     209          38 :       stars%kmxxc_fft = 0
     210       65485 :       DO istr = 1 , stars%nxc3_fft
     211       65485 :          stars%kmxxc_fft = stars%kmxxc_fft + stars%nstr(istr)
     212             :       ENDDO
     213             : 
     214          38 :       IF ( stars%kmxxc_fft .gt. stars%kxc1_fft*stars%kxc2_fft*stars%kxc3_fft ) then
     215             :          WRITE (6,'('' array dimensions in later subroutines too'',&
     216           0 :      &             '' small'')')
     217             :       ENDIF
     218             : 
     219             :  2100 format(/,1x,'old gmaxxc  =',f10.5, '(a.u.)**(-1) ==>  new gmaxxc'&
     220             :      &      ,'  =',f10.5,'(a.u.)**(-1) ',/,t38,'==>  new e_cut(xc)   =',&
     221             :      &            f10.5,' ry')
     222             : 
     223          38 :       END SUBROUTINE prp_xcfft
     224             :       END MODULE m_prpxcfft

Generated by: LCOV version 1.13